Yosra Ellili
email: yousraellili91@gmail.com

Brendan Philip Malone

Didier Michot

Budiman Minasny

Sébastien Vincent

Christian Walter

Blandine Lemercier

Comparing three approaches of spatial disaggregation of legacy soil maps based on 1 DSMART algorithm

Keywords: digital soil mapping, soil landscape relationships, spatial disaggregation, DSMART

1) Introduction

Characterizing soil variability especially over large areas, remains a crucial challenge to foster sustainable management of agronomic and environmental issues and help stakeholders to design regional projects [START_REF] Chaney | POLARIS: A 30-meter probabilistic soil series map of the contiguous United States[END_REF]. At the regional as well as country level, soil maps are often available at coarse spatial resolution [START_REF] Bui | Disaggregation of polygons of surficial geology and soil maps using spatial modelling and legacy data[END_REF] which limits their ability to depict accurate soil information. For instance, the finest soils maps covering France were elaborated by administrative region at 1:250,000 scale, via a set of polygons, called Soil Map Units (SMU) with crisp boundaries. The delineation of SMU is based on soil survey programmes involving pedologists' expertise. In a coarse scale map, each polygon includes one or several Soil Type Unit (STU), which are not explicitly mapped, but their proportions and their environmental conditions, as well as soil characteristics, are provided in a detailed database (Le [START_REF] Bris | Organisation des sols d'Ille-et-Vilaine. Version 1.1[END_REF].

To improve soil variability knowledge and overcome the limitation of a coarse mapping scale, several methods have emerged in the field of Digital Soil Mapping (DSM). These methods offer useful tools to predict soil spatial pattern from scarce or limited soil datasets by exploiting the availability of model based methods and an extensive array of spatialise (and more often than not gridded) environmental variables. In recent decades, DSM techniques have been increasingly used to downscale soil information and improve their spatial resolution. Depending on the quality of data and the complexity of soil cover, [START_REF] Minasny | Methodologies for Global Soil Mapping[END_REF] supply a workflow that outlines different models that can be explored. In general, two main pathways can be distinguished: point based DSM approaches and map disaggregation approaches (Odgers et al., 2014;Holmes et al., 2015). Point DSM approaches used legacy soil profiles, which are irregularly distributed and collected according to specific objectives rather than to optimise a statistical criterion (Holmes et al., 2015). The spatial distribution of soil properties can be estimated by fitting geostatistical models such as ordinary kriging (Odgers et al., 2014;Holmes et al., 2015;[START_REF] Chaney | POLARIS: A 30-meter probabilistic soil series map of the contiguous United States[END_REF][START_REF] Vincent | Spatial disaggregation of complex Soil Map Units at the regional scale based on soil-landscape relationships[END_REF]Chen et al., 2018) or cokriging, which takes into account the spatial interrelations among several soil properties (Webster and Oliver, 2007). Additionally, McBratney et al. (2003) developed the SCORPAN soil landscape model. It is an empirical quantitative function of environmental covariates, allowing predicting soil attributes (soil type or soil property) based on correlative and statistical relationships with predictor variables.

The second approach, known as spatial disaggregation, attempts to downscale the soil map unit information to delineate unmapped STUs [START_REF] Bui | Disaggregation of polygons of surficial geology and soil maps using spatial modelling and legacy data[END_REF]Odgers et al., 2014;Holmes et al., 2015). Alternatively, it can be defined as the process that allows estimating soil properties at a finer scale than the initial soil map. Several techniques have been demonstrated through soil science literature and tested in different case studies around the world. For instance, [START_REF] Kempen | Updating the 1:50,000 Dutch soil map using legacy soil data: A multinomial logistic regression approach[END_REF] have explored the use of multinomial logistic regression (MLR) for digital soil mapping. Other techniques have also been applied as decision trees using rule based induction [START_REF] Bui | Disaggregation of polygons of surficial geology and soil maps using spatial modelling and legacy data[END_REF], Bayesian techniques (Bui et al., 1999) and an area to point kriging method [START_REF] Kerry | Disaggregation of legacy soil data using area to point kriging for mapping soil organic carbon at the regional scale[END_REF].

In the DSM field, machine learning techniques are increasingly used to elucidate the spatial distribution of both soil type and soil properties across a large range of scale [START_REF] Bui | Disaggregation of polygons of surficial geology and soil maps using spatial modelling and legacy data[END_REF][START_REF] Scull | The application of classification tree analysis to soil type prediction in a desert landscape[END_REF]Lacoste et al., 2011;[START_REF] Lemercier | Extrapolation at regional scale of local soil knowledge using boosted classification trees: A two-step approach[END_REF]Nauman and Thompson, 2014;Holmes et al., 2015;[START_REF] Vaysse | Evaluating Digital Soil Mapping approaches for mapping GlobalSoilMap soil properties from legacy data in Languedoc-Roussillon (France)[END_REF][START_REF] Ellili | Mapping soil organic carbon stock change by soil monitoring and digital soil mapping at the landscape scale[END_REF]. They were also applied to disaggregate superficial geology maps available at 1: 250 000 scale in Australia [START_REF] Bui | Disaggregation of polygons of surficial geology and soil maps using spatial modelling and legacy data[END_REF]. The main advantage of these approaches is they allow handling both quantitative and categorical (ordinal or nominal) soil and environmental variables, as explanatory covariates [START_REF] Bui | Disaggregation of polygons of surficial geology and soil maps using spatial modelling and legacy data[END_REF]. Odgers et al. (2014) have developed a machine learning algorithm entitled Disaggregation and Harmonisation of Soil Map Units Through Resampled Classification Trees (DSMART) to predict STU as a function of the high resolution environmental data supplied over different study areas in Australia. The DSMART algorithm is based on a calibration dataset derived from a random selection of a fixed number of sampling points within each soil polygon. Each sampling point is then assigned to one soil type following a weighted random allocation procedure based on the proportions informed in the soil map database. The same procedure was applied by [START_REF] Chaney | POLARIS: A 30-meter probabilistic soil series map of the contiguous United States[END_REF] to spatially disaggregate the soil map of the contiguous United States at a 30m spatial resolution using petascale High Performance Computer (HPC). Because integration of pedological knowledge has been recognized as an effective way to improve digital soil mapping approaches [START_REF] Cook | Use of airborne gamma radiometric data for soil mapping[END_REF][START_REF] Walter | Integrating pedological knowledge into digital soil mapping[END_REF]Stoorvogel et al., 2017;[START_REF] Machado | Spatial Disaggregation of Multi-Component Soil Map Units Using Legacy Data and a Tree-Based Algorithm in Southern Brazil[END_REF][START_REF] Møller | Improved disaggregation of conventional soil maps[END_REF], [START_REF] Vincent | Spatial disaggregation of complex Soil Map Units at the regional scale based on soil-landscape relationships[END_REF] have applied the DSMART algorithm with additional expert soil landscape rules describing soil distribution in the local context of the Brittany region (France). By adding supplement sampling points to the calibration dataset selected according to soil parent material, soil redoximorphic conditions and topographic features, and by integrating soil landscape relationships in the DSMART sample allocation scheme, the authors obtained a coherent soil spatial distribution observing soil organisation along hillslopes and occurrence of intensely waterlogged soils in the stream neighbourhood, as observed in Brittany.

This study aimed to test the hypothesis that adding soil landscape relationships in the disaggregation procedure improved the accuracy of produced disaggregated soil maps. This involves assessing the contribution of soil landscape relationships implemented in the DSMART algorithm by [START_REF] Vincent | Spatial disaggregation of complex Soil Map Units at the regional scale based on soil-landscape relationships[END_REF]. To achieve this objective, we compared disaggregated soil maps either derived from the original DSMART algorithm, the DSMART algorithm with extra soil observations and the DSMART algorithm fed by soil landscape relationships over an area of 6,775 km² in the eastern part of Brittany, France.

2) Materials and methods

2.1) Study area

The Ille et Vilaine department covers an area of 6,775 km² and is located at the eastern part of 2.2.1) Regional soil database at 1:250 000 scale In Brittany, soils are represented through a regional geographic database called "Référentiel Régional Pédologique (RRP)" available at 1:250,000 scale (INRA Infosol, 2014).This regional database identifies soils within Soil Map Units (SMUs), each containing one to several soil types called Soil Type Units (STUs). STUs are defined as areas with homogeneous soil forming factors, such as morphology, geology, and climate. In the study area, 96 SMU and 171 STU have been distinguished and represented by a spatial coverage of 479 polygons. In the regional database, SMUs were spatially delimited with crisp boundaries, while STUs were not explicitly mapped, but their proportion in each SMU as well as associated environmental and soil characteristics were accurately described in a semantic database (Le Bris et al., 2013;INRA Infosol, 2014).

2.2.2) Soil validation data

To assess the quality of disaggregated soil maps, three validation datasets were used (Fig. 1):

 135 soil profiles chosen following a stratified random sampling design and specifically described and sampled from March to May 2017 for independent validation purposes in the framework of the Soilserv research project (Ellili et al., 2019, submitted).

 755 legacy soil profiles collected between 2005 and 2008 during the "Sols de Bretagne" programme (INRA Infosol, 2014).These profiles were sampled to characterize hydromorphic soil conditions and soil landscape heterogeneity.

Existing detailed soil maps (1:25,000) covering 87,150 ha, surveyed according to [START_REF] Rivière | Méthode Tarière Massif Armoricain[END_REF] and revised later to adapt to the STU typologies developed in the RRP (Le [START_REF] Bris | Organisation des sols d'Ille-et-Vilaine. Version 1.1[END_REF].

All soil profiles were allocated after description and analysis by an expert to a suitable STU. Both legacy soil profiles and detailed maps were converted to raster format to perfectly meet the prediction raster at 50m spatial resolution.

2.3) Environmental covariates

The SCORPAN concept (McBratney et al., 2013) allows one to predict STU as a function of a set of covariates describing seven soil forming factors, namely soil properties (s), climate (c), organisms (o), relief (r), parent material (p), age (a) and geographic position (n). In this study, ten environmental variables (Table 1) were considered as covariates in the disaggregation process at a 50m spatial resolution. Terrain attributes included elevation, slope, Compound Topographic Index (CTI) (Beven andKirkby, 1979, Merot et al., 1995) and Topographic Position Index (TPI) [START_REF] Vincent | Spatial disaggregation of complex Soil Map Units at the regional scale based on soil-landscape relationships[END_REF]) that together were derived from a 50m resolution Digital Elevation Model [START_REF] Ign | World reference base for soil resources 2006, first update 2007[END_REF]. These attributes were computed using ArcGIS 10.1 (ESRI, 2002) and MNT surf software [START_REF] Squividant | MNTSurf: Logiciel de traitement des modèles numériques de terrain[END_REF].

Environmental attributes describing soil parent material (Lacoste et al., 2011) and hydromorphic soil conditions via waterlogging index [START_REF] Lemercier | Apport de la cartographie numérique des sols pour prédire l'hydromorphie et l'extension des zones humides potentielles à l'échelle régionale[END_REF] were obtained using decision tree methods. Waterlogging index derives from a natural soil drainage prediction. Four classes were distinguished: well drained, moderately drained, poorly drained and very poorly drained. Aeolian silt deposits and Soil Map Units boundaries are environmental covariates also obtained via expert knowledge from soil scientists.

Landscape units reflecting vegetation, land use, and relief attributes were derived from a MODIS imagery by supervised classification (Le Du Bayo et al., 2008). The Airborne gamma ray spectrometry variable (K:Th ratio) [START_REF] Messner | Apport de la Spectrométrie Gamma Aéroportée pour la cartographie numérique des sols[END_REF], characterizing the degree of weathering of the geological material, was also taken into account.

All soil environmental covariates were converted to raster format at 50 m spatial resolution. The open source DSMART algorithm (Odgers et al., 2014) was applied to spatially disaggregate the existing legacy soil map at 1:250,000 scale. DSMART algorithm uses machine learning classification trees implemented in C5.0 [START_REF] Quinlan | C4.5: Programs for Machine Learning[END_REF] to build a decision tree from a target variable (STU) and the environmental covariates supplied. The DSMART algorithm was written in the Python programming language by Odgers et al. (2014) and was recently translated in the R programming language.

Running DSMART algorithm requires four main steps (Fig. 2): 1) Polygon sampling by a random selection of a fixed number of sampling points (n=30) within each polygon. This procedure allowed to select a total of 14,370 sampling points, per iteration, covering the study area and ensured that all polygons were sampled.

2) Soil Type Unit (STU) assignment to each sampling point following a weighted random allocation method. This step was based on the proportion of each STU informed in the RRP database.

3) Decision tree generation: the full set of sampling points were spatially intersected with the selected environmental covariates. This georeferenced dataset was then used as a calibration dataset to build the decision tree allowing the prediction of an STU as a function of environmental covariates. C5.0 created explicit models, which were applied to the covariates rasters to generate a realisation of STU distribution over the study area at 50 m resolution.

These three steps were repeated 100 times to generate 100 realisations of the potential soil type distribution over the study area at 50 m of resolution. 4) Computing the probabilities of occurrence: the 100 realisations were stacked to calculate the probability of occurrence of each predicted STU by counting the frequency of each STU at each pixel. This procedure led to a set of 171 rasters depicting the probability of occurrence of 171 STU.

2.4.2) Original DSMART algorithm + soil observations (method 2)

This disaggregation approach is similar to the original DSMART algorithm. However, the main difference is that 755 additional soil profiles, spatially collocated, were added to the calibration dataset to build decision trees. These soil profiles make it possible to incorporate real field observations with established soil landscape relationships. For each realisation, a calibration dataset (15, 125 samples) including virtual samples randomly selected from polygon units, as well as soil observations were used to model soil type with environmental covariates. The model was then extrapolated over the study area.

2.4.3) Original DSMART algorithm + expert rules (Method 3)

Including soil landscape relationships in the disaggregation process was explored by [START_REF] Vincent | Spatial disaggregation of complex Soil Map Units at the regional scale based on soil-landscape relationships[END_REF] in a specific regional pedoclimatic context in Brittany (France). Expert soil landscape relationships were used to assign STU to sampling points. These relationships were based on expert pedological knowledge, which takes into account soil parental material as well as topography and waterlogging in the UTS allocation procedure. This approach combines two sources of the dataset to calibrate the model. The first one was derived from semantic information for each SMU/STU combination. It consists in attributing a barcode to each SMU/STU combination, derived from a concatenation of four features contained in the RRP database (parent material, SMU identifier, TPI and waterlogging index), and to compare these barcodes to a stack of regional covariates representing the same four features, to assign each pixel of the study area to a suitable STU. This procedure allowed matching soil exhibiting specific features with their potential spatial distribution. For instance, hydromophic soils occur with slope sequences and valley positions, while well drained soils occur in upslope or middle slope positions. Using a random sampling stratified by SMU's area, a set of sampling points was selected with a proportion of one sample for every 5 hectares and a minimum of five samples per polygon unit.

The second dataset was derived from a random sampling of a fixed number of sampling points in each polygon unit. This procedure ensured that all polygons had been sampled. STU allocation was based on the soil map unit proportions. The full set of each realisation (18, 320 samples) combining expert calibration dataset as well as dataset derived from random sampling procedure was spatially intersected with existing environmental covariates and used as a unique calibration dataset to build decision trees.

2.4.4) Prediction of the most probable STUs

From all soil type probability rasters obtained, only the three most probable STUs (with the highest probability of occurrence) were considered: for each pixel, the final prediction was the combination of the three most probable predicted STUs (1 st STU, 2 sd STU, and 3 rd STU) and their associated probability of occurrence.

The classification confusion index (CI) between the first most probable STU and the second most probable STU was calculated following Eq.1:

CI = 1-(P 1 st STU -P 2 sd STU) [1]
Where P1 st STU and P2 sd STU denote respectively the highest probability of occurrence for 1 st STU and the second highest probability of occurrence for 2 sd STU, calculated at each pixel [START_REF] Burrough | Continuous classification in soil survey: spatial correlation, confusion and boundaries[END_REF]Odgers et al., 2014).

This index was considered as an indicator of certainty assessment about the most probable predicted soil class and is ranging between 0 and 1. It tends to 1. When the 1 st STU and 2 sd STU are predicted with similar probability of occurrence and zero when the probability of occurrence of the 2 sd STU is close to zero.

2.5) Validation of disaggregated soil maps

The quality of soil maps resulting from the three DSMART algorithm based approaches was assessed by combining both spatial and semantical validation methods. Spatial validation is divided into 2 sub approaches ("pixel to pixel" and "window of 3x3 pixels"). For detailed soil maps and accurate soil profiles, "pixel to pixel" validation consists in checking, at each pixel, if the predicted STU respects the observed STU value [START_REF] Heung | Predictive soil parent material mapping at a regional-scale: A Random Forest approach[END_REF]Nauman et al., 2014;[START_REF] Chaney | POLARIS: A 30-meter probabilistic soil series map of the contiguous United States[END_REF][START_REF] Møller | Improved disaggregation of conventional soil maps[END_REF]. The "window of 3x3 pixels" validation assumes that, for each pixel, the predicted STU respects the observed STU value if it matches at least one of its 9 surrounding neighbours [START_REF] Heung | Predictive soil parent material mapping at a regional-scale: A Random Forest approach[END_REF][START_REF] Chaney | POLARIS: A 30-meter probabilistic soil series map of the contiguous United States[END_REF]. This method provides some flexibility by compensating spatial referencing error of soil maps and avoids the impact of fine scale spatial noise.

The semantical validation was also performed considering either each STU or a group of STUs sorted by expert on the basis of similar pedogenesis factors and similar diagnostic horizons [START_REF] Vincent | Spatial disaggregation of complex Soil Map Units at the regional scale based on soil-landscape relationships[END_REF][START_REF] Møller | Improved disaggregation of conventional soil maps[END_REF]. From the initial 171 STUs described in the soil database, the sorting procedure led to 78 groups and 11 STU remained single.

2.6) Pair wise comparisons of disaggregated soil maps

To compare the soil type rasters derived from the three DSMART based approaches, pairwise comparisons were performed using Vmeasure method implemented as open source software in an R package called Spatial Association Between REgionalisations (SABRE) [START_REF] Rosenberg | V-Measure: A Conditional Entropy-Based External Cluster Evaluation Measure[END_REF]. This is a spatial method developed to compare maps in the form of vector objects and it was commonly used in computer science to compare (non spatial) clustering. , the total area of a zone zj is

𝐴𝑗 = ∑ 𝑎 𝑖𝑗 𝑛 𝑖,1
and the total of the domain is

𝐴 = ∑ ∑ 𝑎 𝑖𝑗 𝑛 𝑖=1 𝑚 𝑗=1
.

The SABRE package calculates a degree of spatial agreement between two regionalizations using an information theoretical measure called the V measure. V measure provides two intermediate metrics: homogeneity and completeness. Homogeneity is a measure of how well regions from the first map fit inside zones from the second map (Eq 2). Completeness measures how well zones from the second map fit inside regions from the first map (Eq 5). The final value of V measure is calculated as the weighted harmonic mean of homogeneity and completeness (Eq 8). All metrics range between 0 and 1, where larger values indicate better spatial agreement. V measure, homogeneity, and completeness are global measures of association between the two regionalizations.

Additional indicators of disaggregation quality were calculated using Shannon entropy index of regions and zones [START_REF] Shannon | A mathematical theory of communication[END_REF]Nowosad and Stepinskie, 2018). These indicators qualify local associations by highlighting the region's inhomogeneities (Eq 3, Eq 4), or zone's inhomogeneities (Eq 6, Eq 7). Two normalized Shannon entropy was also computed using the ratios (𝑆 𝑗 𝑅 /𝑆 𝑅 ) and

(𝑆 𝑖 𝑍 /𝑆 𝑍 ) to derive maps of local spatial agreement between the two regionalizations R and Z. These measures have a range between 0 and 1.

When 𝑆 𝑗 𝑅 (Eq3) is close to zero, this denotes that the zone j is homogenous in terms of regions (each zone is within a single region). However, when 𝑆 𝑗 𝑅 value increases the zone is increasingly inhomogeneous in terms of regions (it overlays an increasing number of regions). Therefore, 𝑆 𝑗 𝑅 (Eq 3) assesses the degree of this inhomogeneity or a variance of region in zone j. A global indicator that measures a homogeneity of a given zone in terms of regions is given via Eq 2.

Analogous to homogeneity but with the roles of regions and zones reversed, the dispersion of zones over the entire area is also computed using Shannon entropy (Eq 4 and Eq 7), and a global indicator C (Eq 5) measures a homogeneity of a given region in terms of zones. ) log (

𝑎 𝑖 , 𝑗 𝐴 𝑗 ) [3] 𝑆 𝑅 = -∑ ( 𝐴 𝑖 𝐴 𝑛 𝑖=1
) log (

𝐴 𝑖 𝐴 ) [4] 𝑐 = 1 -∑ ( 𝐴 𝑖 𝐴 𝑛 𝑖=1 ) ( 𝑉𝑎𝑟𝑖𝑎𝑛𝑐𝑒 𝑜𝑓 𝑧𝑜𝑛𝑒𝑠 𝑖𝑛 𝑟𝑒𝑔𝑖𝑜𝑛 𝑖 =𝑆 𝑖 𝑍 𝑉𝑎𝑟𝑖𝑎𝑛𝑐𝑒 𝑜𝑓 𝑧𝑜𝑛𝑒𝑠 𝑖𝑛 𝑡ℎ𝑒 𝑑𝑜𝑚𝑎𝑖𝑛=𝑆 𝑍 ) [5] 𝑆 𝑖 𝑍 = -∑ ( 𝑎 𝑖,𝑗 𝐴 𝑖 𝑚 𝑗=1
) log (

𝑎 𝑖 , 𝑗 𝐴 𝑖 ) [6] 𝑆 𝑍 = -∑ ( 𝐴 𝑗 𝐴 𝑚 𝑗=1
) log (

𝐴 𝑗 𝐴 ) [7] 𝑉 β = (1+β)ℎ𝑐 (βh)+c [8]
β is a coefficient that allows promoting the first or the second regionalization, and by default, β equals 1. 𝑉 β has a range between 0 and 1. It equals 0 in case of no spatial association and 1 in case of perfect association.

The V measure method was applied in two main situations (DSMART+expert rules, Original DSMART) and (DSMART + expert rules, DSMART+extra soil observations). The reference map is always the map derived from DSMART algorithm with expert soil landscape relationships.

3) Results

3.1) Disaggregated soil maps

Applying DSMART based approaches yielded a set of soil maps and associated probability of occurrence rasters. The original DSMART approach allowed to disaggregate the 96 SMUs into 108 STUs while DSMART with expert rules approach yielded 158 STUs and DSMART with extra soil observations approach yielded 172 STUs with respect to the first most probable STU map. A total of 171 STUs were identified in the Ille et Vilaine department within the RRP database.

Unpredicted STUs correspond mainly to rare STUs with low proportions ranging between 2 and 10% within the SMUs containing them.

Figure 3 shows the three maps of the 1 st most probable STU derived from each approach as well as the original soil map. Overall, the three most probable STUs maps captured the main pattern of The quality of maps resulting from DSMART based approaches was quantified via the probabilities of occurrence of each STU predicted and the confusion index maps (Fig. 5). The latter measure indicated areas where the probability of occurrence of the two most probable soil types was close.

Over the study area, the average probability of occurrence of the most probable soil type achieved respectively 0.41 for DMSART map, 0.68 for DMSART with expert rules and 0.28 for DSMART with extra soil observations maps. Meanwhile, the average confusion index reached 0.8 for the original DSMART approach while DSMART with extra soil observations and DSMART with expert rules achieved 0.9 and 0.43, respectively. Although the most probable soil classes provide plausible maps of soil distribution, there is a significant prediction uncertainty as depicted by these measures.

In regions where disaggregated soil maps showed low confusion index, particularly in northwest and north coast areas of Ille et Vilaine department, high confidence in predictions was recorded.

These areas were predominantly deep loamy soils or developed in alluvial and colluvium deposits.

Figure 6 compares the cumulative area of the STUs estimated from the three disaggregated maps and that derived from the regional soil database. For each STU, its relative predicted area was estimated by counting the number of pixels where it was predicted. For the regional soil database, each STU area was computed from total SMU area multiplied by the proportion of the STU. This comparison shows that some STUs were overestimated by the disaggregation approaches when comparing to the soil database. DSMART with extra soil observations and original approaches showed similar cumulative STU areas under the curve whereas DSMART with expert rules had a shape similar to the regional soil database.

The most abundant STU in the database (431: Stagnic Fluvisol developed from alluvial and colluvium deposits) was predicted as the most frequent STU by DSMART with extra soil observations and DSMART with expert rules, and it was predicted as the second most abundant STU by the original DSMART algorithm. The 10 most abundant STUs in the soil database covers almost 43% of the study area. Of them, 7 belong to the 10 STU most predicted by the three disaggregation approaches (Table 2).

3.2) Covariates importance in the decision trees

Figure 7 gives the relative importance of the covariates used in DSMART based approaches. Soil parent material and SMU boundaries were used systematically in condition rules regardless of the disaggregation method. This was consistent with the contrasting pattern of geology and the dependence relationship between SMU and its soil components. Considering the original DSMART approach (Fig. 7.a), distribution functions of Aeolian silt deposits, airborne gamma ray spectrometry variable (K:Th ratio) and elevation contributions were more dispersed according to the STU considered than those of other covariates. For instance, Aeolian silt deposits contribution varied between 20 and 80% with a median value of 42%, whereas slope contribution ranged between 20 and 40 % with a median value of 28%. Aeolian silt deposits have an important weight in STU predictions, due to its ability to represent soils inherited from this superficial parent material, which is poorly represented in lithological maps.

DSMART with soil landscape relationships (Fig.

3.3) Validation of disaggregated soil maps

The validation procedure was performed for each DSMART based approach applied, considering the three most probable soil types and using both semantic objects (STU or soil group) and spatial neighbourhood (per pixel or 3x3 window of pixels).

Considering 755 legacy soil profiles prospected in the framework of "Sols de Bretagne" project, per pixel validation accuracy reached 27%, for original DSMART maps and 34 % for DSMART with expert rules (Table 3). A similar comparison using 135 validation sites derived from Soilserv project showed that 18.1 % of soil profiles match DSMART maps, 19.8 % match DSMART with expert rules maps and only 16.9 % match DSMART with extra soil observations maps (Table 3).

Using a 3 x 3 window of pixels markedly improves the global accuracies, which increased for the two validation datasets (Table 3). DSMART with soil landscape relationships remained the best performing method.

When compared to accurate soil maps (1:25,000), the validation procedure showed that DSMART with extra soil observations as well as DSMART with soil landscape expert rules had almost the same performance (37% and 38%) while best accuracy (44%) was observed for Original DSMART maps (44%) (Table 3). These scores were clearly improved by considering soil groups and 3x3 pixels neighbourhood. For instance, the accuracy of DSMART with expert rules maps using soil group reached 45.9% and increased to 62.1 % when considering 3x3 pixels windows (Table 3). and [START_REF] Møller | Improved disaggregation of conventional soil maps[END_REF] in Denmark (18%) using DSMART algorithm (Table 4). In contrast to the latter studies, a large number of STU (171 STU) compose our soil dataset. This could certainly decrease the chance of predicting the right STU, even through mobilizing relevant geographic dataset to implement soil landscape relationships.

3.4) Comparing disaggregated maps

When considering a window of 3x3 pixels, the overall accuracy increased considerably for the three DSMART based approaches maps, but DSMART with expert soil landscape relationships achieved the highest accuracy scores. [START_REF] Chaney | POLARIS: A 30-meter probabilistic soil series map of the contiguous United States[END_REF] highlighted a high degree of spatial noise in the predictions by including pixel validation neighbours. Overall, prediction accuracy increased twofold with a 3x3 pixel validation window and when grouping soils to a coarser level of soil classification (171 versus 89 soil group). This was recorded for all disaggregated maps regardless of the disaggregation procedure and suggests that fine soil taxonomic dissimilarities can not be accurately mapped by disaggregation processes.

4.2) Legacy soil data

Legacy soil data used in this study provide an overall representation of soil over large areas (1: 250,000 scale). This database was derived from several soil surveys and pedological expert knowledge. SMUs were spatially delineated, and their spatial organisation, as well as STUs features, were described according to available soil data and pedological expertise. STUs and their associated landscape characteristics were identified as accurately as possible using legacy soil profiles collected according to a not probabilistic sampling design between 1968 and 2012. Hence, differences in survey methods covering a large area over a long sampling period could lead to errors in the STU definition or uncertainties in the estimation of their area in a given SMU.

Moreover, soil survey intensity was not uniform within SMUs. Thus, SMU components may be derived from the unequal representation of soil samples across SMUs.

Harmonising soil data to reduce the number of STU is a great challenge by itself. Grouping some STUs regarding their pedological similarities such as sharing comparable morphological criteria, having similar pedogenic horizons and occurring in analogous environmental conditions is worthwhile to be investigated. More importantly, unifying soil data according to more functional aspects such as soil agricultural potential allows also to generate a relevant regional soil database easily handled by soil users to satisfy their needs. Many countries around the world have already harmonized their soil databases such as Denmark and Australia, where high pedological complexity was captured with a reasonable STU number, with not exceeding 23 soil groups in Denmark [START_REF] Møller | Improved disaggregation of conventional soil maps[END_REF] and 73 soil groups in Australia (Holmes et al., 2015).

4.3) Taxonomic similarities

In the recent DSM literature, DSMART approach is considered as an efficient tool to disaggregate existing coarse soil maps. In this study, we compared variants of the DSMART based approach, which differed by the training dataset used to calibrate the C5.0 model and the allocation procedure.

Modified DSMART algorithms used additional calibration datasets derived from supplement soil observations and expert sampling of polygons. Hence, taxonomic similarities were not taken into account neither in the calibration process nor in the current component assignment scheme. Even if there is a large number of STUs addressing inherent soil landscape heterogeneity, there is most likely a short taxonomic distance between many of them. As a result, these STUs may have similar forming conditions, making it a challenge to suitably constrain the prediction probabilities using DSMART algorithm. This likely explains the high confusion index scores recorded in the present study, particularly for original DSMART and DSMART with extra soil profiles approaches. As demonstrated by [START_REF] Minasny | Spatial prediction of soil properties using EBLUP with the Matérn covariance function[END_REF], including taxonomic distance in decision trees using pedological knowledge is a relevant way to decrease the misclassification error. Therefore, future effort and improvements of the DSMART algorithm should take into account the taxonomic distance between STU in the disaggregation procedure.

4.4) Mapping comparison

A quantitative comparison between disaggregated soil maps was performed using a novel approach called V measure method. This method was commonly used to assess the spatial agreement between land cover maps and thematic biotic and abiotic factors maps, as done by Nowosad and Stepinski (2018) in the United States, but never before for soil maps.

In the present study, V1 (0.53) was larger than V2 (0.47) suggesting that DSMART with expert soil landscape relationships map is much more similar to Original DSMART map than DSMART with extra soil observations map. This might be explained by the allocation procedure for training samples. The original DSMART algorithm tends to promote most abundant STUs with high proportions of occurrence within polygons and penalized STUs with low proportions (comprise between 2 and 10%). Therefore, frequent STUs are more likely to be predicted rather than rare STUs. Meanwhile, by adding supplement soil profiles, preliminarily assigned to a suitable STU to the training dataset, we constrain STUs with low proportions of occurrence predictions.

Major differences between DSMART with expert rules map and DMSART with soil observations were mainly observed in the southern part of the study area and valleys areas. In general, Fluvisol Stagnic soils were overestimated by DSMART with extra soil observations. This was likely due to the purposive sampling design followed to supplement soil observations. The 755 legacy soil profiles were selected to characterize hydromorphic soil conditions and to characterize inherent soil landscape variability supposed to be organized along the hillslope.

4.5) Improvements and future work

Even though this work emphasizes the contribution of pedological knowledge in the disaggregation process, other pathways can also be explored to improve map's accuracy. As recommended by [START_REF] Mulder | National versus global modelling the 3D distribution of soil organic carbon in mainland France[END_REF], compensating the temporal changes and differences in laboratory analytics is a good option to improve the quality of legacy soil data. This suggests harmonising local soil database and regrouping some STUs with similar soil forming factors through statistical modelling.

Moreover, additional environmental covariates with high spatial resolution should be used to capture micro landscape variability [START_REF] Lacoste | High resolution 3D mapping of soil organic carbon in a heterogeneous agricultural landscape[END_REF]Odgers et al., 2014;[START_REF] Chaney | POLARIS: A 30-meter probabilistic soil series map of the contiguous United States[END_REF][START_REF] Møller | Improved disaggregation of conventional soil maps[END_REF]. For example, adding a more detailed Digital Elevation Model allowed to capture small terrain features, where may be particular, STUs occurs. Improving both polygon sampling procedure and current components assignment scheme turned out to be important to reduce uncertainty prediction. This suggests drawing virtual soil samples proportionally to polygons areas and using supplement STU characteristics based on surveyor observations (slope shape, hillslope position, soil texture …) to guide STU allocation procedure [START_REF] Møller | Improved disaggregation of conventional soil maps[END_REF].

Assuming that the decision tree can be built to relate STU descriptors to legacy soil data, this method can replace weighted random allocation procedure and should help minor STU prediction by constraining raster probabilities.

5) Conclusion

We quantitative.

Table 2. Ten most extended STUs according to the regional soil database and their respective rank by area using three DSMART based disaggregation procedures Table 3. Overall accuracies (%) obtained using various external validation approaches for the three most probable STU

  Brittany, France (48°N, 2° W) (Fig 1). It is drained by the rivers Ille and Vilaine and their tributaries. Its climate is oceanic, with a mean annual rainfall of 669 mm and mean annual temperature of 11.3° (Source: Climate Data EU). Main land uses comprise arable land, temporary and permanent grasslands, woodland, and urban areas. In the present study, anthropogenic areas were not considered. Elevation ranges between 0_20 m in the coastal zone and 20_150 m almost everywhere expect in the western part of the department where it tills 256 m. The topography is generally gentle with maximum slopes not exceeding 16%. The Ille et Vilaine department is part of the Armorican Massif with complex geology (BRGM, 2009): intrusive rocks (granite, gneiss and micaschist) in northern and north western zones, sedimentary rocks (sandstone) and metamorphic rocks (Brioverian schist) in the central and southern zones, and superficial deposits (Aeolian loam with decreasing thickness from north to south overlaying bedrock, alluvial and colluvium deposits). According to the World Reference Base of Soil Resources, soils occurring in Ille et Vilaine include Cambisols, Luvisols Stagnic Fluvisols, Histosols, Podzols, and Leptosols (IUSS Working Group WRB, 2014). 2.2) Soil data 5 https://doi.org/10.5194/soil-2019-36 Preprint. Discussion started: 7 June 2019 c Author(s) 2019. CC BY 4.0 License.

2. 4 )

 4 Disaggregation procedure: DSMART algorithm 2.4.1) Original DSMART algorithm (Method 1)

  We divide the entire study area into 2 different sets of regions, referred to as regionalizations R and Z. The first regionalization R divides the domain into n regions ri (i=1 to n) and the second regionalization Z divides the domain into m zones zj (j=1 to m). Superposition of the 2 regionalization R and Z divides the domain into n x m segments having aij area. The total area of a region ri is 𝐴 𝑖 = ∑ 𝑎 𝑖𝑗 𝑚 𝑗,1

  11 https://doi.org/10.5194/soil-2019-36 Preprint. Discussion started: 7 June 2019 c Author(s) 2019. CC BY 4.0 License.

  soil distribution of the coarse soil map. As one could expect according to the geological parent material map(Lacoste et al., 2011), extensive areas of deep silty soils are developed in Aeolian loam deposits encountered in the north east as well as in the north central parts of the study area.Colluvial and alluvial soils were mainly predicted in the north coast part and large valleys zones. The visual comparison of disaggregated soil maps highlighted global similarities in the soil spatial distribution markedly affected by SMU boundaries. The three approaches distinguished very well soils developed in marsh parent material in the coastal part (north) of the study area. However, DSMART with soil landscape expert rules map as well as DSMART with extra soil observations map remained more detailed and underlined a clear internal disaggregation of SMUs especially in the north and the central parts of the Ille et Vilaine department. Visual inspection of the obtained DSMART with extra soil observations map as well as DSMART with expert rules map showed an increase in soil heterogeneity when compared to Original DSMART map. More importantly, legacy soil profiles made it possible to take into account some rare soil types with low probability to be predicted. Therefore, adding supplement sampling points via the expert calibration dataset and the 755 extra soil profiles allowed to predict STUs characterized in the soil database with a low spatial extent. Nevertheless, the three DSMART based approaches spatially disaggregated the most frequent components disregarding the less frequent ones.

Figure 4

 4 Figure4shows maps of the global probability of redoximorphic soils across the study area. STU probability rasters, depicting hydromorphic soils, were added together to produce continuous maps of hydromorphic soil probability. Visual inspection of three maps highlighted global similarities, but local differences were recorded along the hydrographic network and in the southern part of the study area. As could be expected, DSMART with expert rules well predicted hydromorphic soils in valleys and coastal areas, with a probability of occurrence exceeding 80%. Adding soil landscape relationships in the allocation process constrained hydromorphic soil predictions in specific landscape positions. The same trend characterized DSMART with extra soil observations map, particularly in the central part of the study area. Therefore, including 755 soil profiles had an important role in the disaggregation process in the northern and the central parts where these profiles were located.

  7.b) showed almost the same distribution function of all covariates except for elevation where its distribution function was more dispersed. Since a part of training samples was chosen with expert knowledge based on three environmental covariates: TPI, a waterlogging index and soil parent material, we would expect the prominent role of waterlogging index and TPI to constrain hydromorphic soils predictions and to achieve STU distribution in the appropriate order along the toposequence. This most likely explains the dominance of Fluvisol Stagnic in valleys areas followed by a transition to Cambisols commonly found at upslope and midslope positions along the toposequences. Analogous to the original DSMART algorithm, DSMART with extra soil observations (Fig. 7.c) highlighted almost the same distribution of use of soil environmental covariates in the decision trees, except for aeolian silt deposits, K:Th ratio and elevation. The latter covariates contributions remained less dispersed compared to the original DSMART approach.

Figure 8

 8 Figure8shows inhomogeneity maps measured by Shannon entropy. The map derived from DSMART with soil landscape relationships was chosen as a reference map. This map deeply disaggregates the initial SMUs into 120,653 regions with irregular shapes. By contrast, Original DSMART map remained very similar to the original map and delineated the study into 40,459 regions. Both disaggregated maps reflect the main pattern of soil distribution over the study area despite the difference in the disaggregation process. Visual inspection of maps DSMART with soil landscape rules map and Original DSMART map revealed an overall similarity between disaggregated maps, but local differences between them were depicted.We calculated h1= 0.49, c1 = 0.58 and V1 = 0.53 as global measures of spatial agreement between the two maps (DSMART+expert rules and Original DSMART). The average homogeneity of the DSMART with soil landscape rules map with respect to the Original DSMART map was qualified via h homogeneity index. Similarly, the average homogeneity of the Original DSMART map with respect to the DSMART with soil landscape rules map was qualified via c completeness index.Visually, the Fig.8.b map seemed to be more homogeneous than the map Fig.8.a in agreement with the statistical assessment c > h. The large number of DSMART with soil landscape rules map regions, which was three times higher than Original DSMART map zones, might explain this difference. It is more likely that DSMART with soil landscape rules map regions cross through multiple Original DSMART map zones than vice versa. However, two disaggregated maps remained spatially associated according to the high V1 score. The two inhomogeneity maps (Figs.

  8a and 8b) highlighted the locations of greatest differences between two maps, mainly along the hydrographic network.When comparing disaggregated soil maps derived from modified DSMART algorithm (DSMART with soil landscape rules and DSMART with supplement soil observations), we note that the DSMART with extra soil observations map delineated the study area into 132,942 regions. For both maps, internal disaggregation was well pronounced expect for DSMART with extra soil observations map in the southern part of the study area. Visual inspection of selected maps showed high spatial agreement and highlighted some locations of greatest differences, particularly in the southern part of the Ille et Vilaine department. Even if the hydrographic network was well detailed in both maps, it appeared more developed in DSMART with extra soil observations soil map.Applying V measure method for assessing the spatial similarity between DSMART with soil landscape rules map and DSMART with supplement soil observations map provided similar information theoretical measures h2 = 0.47, c2 =0.48, and V2 = 0.47. Visual comparison of soil inhomogeneity maps revealed constant variance measured by normalized Shannon entropy. This was in agreement with the quantitative assessment c = h. Overall, the two disaggregated maps were spatially correlated, as indicated by the global spatial agreement measure V2.4) Discussion4.1) Performance of the disaggregation proceduresProduced disaggregated soil maps closely resemble the abundant soils in the original soil map(Holmes et al., 2015; Fig.3). The 1 st most probable STU map derived from DSMART based approaches captured the main spatial pattern of soil distribution across the study area. More internal variation within SMUs was found when using DSMART with added point observations and DSMART with soil landscape relationships. Local soil heterogeneity reflecting inherent pedological complexity was depicted by the 1 st STU maps which deliver a deterministic soil landscape distribution, continuously varying with landscape features.External validation was performed to assess the quality of disaggregated soil maps. Using 135 independent soil profiles and a per pixel validation approach, the overall accuracy reached 18.1% for DSMART algorithm 1 st STU map, 19.8% for DSMART with expert rules 1 st STU map and 16.9% for DSMART with extra soil profiles 1 st STU map. In the DSM literature, researchers who applied classification tree decision methods founded similar validation results. For instance, by applying DSMART algorithm in eastern Australia and using 285 legacy soil profiles, Odgers et al.(2014) achieved an overall accuracy of 23%. Similarly,Nauman and Thompson (2014) explored the use of expert rules for soil landscape relationships in the United States and achieved global accuracy ranged between 22% and 24%. Similar disaggregation performance was recorded byHolmes et al. (2015) in Western Australia (20%),[START_REF] Chaney | POLARIS: A 30-meter probabilistic soil series map of the contiguous United States[END_REF] in the United States (17%)

  applied three DSMART based approaches, including original DSMART algorithm, DSMART with extra soil observations and DSMART with soil landscape relationships, to disaggregate legacy soil polygons over a large area in Brittany (France). Regardless of the disaggregation approach, the produced soil maps at 50 m spatial resolution successfully address the main soil spatial pattern regarding prior pedological knowledge of our study area. Performance assessed against 135 independent soil profiles, 755 legacy soil profiles, and accurate 1:25,000 soil maps highlighted that DSMART with expert rules maps achieved highest validation measures. Overall, modified DSMART algorithms allowed minor STUs prediction, whereas original DSMART algorithm promoted abundant STUs prediction with poor spatial structure improvement. Adding pedological knowledge as well as extra soil observations in the prediction process constrained STU probabilities, even STUs with low proportions. However, some particular STUs reflecting hydromorphic soils or loamy soils were greatly overestimated for all the three DSMART based approaches.Soil maps produced using the original DSMART and DSMART with expert rules have a high spatial agreement, but the latter map appeared more detailed and provided a spatially continuous and consistent STU's prediction. Therefore, generalizing soil landscape relationships taken to account several STU descriptors and landscape features should be implemented in the future version of DSMART algorithm to capture soil landscape heterogeneity and consequently guarantee coherent variability of soil properties.

Figure captions Figure 1 :Figure 2 :Figure 4 :Figure 5 :Figure 6 :Figure 7 :

 captions124567 Figure captions

Figure 1 :

 1 Figure 1: Location of the study area and the validation datasets

Figure 2 :

 2 Figure 2: Schematic of the DSMART based approaches algorithm. The steps in DSMART are: 1) construct the calibration dataset; 2) train C5.0 model; 3) estimate STU maps and their associated probabilities of occurrence

Figure 3 :

 3 Figure 3: Digital soil map of the most probable STU and their associated probability of occurrence for the whole study area and for a focus zone, a) Legacy soil map: most probable STU for each SMU, b) original DSMART approach; c) DSMART with expert rules; d) DSMART with extra soil observations

Figure 4 :

 4 Figure 4: Global probability of hydromorphic soils over the study area derived from a) original DSMART, b) DSMART with soil landscape relationships and c) DSMART with extra soil observations. The probabilities of the three STU with highest prediction occurrence are summed if they are hydromorphic.

Figure 5

 5 Figure 5 Confusion index maps for a) Classic DSMART approach; b) DSMART with expert rules; c) DSMART with extra soil observations

Figure 7 :

 7 Figure 7: Violin plots of the relative importance of each environmental covariate used in a) Original DSMART approach; b) DSMART with expert rules; c) DSMART with extra soil observations

Figure 8 :

 8 Figure 8: Spatial association between disaggregated maps of Ille et Vilaine department. a) map of inhomogeneity of DSMART with soil landscape relationships map in terms of original DSMART map b) map of inhomogeneity of original DSMART map in terms of DSMART with soil landscape relationships map c) map of inhomogeneity of DSMART with soil landscape relationships map in terms of DSMART with extra soil observations map d) map of inhomogeneity of DSMART with extra soil observations map in terms of DSMART with soil landscape relationships map. Inhomogeneity (variance) is measured by normalised Shannon entropy

Table headings Table 1 .

 headings1 Description of the environmental covariates selected. Summary of environmental

	covariates. P: parent material; S: soil properties; R: relief; O: Organisms; C: categorical; Q:

Table 4 :

 4 Comparison between the size areas covered, number of soil map units, soil type units of the original legacy soil maps and the accuracy achieved in other studies using DSMART algorithm

24 https://doi.org/10.5194/soil-2019-36 Preprint. Discussion started: 7 June 2019 c Author(s) 2019. CC BY 4.0 License.

Table 1 .

 1 Description of the environmental covariates selected Summary of environmental covariates. P: parent material; S: soil properties; R: relief; O: Organisms; C: categorical; Q: quantitative.

Table 3 .

 3 Overall accuracies (%) obtained using various external validation approaches for the three most probable STU

		Pixel to pixel validation of STU
		DSMART approach	Most probable STU Second most	Third most	Total
			probable STU	probable STU
		Original DSMART	23	13	8	44
	Soil maps (87 150 ha)	DSMART with expert rules DSMART with extra soil	19 22	11 9	7 7	37 38
		observations
		Original DSMART	11	5	3.8	18.1
	Independent soil profiles (n=135)	DSMART with expert rules observations DSMART with extra soil	10 8.2	4.4 6	3.7 2.7	19.8 16.9
		Original DSMART	14	7	6	27
	Legacy soil profiles (n=755)	DSMART with expert rules DSMART with extra soil observations	18	9	7	34

factor Type Unit or number of classes Terrain attributes derived from the digital elevation model

  doi.org/10.5194/soil-2019-36 Preprint. Discussion started: 7 June 2019 c Author(s) 2019. CC BY 4.0 License. doi.org/10.5194/soil-2019-36 Preprint. Discussion started: 7 June 2019 c Author(s) 2019. CC BY 4.0 License.

		Pixel to pixel validation of STU group
		DSMART	Most probable STU	Second most	Third most	Total
		approach				probable STU	probable STU
		Original DSMART			26		13	9	48
	Soil maps (87 150 ha)	DSMART with expert rules DSMART with extra soil		22.5 25		13.7 10	9.7 7	45.9 42
	Elevation Slope Compound Topographic Index (TPI) R R observations R Original DSMART Topographic Position Index R Pedology and geology DSMART with expert rules Independent soil profiles observations (n=135) DSMART with extra soil	16 18 15	Q Q Q C	m 7 % Log (m 3 ) 5 classes 8.4 8	4.6 5.2 3.8	27.6 31.6 26.8
	Soil parent material Original DSMART		P	19	C	22 classes 12	9	40
	Soil Map Units Aeolian silt deposits Waterlogging index DSMART with expert rules Legacy soil profiles DSMART with extra soil (n=755) observations	R P S	23.4	C C C	96 classes 2 classes 4 classes 15	11.8	50.2
	Organism
	Landscape units			O		C	19 classes
	Gamma ray spectrometry from 250 m airborne geophysical survey interpolations
	K:Th ratio	P Neighbourhood of 3 x 3 validation of STU Q
		DSMART approach		Most probable STU Second most	Third most	Total
						probable STU	probable STU
		Original DSMART			31		16	14	61
	Soil maps (87 150 ha)	DSMART with expert rules DSMART with extra soil		29.6 28		19.4 11	13.1 9	62.1 48
		observations
		Original DSMART			15		6	4.3	25.3
	Independent soil profiles (n=135)	DSMART with expert rules observations DSMART with extra soil		17 11		6.7 7	4.8 3	28.5 21
		Original DSMART			19		10	7	36
	Legacy soil profiles (n=755)	DSMART with expert rules DSMART with extra soil observations		27.9		15	11.9	54.8

37 https://38 https://

Table 4 :

 4 Comparison between the size areas covered, number of soil map units, soil type units of the original legacy soil maps and the accuracy achieved in other studies using DSMART algorithm

	Study	Area (km²)	Map units	Soil type unit	Accuracy
	Odgers et al (2014)	68,000	1,110	72	23

https://doi.org/10.5194/soil-2019-36 Preprint. Discussion started: 7 June 2019 c Author(s) 2019. CC BY 4.0 License.

https://doi.org/10.5194/soil-2019-36 Preprint. Discussion started:

June 2019 c Author(s) 2019. CC BY 4.0 License.

https://doi.org/10.5194/soil-2019-36 Preprint. Discussion started: 7 June 2019 c Author(s) 2019. CC BY 4.0 License.

https://doi.org/10.5194/soil-2019-36 Preprint. Discussion started: 7 June 2019 c Author(s) 2019. CC BY 4.0 License.

https://doi.org/10.5194/soil-2019-36 Preprint. Discussion started: 7 June 2019 c Author(s) 2019. CC BY 4.0 License.

https://doi.org/10.5194/soil-2019-36 Preprint. Discussion started: 7 June 2019 c Author(s) 2019. CC BY 4.0 License.

https://doi.org/10.5194/soil-2019-36 Preprint. Discussion started: 7 June 2019 c Author(s) 2019. CC BY 4.0 License.

https://doi.org/10.5194/soil-2019-36 Preprint. Discussion started: 7 June 2019 c Author(s) 2019. CC BY 4.0 License.

Acknowledgments

The authors gratefully acknowledge all farmers at the Ille et Vilaine site involved in our research.

We thank the technical staff who actively participated in field sampling and laboratory analysis.

This research was performed in the framework of the INRA "Ecoserv" metaprogram. This work was also supported by Sols de Bretagne project and Soilserv program funded by ANR (Agence Nationale de la Recherche) (ANR-16-CE32-0005-01).


B Cotté
email: benjamin.cotte@ensta-paristech.fr

B. Cotté

Wind turbine noise extended

I. INTRODUCTION

Wind turbine noise can be perceived at distances greater than one kilometer and is characterized by amplitude modulations at the receiver [START_REF] Larsson | Amplitude modulation of sound from wind turbines under various meteorological conditions[END_REF]Zajamsek et al., 2016). As noise restrictions limit the areas where onshore wind farms can be built, an accurate prediction of the far-field noise is needed in order to improve the placement of the turbines at a given site, as well as to develop noise mitigation methods. This requires to model the dominant broadband aerodynamic noise sources as well as the main outdoor sound propagation effects that occur between the wind turbines and the receivers. The main aerodynamic noise sources are generally considered to be turbulent inflow noise, corresponding to the interaction of atmospheric turbulence with the blade leading edge, and trailing edge noise, corresponding to the scattering of the turbulent boundary layer at the blade trailing edge. As shown in the experimental campaign of Buck et al. (2016), turbulent inflow noise is generally dominant at low frequencies, typically below 300-400 Hz for a modern upwind turbine, while trailing edge noise dominates at higher frequencies, as already shown by [START_REF] Oerlemans | Prediction of wind turbine noise and validation against experiment[END_REF].

To model aerodynamic noise sources, the state-of-the-art approach is to divide the wind turbine blades into radial segments, and to sum incoherently the noise contributions from each segment at the receiver locations [START_REF] Oerlemans | Prediction of wind turbine noise and validation against experiment[END_REF]Zhu et al., 2005).

To model atmospheric propagation, however, this approach is rarely used, and it is more common to model the wind turbine as a point source of specified power located at the rotor center [START_REF] Lee | Prediction of far-field wind turbine noise propagation with parabolic equation[END_REF][START_REF] Prospathopoulos | Application of a ray theory model to the prediction of noise emissions from isolated wind turbines and wind parks[END_REF]. Recently, several methods have been proposed to include an extended source model in wind turbine noise propagation calculations. [START_REF] Mcbride | A comprehensive hamiltonian ray tracing technique for wind turbine noise propagation under arbitrary weather conditions[END_REF] and [START_REF] Heimann | Sound propagation through the wake flow of a hilltop wind turbine -a numerical study[END_REF] have considered extended models in ray-based models. [START_REF] Mcbride | A comprehensive hamiltonian ray tracing technique for wind turbine noise propagation under arbitrary weather conditions[END_REF] have kept all the blade segments used in the aeroacoustic source model in their 3D ray-tracing approach, while [START_REF] Heimann | Sound propagation through the wake flow of a hilltop wind turbine -a numerical study[END_REF] consider 24 fixed point sources distributed over the rotor disk with an identical sound power level (no source model used). One of the known weaknesses of these ray-based models is the treatment of diffraction, for instance in the presence of an acoustic shadow zone. Other authors have proposed methods based on the parabolic equation (PE), that is able to treat diffraction effects accurately. Barlas et al. (2017) have considered a PE model considering only one point source per blade. This point source is located at the segment location where the maximum noise level is calculated by their aerodynamic noise source model for each frequency. Cotté (2018) has kept several segments along the blade, and has used the backpropagation method to preserve the directivity of the noise sources, which makes the approach very computationally demanding.

In this study, two methods are compared to include extended aeroacoustic source models in a parabolic equation code for acoustic propagation in an inhomogeneous atmosphere. The source model is based on Amiet's theory [START_REF] Roger | Extensions and limitations of analytical airfoil broadband noise models[END_REF][START_REF] Tian | Wind turbine noise modeling based on Amiet's theory: Effects of wind shear and atmospheric turbulence[END_REF], but the methods could be applied to other source models, such as the so-called BPM semiempirical model that is widely used in wind turbine noise prediction studies [START_REF] Oerlemans | Prediction of wind turbine noise and validation against experiment[END_REF]Zhu et al., 2005). In the first method, that was recently proposed (Cotté, 2018), an initial starter for the PE model is obtained for each segment of the blade using the backpropagation approach. In the second method, that is introduced in the present study, the blade segments are viewed as moving monopole sources, and only a limited number of parabolic equation simulations are needed which strongly reduces the computation time.

The two methods are compared to the point source approximation, first in a homogeneous medium for validation purposes and second in a stratified atmosphere.

The paper is organized as follows. In Section II, the source and propagation models are described, as well as the proposed extended source methods. Then, the models are compared in Section III, first in a homogeneous atmosphere to be validated against an analytical solution, and then in an inhomogeneous atmosphere to take into account refraction effects in different directions of propagation.

II. DESCRIPTION OF THE EXTENDED SOURCE MODELS

A. Description of the broadband noise sources using strip theory It is common in the literature to calculate the noise spectrum of a full blade using strip theory (Christophe et al., 2009;[START_REF] Rozenberg | Rotating Blade Trailing-Edge Noise: Experimental Validation of Analytical Model[END_REF][START_REF] Sinayoko | Trailing edge noise theory for rotating blades in uniform flow[END_REF]. This theory consists in dividing the blade into M small segments or strips along the radial direction in order to take into account the variation of the blade geometry and the incident flow, as schematically shown in Fig. 1 For each segment at each angular position β, the power spectral density (PSD) of the acoustic pressure p for the rotating airfoil at angular frequency ω is written [START_REF] Sinayoko | Trailing edge noise theory for rotating blades in uniform flow[END_REF][START_REF] Tian | Wind turbine noise modeling based on Amiet's theory: Effects of wind shear and atmospheric turbulence[END_REF]:

S R pp (x T R , ω, β) = ω e ω S F pp (x B R , ω e , β), (1)
with ω e the emission angular frequency, x T R the receiver coordinates in the wind turbine reference system, x B R the receiver coordinates in the blade reference system, and S F pp the PSD for an airfoil that is fixed relative to the receiver. The expression for the Doppler factor ω/ω e is given in [START_REF] Sinayoko | Trailing edge noise theory for rotating blades in uniform flow[END_REF].

In order to calculate the PSD of acoustic pressure S F pp for an airfoil that is fixed relative to the receiver, various methods have been proposed in the literature. In this study, a model of trailing edge noise and turbulent inflow noise for wind turbines based on Amiet's theory is used, that is detailed in [START_REF] Tian | Wind turbine noise modeling based on Amiet's theory: Effects of wind shear and atmospheric turbulence[END_REF]. The model is valid for an aspect ratio L m /c m ≥ 3, m = 1..M , and the receiver is supposed to be in the far-field.

B. Acoustic propagation model based on the parabolic approximation

The acoustic propagation model considered here is a parabolic equation with fractional steps, called split-step Padé, based on higher order Padé approximants and solved with the method of Collins (1993). It has been shown in Cotté (2018) for a typical wind turbine configuration that it is more computationally effective than a classical wide-angle parabolic equation based on a Padé (1,1) approximation of the propagation operator. In this study, the effective sound speed approximation is used, which allows one to take into account the refraction effects due to the vertical wind gradients in the equation for a medium at rest:

c eff (z) = c(z) + U (z) cos τ = γ 0 rT (z) + U (z) cos τ, (2)
with z the height above ground, γ 0 the specific-heat ratio, r the specific gas constant, U (z)

and T (z) the mean vertical profiles of wind speed and temperature, and τ the angle between the wind direction and the propagation direction from the source to the receiver noted as x, as shown in Fig. 1(b). In order to introduce the notations needed for the extended source models described in Secs. II C and II D, the main equations of the model are briefly reminded below.

Using the axisymmetric approximation, the three-dimensional Helmholtz equation can be reduced to the following two-dimensional equation in the far-field:

∂ 2 ∂x 2 + ∂ ∂z 2 + k 2 q c = 0, (3)
where q c = p c √ x connects the q c variable to the complex pressure p c , and k is the acoustic wavenumber. This wavenumber can be written as k 2 = k 2 0 n 2 = k 2 0 (1 + ), where n(z) = c 0 /c eff (z) is the index of refraction and k 0 = ω/c 0 is the value of the acoustic wavenumber at the reference sound speed c 0 . Introducing the propagation operator

Q = 1 + + 1 k 2 0 ∂ ∂z 2 1/2 = (1 + L) 1/2 , (4)
that is independent of x in range-independent media, Eq. ( 5) becomes [START_REF] Gilbert | Application of the parabolic equation to sound propagation in a refracting atmosphere[END_REF]:

∂ ∂x + ik 0 Q ∂ ∂x -ik 0 Q q c = 0. ( 5
)
We can decouple Eq. ( 5) into two equations characterizing a wave propagating in the positive

x direction, denoted as q + (propagating wave), and a wave propagating in the negative x direction denoted as q -(backpropagating wave). Using the notation γ = ±1, one obtains from Eq. ( 5):

∂ ∂x -iγk 0 Q q γ = 0. ( 6
)
Introducing the variable φ γ corresponding to the envelope of the pressure:

q γ (x, z) = φ γ (x, z) exp(iγk 0 x), (7)
and substituting Eq. ( 7) into Eq. ( 6), the following equation is obtained:

∂φ γ ∂x = iγk 0 (Q -1) φ γ . (8)
To solve Eq. ( 8), the domain is discretized using a rectangular mesh of size ∆x and ∆z along the x and z-axis respectively. The split-step Padé (N,N) method is used to advance the field from x to x + ∆x for γ = 1, or from x to x -∆x for γ = -1 (Collins, 1993;Dallois et al., 2002).The angular validity increases with the order N of the development and depends on the mesh size ∆x chosen (Dallois et al., 2002). As shown in Cotté ( 2018), accurate results are obtained in a typical wind turbine configuration with N = 2 and mesh sizes ∆x = 2λ and ∆z = λ/10, where λ is the acoustic wavelength. Along the vertical direction, the domain is bounded by a ground impedance condition at z = 0, and by an absorbing layer at the top of the domain to obtain non-reflecting boundary conditions [START_REF] Salomons | Computational Atmospheric Acoustics[END_REF].

C. Extended source model based on the backpropagation method (Amiet-PE model)

The first extended source model considered in this study is based on the parabolic equation property to decouple forward and backward-propagating waves, as shown in Eqs. ( 6) and ( 8). It was proposed by Cotté (2018) and is called Amiet-PE model. The basics of the method are summarized in this section.

For each segment m, each angular position β of the blade and each angular frequency ω, a parabolic equation calculation is performed for which an initial condition at x = 0 is needed. This initial condition is obtained numerically using the backpropagation method, whose principle is illustrated in Fig. 2(a). It consists first in back-propagating a known pressure field, that is noted "initial solution" in Fig. 2(a), at x = x is to x = 0, taking γ = -1 in the equations. Then, in a second step, the starter at x = 0 is forward-propagated to the desired distance using the "classical" parabolic equation with γ = +1.

In the backpropagation method, the initial solution at x = x is is obtained for heights z is,p = p∆z, p = 0..P , from the expression (1) for the PSD of acoustic pressure of a rotating Color online.

blade. The initial solution thus includes the source directivity as viewed by this vertical line of receivers in this specific direction. In the presence of ground, the initial solution is written:

q c (z is,p ) = S R pp (x T R , ω, β) √ x S e ik 0 R 1,p × 1 + Q R 1,p R 2,p e ik 0 (R 2,p -R 1,p ) , (9)
where

R 1,p = r 2 is + (z S -z is,p ) 2 and R 2,p = r 2 is + (z S + z is,p ) 2
are respectively the distance between the segment at (0, y S , z S ) or the image segment at (0, y S , -z S ) and the p th initial starter point, with r is = x 2 is + (y S -y is ) 2 , and Q is the spherical wave reflection coefficient.

During the backpropagation calculation, the atmosphere is supposed homogeneous and the ground is taken as rigid (Q = 1). Then, the starter at x = 0 can be propagated using any ground impedance and any sound speed profile. Note also that each calculation is performed in a slightly different plane that crosses the far-field receiver at x = x R , as shown in Fig. 2(b). This method is therefore strictly exact only at this distance. For x = x R , the total acoustic pressure is obtained by summing contributions with different y values. Since the radius of the rotor (typically 50 m) is generally small compared to the propagation distances considered, the method remains valid over a wide range of distances. The computational cost of this method is quite high, since M N β PE calculations per frequency and per propagation direction need to be performed, where N β is the number of angular positions used to discretize the rotor plane.

D. Extended source model based on moving monopoles (MM model)

The second extended source model considered in this study represents each segment of the blade as a monopole rotating at angular velocity β. It is called the moving monopoles (MM) model. Compared to the first method, it does not rely on the parabolic approximation and can be applied to any propagation model.

In the MM model, the sound pressure level (SPL) at the receiver is calculated for a segment m at angular position β using the point source approximation [START_REF] Salomons | Computational Atmospheric Acoustics[END_REF]:

SPL(ω, β) = SWL(ω, β) -10 log 10 (4πR 2 1 ) + ∆L(ω, β) -α(ω)R 1 , (10)
where SWL(ω, β) is the angle-dependent sound power level (SWL),

R 1 = x 2 + y 2 S + (z S -z) 2
is the distance between the segment at (0, y S , z S ) and the receiver at (x, 0, z), ∆L is the sound pressure relative to the free field, and α is the absorption coefficient in dB/m.

The angle-dependent SWL can be obtained from the free-field SPL calculated using Amiet's model. Assuming free-field conditions (∆L = 0) and no absorption in the medium, Eq. ( 10) becomes:

SWL(ω, β) = SPL F F (ω, β) + 10 log 10 (4πR 2 1 ) = 10 log 10 S R pp (x T R , ω, β) p 2 ref + 10 log 10 (4πR 2 1 ), (11)
with SPL F F the free-field SPL and p ref = 20 µPa the reference pressure. From Eqs. ( 10) and ( 11), the following equation for SPL(ω, β) is obtained:

SPL(ω, β) = 10 log 10 S R pp (x T R , ω, β) p 2 ref + ∆L(ω, β) -α(ω)R 1 . (12)
In Eq. ( 12), the main unknown is the relative sound pressure level ∆L(ω, β). For the propagation over a finite impedance ground in a homogeneous atmosphere at rest, it can be calculated analytically [START_REF] Salomons | Computational Atmospheric Acoustics[END_REF]:

∆L = 10 log 10 1 + Q R 1 R 2 e ik 0 (R 2 -R 1 ) 2 , ( 13
)
with R 2 = x 2 + y 2 S + (z S + z) 2 the distance between the image-source and the receiver.

In order to include refraction effects, ∆L(ω, β) can be calculated using the parabolic approximation method described in Sec. II B. The initial starter corresponding to the monopole source is calculated numerically using the backpropagation method in order to preserve the angular validity of the split-step Padé (2,2) method [START_REF] Galindo | Approximations in the PE method. Phase and level errors in a downward refracting atmosphere[END_REF]. In order to limit the number of PE calculations to perform, a set of N h source heights distributed along the rotor plane are considered:

H n = H min + n∆H, n = 0, .., N h -1, (14)
with ∆H the height step given by:

∆H = H max -H min N h -1 , ( 15
)
where H min and H max are respectively the minimum and maximum heights to consider. The relative sound pressure level ∆L(ω, β) in Eq. ( 12) is then obtained using a nearest-neighbor interpolation. As an example, the monopole sources are represented for the three blades in β using Eq. ( 12), which will be necessary to calculate the amplitude modulation in Sec. III.

This means that the point source approximation is only used to account for propagation effects in the present study.

III. RESULTS USING BOTH EXTENDED SOURCE MODELS A. Configurations studied

In this study, the same 2.3 MW wind turbine as in [START_REF] Tian | Wind turbine noise modeling based on Amiet's theory: Effects of wind shear and atmospheric turbulence[END_REF] and Cotté ( 2018) is considered, with a diameter of 93 m, a hub height of 80 m and three blades of length 45 m.

As justified in [START_REF] Tian | Wind turbine noise modeling based on Amiet's theory: Effects of wind shear and atmospheric turbulence[END_REF], each blade is decomposed into M = 8 segments to respect the constraint on the aspect ratio L m /c m ≥ 3, m = 1, .., M mentioned in Sec. II A.

The rotation of the blade is divided into N β = 30 angular positions (resolution of 12 o ). The Two test cases are considered to evaluate the accuracy of the Amiet-PE and MM models.

In the first case, only trailing edge noise is included, and the wind speed profile is assumed to be constant in the source model (no wind shear). The propagation conditions are assumed to be homogeneous (c(z) = c 0 ), with a finite impedance ground. The absence of refraction effects makes it possible to compare the results of the coupled model with the analytical solution in a homogeneous atmosphere based on Equations ( 12) and ( 13). In the second test-case, both trailing edge and turbulent inflow noise sources are considered, and the atmosphere is supposed to be neutral. Using the Monin-Obukhov similarity theory, this means that the vertical profiles of the mean wind speed U (z) and of the temperature T (z) are given by [START_REF] Salomons | Computational Atmospheric Acoustics[END_REF]:

U (z) = u * κ ln z z 0 , (16)
T (z) = T 0 + α 0 z, (17)
where u * = 0.49 m/s is the friction velocity, z 0 = 0.1 m is the surface roughness length, III).

T 0 =
In the Amiet-PE model, the initial starter is computed at a distance x is = 100 m, and the far-field receiver is placed at x R = 1000 m; see Fig. order of magnitude of the computation time, a set of PE calculations for the 49 frequencies takes approximately 8 minutes to run on one core of a PC equipped with an Intel Xeon X5650 processor at 2.66 GHz. For each direction τ , the computation time of the MM model is thus between approximately 24 minutes with N h = 3 and 2h30' with N h = 19, and the computation time of the Amiet-PE model is greater than 30 hours.

B. Validation in a homogeneous atmosphere

First, the third octave band spectra of SPL averaged over one rotation are plotted in Fig. 5 for a receiver at a height of 2 m and at a distance of 500 m or 1000 m downwind (τ = 0 o ). The results with the point source approximation, the Amiet-PE model, and the MM model with three and seven source heights are compared to the analytical solution.

Using the point source approximation, there are fluctuations due to ground interference dips that are much reduced using an extended source model. These fluctuations are still visible in the MM model with 3 source heights. Excellent agreement is found with the analytical solution using either the Amiet-PE model or the MM model with 7 source heights. Then, the overall sound pressure level (OASPL) averaged over one rotation and the amplitude modulation (AM) are plotted as a function of x in Fig. 6 for a receiver at a height of 2 m crosswind (τ = 90 o ). The amplitude modulation is defined as the difference between the maximum and the minimum of the OASPL over one rotation. On the one hand, all the OASPL predictions are within 1 dB(A) from the analytical calculation, even with the point source approximation. On the other hand, AM is seen to be much more sensitive to the source model used. As explained in Sec. II D, the point source approximation is only used to account for propagation effects in the present study, so it is theoretically possible to calculate AM using this source model, although the predicted value does not agree with the analytical solution. Using the Amiet-PE model, the AM predictions are accurate only for distances larger than 500 m approximately, which can be attributed to the fact that receivers at short ranges are far from the point at x R = 1000 m where all the propagation planes cross, Finally, the directivities of OASPL and AM at a distance 1000 m and a height of 2 m are plotted in Fig. 7, using an angular step ∆τ = 10 o . The Amiet-PE calculations are not shown because it would be too computationally expensive with such a small value of ∆τ . In the OASPL directivity plot, some differences compared to the analytical solution are obtained with the point source approximation, that remain smaller than 1 dB(A) except close to the interference dips, while very accurate results are obtained using the MM model with 7 heights. The AM predictions using the point source approximation are completely off, while the ones obtained with the MM model are quite accurate, especially with 7 heights.

The MM model predictions with 10 and 19 heights are not shown are they are almost identical to the predictions with 7 heights.

C. Results in a neutrally stratified atmosphere

In a neutral atmosphere, wind turbine noise propagation is completely different downwind and upwind, due to the presence of a shadow zone in the latter case. This is clearly seen To quantify the error made using various methods, let us define the maximum difference over a quantity A(x) as:

MaxDiff(A) = max 500 m≤x≤1200 m |A(x) -A ref (x)| , (18)
where A ref is a reference calculation. In Table I, the maximum difference MaxDiff(OASPL)

is given using various models considering the Amiet-PE model as a reference. Note that the maximum difference is only calculated between 500 m and 1200 m in Eq. ( 18) as the validity of the Amiet-PE model is questionable at short ranges, as discussed in Sec. III B, and because the dwellings are generally located at least 500 m from the closest wind turbine.

Table I shows that the point source approximation yields maximum differences greater than 1.0 dB in the crosswind and upwind directions. The MM model yields accurate results in all directions if at least 10 source heights are considered.

The evolution of AM with distance is plotted in Fig. 8(b) for the same four directions.

The AM remains smaller than 0.2 dB(A) downwind, with similar results for all models. In the other directions, the AM is much higher, and is not well predicted using the point source approximation. The high values of the AM in upward-refracting conditions (τ = 110 o and τ = 180 o ) for distances greater than 400 m are attributed to the fact that the receiver will enter and leave the illuminated region during the blade rotation when it is close to the limit of the shadow zone, as shown in Barlas et al. (2017) andCotté (2018). Since the shadow zone starts at a shorter range for higher frequency, the highest AM is encountered at different frequencies depending on the receiver positions (Cotté, 2018). In Table II, the maximum difference MaxDiff(AM) is given using the Amiet-PE model as a reference. It can be seen Color online.

that at least 10 source heights are needed in the MM model in order to obtain a maximum difference smaller than 1.1 dB in all directions.

To confirm that the MM model predictions converge with increasing source heights in all propagation directions, the directivities of OASPL and AM are plotted in Fig. 9 at a distance of 1000 m and a height of 2 m. Using the MM model with 19 source heights as the reference calculation, the difference between the OASPL predictions are observed in the approximately, refraction effects are small and the OASPL horizontal directivity keeps a dipole shape, as classically measured and predicted at short range (Buck et al., 2016;Oer-lemans and Schepers, 2009;Zhu et al., 2005). At longer ranges, the OASPL directivity changes to an asymmetric shape, with small levels upwind, as seen for instance in Barlas et al. (2017) and [START_REF] Mcbride | A comprehensive hamiltonian ray tracing technique for wind turbine noise propagation under arbitrary weather conditions[END_REF]. The AM directivity shows some peaks in various upwind directions depending on the propagation distance for x > 400 m. This can be attributed to the influence of the acoustic shadow zone, as explained previously. Note that significant AM values have also been reported in the downwind directions in other studies.

IV. CONCLUSION

In this study, two methods have been tested to include extended aeroacoustic source models in a parabolic equation code for wind turbine noise propagation in an inhomogeneous atmosphere. These two methods have been compared to the point source approximation that is classically used in wind turbine noise propagation studies. The source model is based on Amiet's theory, and the parabolic equation code uses a split-step Padé approximant. In the first method, called Amiet-PE, an initial starter is obtained for each segment of the blade using the backpropagation approach. This method enables one to accurately model the directivity of the noise sources but is very computationally intensive. In the second method, the blade segments are viewed as moving monopole sources (MM model), and only a limited number of parabolic equation simulations are needed depending on the number of source heights considered to discretize the rotor plane.

The various models are first validated using an analytical reference solution in a homogeneous medium. The Amiet-PE model and the MM model with at least 7 source heights (∆H ≤ 15 m) are in excellent agreement with the reference solution, in terms of spectra, OASPL and AM. The point source approximation is relatively accurate to predict the OASPL, but it is unable to predict the AM, and tends to exaggerate the ground interference dips in the spectra, even at large distances from the source.

The models are then compared in a neutrally stratified atmosphere, characterized by a logarithmic velocity profile. The most challenging propagation conditions are encountered upwind, where an acoustic shadow zone appears for propagation distances greater than approximately 400 m. The point source approximation fails to calculate the correct OASPL in these directions, because it predicts a shadow zone that starts too close to the wind turbine. In order to correctly capture the AM behavior upwind, the Amiet-PE model and the MM model with at least 10 source heights (∆H ≤ 10 m) are shown to yield accurate results. The MM model is much more computationally effective than the Amiet-PE model, with a ratio M N β /N h ≈ 24 between the two models with N h = 10 source heights.

The MM model proposed in this article could be used in the future to study the effect of strong wind speed gradients, that is potentially the source of amplitude modulation at night [START_REF] Van Den Berg | Wind turbine power and sound in relation to atmospheric stability[END_REF]Zajamsek et al., 2016), the influence of the wind turbine wake on propagation (Barlas et al., 2017), or the combined influence of topography and meteorology, using for instance the rotated PE approach described in [START_REF] Lihoreau | Outdoor sound propagation modeling in realistic environments: Application of coupled parabolic and atmospheric models[END_REF].

  (a). Each segment is represented as an airfoil of chord c m and span L m , m = 1, .., M . The different segments are supposed to be uncorrelated, so that the noise contributions from all blade segments can be summed at the receiver. This assumption is one of the main limitation of strip theory at low frequencies, as discussed byChristophe et al. (2009).

  FIG. 1. Notations for (a) the rotor plane with blades divided into segments, and (b) the wind

  FIG. 2. Schematics for (a) the backpropagation method used to obtain the starter at x = 0 from

Fig. 3

 3 Fig. 3 with M = 6 segments per blade. Using N h = 5 source heights in the MM model, the

  FIG. 4. Vertical profiles of wind speed U (z) and temperature T (z) in a neutral atmosphere. The

  10 o C is the ground temperature, α 0 = -0.01 K/m is the dry adiabatic lapse rate, and κ = 0.41 is the von Kármán constant. The value of the friction velocity is chosen so that U (z = 80 m) = 8 m/s. The vertical profiles of U (z) and T (z) are plotted in Fig.4. In both test-cases, the scattering effect of turbulence is not included in the model, which means that the SPL might be underestimated when a shadow zone is present (Cotté, 2018, Section 4.4).The propagation domain has a size of 1200 m along x and 300 m along z. PE calculations are performed for 49 frequencies in order to predict the third octave band spectra between 100 Hz and 2000 Hz(Cotté, 2018). The ground impedance is calculated with a two-parameter variable porosity model, which is physically admissible and yields a better agreement with measurements than commonly used one-parameter models (e.g. Delany-Bazley or Miki), as shown by[START_REF] Dragna | On the inadvisability of using single parameter impedance models for representing the acoustical properties of ground surfaces[END_REF]. The effective resistivity is σ e = 50 kNs/m 4 and the rate of change of the porosity is α e = 100 m -1 , that are typical values for a natural soil(Dragna et al., 2015, Table

  FIG. 5. Third octave band spectrum of the SPL downwind (τ = 0 o ) at z = 2 m and (a) x = 500 m or (b) x = 1000 m: analytical solution ( ), point source approximation (--), Amiet-PE ( ), MM with 3 heights ( ) or 7 heights ( ). Color online.

  FIG. 6. OASPL and AM versus x at z = 2 m crosswind (τ = 90 o ): analytical solution ( ), point

  FIG. 7. Directivity of (a) OASPL and (b) AM at x = 1000 m and z = 2 m: analytical solution ( ), point source approximation (--), MM with 3 heights ( ) or 7 heights ( ). Color online.

  FIG. 8. OASPL and AM with respect to x at z = 2 m in a neutral atmosphere at τ = 0 o (downwind), τ = 80 o , τ = 110 o , and τ = 180 o (upwind): Amiet-PE (-), point source approximation (--), MM with 19 heights (•). Color online.

  FIG. 9. Directivity of (a) OASPL and (b) AM at x = 1000 m and z = 2 m in a neutral atmosphere: MM with 19 heights ( ), MM with 10 heights ( ), MM with 3 heights ( ), and point source approximation (--).

  upwind directions (100 o ≤ τ ≤ 260 o ), with differences up to 8.3 dB(A) for the point source approximation, 2.4 dB(A) for the MM model with 3 heights, and only 0.3 dB(A) for the MM model with 10 heights. The same behavior is observed in the AM directivities, with differences up to 4.6 dB(A) for the point source approximation, 3.1 dB(A) for the MM model with 3 heights, and only 0.5 dB(A) for the MM model with 10 heights. The movie Mm. 1 shows how the OASPL and AM horizontal directivities vary for distances between 200 m and 1200 m every 10 m. The relative contributions of trailing edge noise and turbulent inflow noise are also plotted, as can be seen in Fig. 10 where two snapshots of the movie corresponding to x = 300 m and x = 1000 m are shown. Up to 300 m

For

  FIG. 10. Directivity of OASPL and AM calculated with the MM model using N h = 10 source

  FIG. 11. Third octave band spectrum of the SPL in a neutral atmosphere at x = 1000 m and

TABLE I .

 I Maximum difference MaxDiff(OASPL) with respect to the Amiet-PE model using the

	point source (PS) approximation and the MM model with 3 heights (MM3), 7 heights (MM7),
	10 heights (MM10) and 19 heights (MM19) for different angles of propagation τ . Boldface values
	correspond to differences strictly greater than 1.0 dB.
	τ	PS	MM3	MM7	MM10	MM19
	0 o	0.6 dB	0.4 dB	0.3 dB	0.3 dB	0.3 dB
	80 o	1.1 dB	1.0 dB	0.9 dB	0.8 dB	0.9 dB
	110 o	5.6 dB	3.6 dB	1.2 dB	0.8 dB	0.8 dB
	180 o	11.0 dB	2.6 dB	0.4 dB	0.4 dB	0.3 dB

TABLE II .

 II Maximum difference MaxDiff(AM) with respect to the Amiet-PE model using the point

	source (PS) approximation and the MM model with 3 heights (MM3), 7 heights (MM7), 10 heights
	(MM10) and 19 heights (MM19) for different angles of propagation τ . Boldface values correspond
	to differences strictly greater than 1.0 dB.
	τ	PS	MM3	MM7	MM10	MM19
	0 o	0.1 dB	0.4 dB	0.1 dB	0.0 dB	0.0 dB
	80 o	2.1 dB	0.2 dB	0.3 dB	0.3 dB	0.4 dB
	110 o	3.4 dB	4.0 dB	1.3 dB	1.1 dB	0.5 dB
	180 o	4.7 dB	3.8 dB	1.4 dB	1.0 dB	0.4 dB

ACKNOWLEDGMENTS

The author would like to thank Benoit Gauvreau and David Ecotière from the Environmental Acoustics Unit (UMRAE), and Tommy Rigall from IMSIA for their useful comments on the manuscript. Barlas, E., Zhu, W., Shen, W., Dag, K., and Moriarty, P. (2017). "Consistent modelling of wind turbine noise propagation from source to receiver," Journal of the Acoustical Society of America 142(5), 3297-3310. Zajamsek, B., Hansen, K., Doolan, C., and Hansen, C. (2016). "Characterisation of wind farm infrasound and low-frequency noise," Journal of Sound and Vibration 370, 176-190. Zhu, W., Heilskov, N., Shen, W., and Sørensen, J. (2005). "Modeling of aerodynamically generated noise from wind turbines," Journal of Solar Energy Engineering 127, 517-528.