Probability mapping of iron pan presence in sandy podzols in South-West France, using digital soil mapping

Probability mapping of iron pan presence in sandy podzols in South-West France, using digital soil mapping

Accepted Manuscript Probability mapping of iron pan presence in sandy podzols in South-West France, using digital soil mapping Anne C. Richer-de-Forg...

985KB Sizes 0 Downloads 36 Views

Accepted Manuscript Probability mapping of iron pan presence in sandy podzols in South-West France, using digital soil mapping

Anne C. Richer-de-Forges, Nicolas P.A. Saby, Vera L. Mulder, Bertrand Laroche, Dominique Arrouays PII: DOI: Reference:

S2352-0094(16)30148-1 doi: 10.1016/j.geodrs.2016.12.005 GEODRS 111

To appear in:

Geoderma Regional

Received date: Revised date: Accepted date:

3 June 2016 16 December 2016 22 December 2016

Please cite this article as: Anne C. Richer-de-Forges, Nicolas P.A. Saby, Vera L. Mulder, Bertrand Laroche, Dominique Arrouays , Probability mapping of iron pan presence in sandy podzols in South-West France, using digital soil mapping. The address for the corresponding author was captured as affiliation for all authors. Please check if appropriate. Geodrs(2016), doi: 10.1016/j.geodrs.2016.12.005

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

ACCEPTED MANUSCRIPT Probability mapping of iron pan presence in sandy podzols in South-West France, using Digital Soil Mapping Anne C. Richer-de-Forgesa, Nicolas P.A. Sabya, Vera L. Mulderb, Bertrand Larochea

b

INRA, InfoSol Unit, US 1106, 45075 Orléans, France Soil

Geography

and

Landscape

Group,

Wageningen

RI

a

PT

& Dominique Arrouaysa

University,

SC

Droevendaalsesteeg 3, P.O. Box 47, 6700 AA, Wageningen, The Netherlands

NU

Corresponding author: [email protected]

MA

Highlights

Iron pan presence was successfully mapped using field observations and ancillary

D

data

PT E

Ancillary data provided added value to the quality of maps

Abstract

CE

Possible uses of these probability maps are described

AC

This work evaluated two different digital soil mapping methods for mapping the presence of iron pans in South-West France. The presence of iron pans limit rooting depth, thereby affecting available water content for plants and increasing vulnerability of trees to storms. In some cases, it may also limit the water infiltration rate and cause anaerobic conditions limiting rooting depth, biological activity and plant growth. This work evaluates the potential of a road-side survey sampling and subsequent digital soil mapping techniques to map the probability of iron pan presence in a

1

ACCEPTED MANUSCRIPT French region. We tested 1) indicator kriging without taking into account ancillary covariates and compared the mapped probability with 2) logistic regression kriging which did account for ancillary covariates. The probability of iron pan occurrence was mapped at a regional level with a reasonable precision but adding ancillary covariates to the model gave more realistic predictions. Likely, the goodness of the

PT

prediction is still hampered by the resolution of the employed DEM and by the spatial

RI

distribution of the observation points. Furthermore, we made some suggestions of the

forestry, agriculture and geotechnical projects.

NU

Keywords

SC

possible usage of the iron pan map; probability maps could help regional planning for

MA

Iron pan presence, soil, digital soil mapping, podzols, France

D

Short running title

AC

CE

PT E

Predicting iron pans presence using DSM

2

ACCEPTED MANUSCRIPT 1. Introduction Iron pans are encountered in many soil horizons. They are particularly common within PODZOLS, forming cemented or compacted spodic horizons (IUSS Working Group WRB, 2014). The presence of iron pans limits rooting depth, thus affecting available water content for plants and vulnerability of trees to storms. In some cases,

PT

it may also limit the water infiltration rate and cause anaerobic conditions. The study

RI

was carried out in a region covering more than 1,150, 200 ha in southwestern

SC

France: the region ‘des Landes de Gascogne’. This region is mainly covered by sandy PODZOLS containing an iron pan which presence is generally not continuous.

NU

Iron pans are geographically distributed as little patches in the landscape, which makes their mapping difficult. Moreover, constraints about accessibility strongly limit

MA

the number and the density of field observations that can be performed. Therefore, we tested the potential of a road-side survey sampling and subsequent digital soil

PT E

D

mapping techniques to map the probability of iron pan presence at a regional level.

The main objectives of this study were 1) assessing the feasibility of mapping the

CE

probability of iron pan presence over this region, 2) evaluating the potential improvements of the maps when adding ancillary data to the model and 3) providing

AC

end-users with accurate data, allowing them to make an informed decision in riskassessments related to iron pan presence.

3

ACCEPTED MANUSCRIPT 2. Material and Method 2.1. Study area The study area is located in southwestern France, approximately between 43.3 and 44.5°N and 1.45°W to 0.10° E (Figure 1). The ‘Landes de Gascogne’ is a nearly flat

PT

(slope <= 0.5%) plain with elevation ranging from 0 to 256 m above sea level .

SC

RI



The region receives about 900 mm year−1 of rain, and has a mean annual

NU

temperature of 12.7°C. The land use is mainly composed of even-aged stands of

MA

maritime pine (Pinus pinaster Ait.) some of which having been deforested and converted into intensively managed irrigated maize (Zea Mais)(Jolivet et al., 1997;

D

2007). The parent material of the soils is known as ‘sables des Landes’, which

PT E

originated from tertiary fluvial alluviums that were spread by successive wind erosion periods, during the Pleistocene (Legigan, 1979). The soils that developed from this

CE

parent material are mainly PODZOLS (IUSS Working Group WRB, 2014). Some of these PODZOLS are characterized by the presence of a very hard iron pan,

AC

cemented by organic compounds and Al and Fe oxydes (Cazenave, 1970). The profile development is directly related to the seasonal variations of the water table (Righi and Wilbert, 1984). The spodic horizons may be cemented or not, depending on the water table regime (Cazenave, 1970; Righi and Wilbert, 1984; Jolivet et al., 2007) (Supplementary material 1 and 2). According to these authors, the depth and seasonal variations of the water table are mainly controlled by: i)

the distance and relative height to small valleys,

4

ACCEPTED MANUSCRIPT ii)

an eolian micro-relief mainly consisting of a succession of low ridges and shallow troughs, ranging from about 0.30 to 1.50 m high and 10-50 m wide.

Therefore, the iron pan is generally not continuous and is distributed as little patches in the landscape which makes its mapping difficult.

PT

Furthermore, close to the Atlantic Ocean, the coastal dunes region developed into foredunes and secondary dunes more or less parallel to the coast. The parent

RI

material of the soils originated from eolian deposits during the Holocene period

SC

(Penin, 1980). On these dunes, soils are mainly ARENOSOLS (IUSS Working Group

NU

WRB, 2014), with little or no soil development. However, as these dunes progressively covered the inland formations, the iron pans may be present and buried

MA

at their basis (Arenosols over Podzols), and are sometimes observed in depressions

2.2. Field observations

D

between the lines of dunes.

PT E

Field observations included 4,102 auger borings and 286 soil pits, corresponding to a mean sampling density of 0.38 observation per km 2 (Figure 1).The sampling strategy

CE

mainly entailed a road-side survey, where the presence/absence of an iron pan was registered. Also, the geographical coordinates of each sampling site were recorded

AC

with an average precision of + 5 m. 2.3. Co-variates

Our hypothesis was that topographic variability was the main controlling factor for the seasonal fluctuations of the water table (Cazenave, 1970; Righi and Wilbert, 1984). In addition, we used land cover as a supplementary co-variate, as irrigated agricultural crops had been preferentially located in areas having a water-table at

5

ACCEPTED MANUSCRIPT shallow depths. Altogether, a set of 12 potential covariates were selected for model calibration (Table 1). 2.4. Modelling and mapping the presence of iron pans

PT

We compared two interpolation methods, i.e. a pure geostatistical approach (indicator

RI

kriging) and a hybrid statistical/geostatistical approach (logistic regression kriging).

SC

Previously, both methods proved successful in mapping the presence or absence of a soil characteristic and both approaches allowed the prediction of a categorical

NU

variable, being either presence (1) or absence (0). Here, we describe only briefly the two methods, see Hengl et al. (2003; 2007?) for more details. For both methods, we

MA

started by transforming the categorical variable presence/absence at sample location

1, 𝐼𝑟𝑜𝑛 𝑃𝑎𝑛 𝑖𝑠 𝑝𝑟𝑒𝑠𝑒𝑛𝑡 𝑎𝑡 𝑙𝑜𝑐𝑎𝑡𝑖𝑜𝑛 𝑠 (1) 0, 𝑜𝑡ℎ𝑒𝑟𝑤𝑖𝑠𝑒

PT E

+ 𝑧𝐼𝑟𝑜𝑛 𝑃𝑎𝑛 (𝑠) = {

D

s into a 0/1 indicator (Bierkens and Burrough, 1993):

CE

2.4.1. A pure geostatistical technique: the indicator kriging approach The spatial distribution of the probability of iron pan presence was modelled following

AC

an indicator kriging approach (Goovaerts, 2009). Indicator kriging provides a nonparametric distribution, estimated directly using fixed thresholds by considering indicator transforms of the conditioning data in the form of cumulative distribution functions with step functions. The method estimates at the unsampled location 𝑢 the probability of iron pan presence by applying kriging to indicator transform (Goovaerts, 2009).

6

ACCEPTED MANUSCRIPT 𝑛 + 𝑧𝐼𝑟𝑜𝑛 𝑃𝑎𝑛 (𝑢)

+ = ∑ 𝜆𝛼 𝑧𝐼𝑟𝑜𝑛 𝑃𝑎𝑛 (𝑠) 𝛼=1

2.4.2. Hybrid statistical/geostatistical approach: logistic regression kriging

PT

Regression kriging allows to take into account ancillary data. With regression kriging,

RI

typically, there is a fixed trend component estimated by predictor variables and a

SC

residual component which is subsequently kriged after the trend component has been removed from the target variable (Hengl et al., 2007). In order to model the

NU

fixed trend present in the data of iron pan presence, a logistic regression model was fit using a logit model (Hengl et al., 2003). This means that the regression modelling

MA

is supplemented with the modelling of the variogram of the residuals, which can then be interpolated and added back to the regression.

PT E

D

The prediction can be done using:

+ + 𝑇 −1 𝑧𝐼𝑟𝑜𝑛 + 𝑒𝐼𝑟𝑜𝑛 𝑃𝑎𝑛 (𝑠) = [1 + exp(−𝛽 . 𝑞0 )] 𝑃𝑎𝑛 (𝑠) (2)

CE

+ (𝑠) the Where are 𝛽 the regression coefficients 𝑞0 are the predictors and 𝑒𝐴𝑙𝑖𝑜𝑠

interpolated residuals. This two-step approach was implemented using R (R

AC

development core team, 2013) and the contributing GSIF R package (Hengl, 2014). We refer to this method in the following text as “RKI”.

2.5. Model calibration and evaluation The strongly clustered points, located along the roads, were not ideal for the model calibration and validation. In order to overcome this issue, we applied an evaluation

7

ACCEPTED MANUSCRIPT exercise (Oreskes et al., 1994; Oreskes, 1998) instead of a formal validation. We developed a specific procedure, based on a stratified scheme, to take into account this configuration in an evaluation procedure which resembles a k-fold validation. Hence, we first divided the region into 314 regular squares (5 x 5 km) and selected randomly one sampled point in each of these squares which was kept aside for

PT

evaluation. Then, we ran the calibration on all the remaining observations. In a

RI

second step, the probability iron pan presence was predicted at the location of the

SC

evaluation sites and for these sites the evaluation criteria (section 2.6) were calculated. We ran this stratification 50 times without duplicate use of the samples

MA

NU

sites and averaged the evaluation criteria.

2.6. Assessment of the quality of predictions

D

We describe the accuracy measures for iron pan presence maps only briefly. For an

PT E

elaborate review, including the estimation of these measures, we refer to Brus et al. (2011) and de Gruijter et al. (2006).

CE

We first considered three map quality measures: overall purity, the purity defined at the level of the two map units (presence referred as p1 and absence referred as p0)

AC

and the class representation for the two map units (presence referred as r1 and absence referred as r0) (Brus et al., 2011). Overall purity is defined as the proportion of the mapped area in which the predicted presence of iron pan, equals the true value. In other words, it is the areal proportion of the map unit correctly classified. The class representation for the soil class u is the proportion of the area where in reality soil class u occurs that is also mapped as u. We also considered the

8

ACCEPTED MANUSCRIPT percentage pc of evaluation sites positively predicted (i.e., presence of iron pans). Finally, we calculated the kappa and Naesset’s index of the maps (Naesset, 1996).

2.7. Spatial predictions for mapping

PT

The final DSMs of iron pan presence are the averaged outcome of the 50 predictions of the probability of iron pan presence, on a 250 meter grid, using either IK or RKI. In

RI

addition, we produced several indicator (binary) maps representing a reclassification

SC

of the final DSM by applying a threshold value on the probability maps (i.e. 0.5 and

NU

0.8). The latter were considered important for specific risk assessments. The 0.5 probability was chosen as a threshold for iron pan to be more probable than

MA

not. The 0.8 threshold shows the reduction in the area and points out those areas where iron pan presence is largely dominant and where the energy and the cost

PT E

D

needed to break this iron pan will be the highest.

AC

3.1. Evaluation

CE

3. Results

The overall accuracy diagnostics were very close among the two models, ranging from about 0.70 to about 0.74 (Figure 2). For IK, it shows a slight decrease when the threshold is set to 0.8 whereas it remains stable for RKI. This is likely attributable to the fact that IK predictions become more dependent on sampling density when the threshold increases. The map unit purity of the absence of iron pans is very close for both models and it decreases when the threshold increases to 0.8. As this last threshold is stricter on the presence of iron pans, it seems logical that we observe

9

ACCEPTED MANUSCRIPT more iron pans in areas predicted as iron pan absence. The map purity of the presence of iron pan (p1) is rather similar (close to 0.6) for IK and RKI when the threshold is fixed to 0.5. It increases when the threshold is fixed to 0.8. It is noticeable that the median map purity of the presence of iron pan (p1) using IK and the 0.8 threshold is close to 1. It means that the areas predicted as having iron pans are

PT

quite pure and likely due to a high density of observed iron pans in the sampling. The

RI

proportion of validation sites positively predicted (pc), i.e. as showing iron pans, falls

SC

from more than 20% to less than 5% when the threshold increases from 0.5 to 0.8. The relative difference between IK and RKI increases with the threshold. It is not

NU

surprising as RKI, using the relation with co-variates, is able to predict iron pans presence where there are few sampling points whereas IK is not. The class

MA

representation of the absence of iron pan (r0) is rather similar for both methods (IK and RKI). It shows a marked increase when the threshold is fixed to 0.8, exhibiting

D

values close to 1. These very high values of class representation are likely due to the

PT E

fact that the area without iron pan is much larger than the area showing iron pans. The class representation of the presence of iron pans is much lower and decreases

CE

when the threshold is set to 0.8.

AC



The results from kappa and tau indexes indicate a rather low performance of predictions (Figure 3). They both decrease when the threshold increases from 0.5 to 0.8. The difference between tau and kappa, especially for the threshold 0.8, can be explained by the fact that most of the study area is dominated by soils without iron pan.


10

ACCEPTED MANUSCRIPT 3.2. Maps In this section we show the maps of the probability of iron pan presence using the two methods. The map produced using RKI (Figure 4) looks contrasted, exhibiting very large areas

PT

where the probability is zero or very low and patches were the probability is very close to 1. Iron pan is mainly predicted to occur in rather flat, interfluvial areas close

RI

to the western coastal lakes and in flat, eastern areas located far from the valleys.

SC

Note that iron pan is also predicted to occur on the northern part of the coastal line and in some small patches in the coastal dunes. In the southern part of the region,

NU

characterized by a denser valley network (see fig. 1), most of the probabilities are

MA

close to zero.


D

The map produced using IK (Figure 4) is smoother, even though some patches close

PT E

to the lakes are still visible and the global patterns of both maps are consistent. However, the coastal dunes exhibit a smooth gradient of iron pan presence from

CE

north to south. This may be attributable to two reasons: i) the number of observation points was very limited in the coastal dunes area, especially in its northern part (see

AC

Fig. 1) and ii) in the adjacent flat plain interfluves, the iron pan is more often observed in the northern part, which may have influenced the kriging predictions. The map of the difference between RKI and IK is shown Figure 4. Large negative differences are present along the small rivers network and on the northern part of the coastal dunes. The pattern of positive differences is smoother and mainly related to the high relative position in the landscape.

11

ACCEPTED MANUSCRIPT Figure 5 shows the areas where the probability of iron pan presence exceeds a given probability value. We observe a very strong decrease in iron pan presence when the probability is set to 0.8. Most of the large patches with high probability are located in interfluvial areas close to the coastal lakes. Smaller patches are more numerous with

PT

RKI, especially along the coastal area.

RI



SC

4. Discussion

NU

4.1. Which maps make sense?

The maps clearly show that there are differences between the employed models.

MA

Indeed, the smoothed pattern shown by the IK map is not reliable, as judged by our knowledge on the drivers of iron pan presence. The RKI map looks much more

D

realistic and is able to capture some sharp discontinuities and local variations. This is

PT E

particularly visible in the coastal area. In this region, iron pans were successfully predicted in the depressions within the costal dunes and even along some beaches.

CE

Indeed, before having been fixed by pine plantation in the 19 th century, the dunes continuously moved toward inland, burying progressively the PODZOLS of the plain.

AC

At the same time, their western edge was eroded by the ocean, so that paleo-podzols (Arenosols under Podzols) now appear between dunes alignments and at the basis of the costal line (supplementary material 3). The pattern of the RKI predictions is closely related to the fixed part of the model which corresponds to the distance to and relative height of small valleys, which was an expected result. This is especially visible in the southern part of the map, where the absence of iron pan was expected due to the rather dense network of valleys

12

ACCEPTED MANUSCRIPT (see fig. 1). Moreover, the RKI model captures more local variability and better reflects the fact that the proportion of soils having an iron pan is highly variable in space (Righi and Wilbert, 1984). The map of the differences between RKI and IK maps reinforces the conclusions that RKI allowed capturing some expected patterns of iron pan absence (i.e., along the valleys), whereas IK was not. We conclude that

PT

the smooth pattern of the IK map is not realistic, given the high local variability of the

RI

controlling factors of iron pan presence in this region (Cazenave, 1970; Righi and

SC

Wilbert, 1984; Jolivet et al., 2007). Moreover, due to the small sampling density of observations in the dunes, IK predicted smooth areas of iron presence in the

NU

northern coastal dunes, whereas RKI was able to locate them more accurately. Finally, both the evaluation results using 50 validation sets and our expert judgement

MA

lead to the conclusion that the RKI predictions are more reliable than the IK ones. 4.2. Limitations

PT E

D

The overall accuracy was comparable and consistent with previous studies (e.g., Collard et al., 2014; Grinand et al., 2008; Kempen et al., 2009, 2012; Lemercier et al., 2012), and in accordance to the recommendation by Marsman and de Gruijter

CE

(1986), who suggested 70% as an acceptable value for the producer. However we

AC

must keep in mind that our mapping exercise includes only two mapping units. The class representation of iron pan suggests that our model was not fully efficient to accurately map the iron pan presence over the entire area, even with a very large number of observations. The kappa and tau index were rather low, indicating a rather poor performance. This suggests that the employed co-variates did not allow capturing all the soil-landscape relationships determining the iron pan presence. This may partly be attributed to the fact that the resolution of the employed DEM was not precise enough to capture very local variations in the micro-relief. Indeed, further

13

ACCEPTED MANUSCRIPT trials should be performed when a higher resolution DEM will become available. Moreover, maps of the water table fluctuations could be useful to incorporate as covariates. Unfortunately, such maps are not yet available with enough accuracy. In addition, increasing the field sampling density in coastal dunes and in areas located

PT

far from the roads should improve the predictions. However, our results did demonstrate that gathering numerous field and relevant

RI

ancillary data helps providing a regional prediction of iron pan presence that makes

SC

sense, especially in places where their presence is very dense. On the other hand, our mapping exercise may suffer from limitations related to the spatial coverage of

NU

the dataset. Indeed, in the less accessible areas where less intensive sampling was

MA

conducted, we may have missed some patterns of iron pan presence or absence. Further work should be conducted, including the use of new potentially relevant covariates such as a more precise DEM and understory vegetation. A more precise

PT E

D

DEM may better capture local elevation variability, and the understory vegetation maps could better reflect the water-regime (i.e., Jolivet et al., 2007). We should also test other modelling approaches to improve the prediction in this area, for instance by

CE

better taking into account the clustered spatial structure of the sampling.

AC

4.3. Usefulness and prospects Iron pan is a strong limiting factor for rooting depth of trees and crops and it is very hard to break. Therefore, from a very practical point of view, the produced iron pan map will help regional planning for forestry and agriculture (e.g. choosing tree species for plantation, estimating irrigation needs of summer crops). Deriving some maps of the probability of exceeding a given probability of iron pan presence (see figure 5) may also have practical usefulness for geotechnical projects. Such

14

ACCEPTED MANUSCRIPT probability threshold maps may be useful to evaluate the efforts, the machinery and the related costs to break the iron pan if necessary.

The probability of iron pan presence was mapped at a regional level, with a

PT

reasonable precision. A next step could be to try to map its upper depth in order to produce a spatial prediction of the effective rooting depth, as recommended by the

RI

GlobalSoilMap project (Arrouays et al., 2014). Indeed, we wanted first to know where

SC

iron pans were present because there is no reason to try to map the depth of a layer

NU

that does not exist. Mapping effective depth will be a second step which is not trivial and can be done using a large range of mapping methods. A much more challenging

MA

step would be to map iron pan thickness, which could be helpful for geotechnical applications such as burying electrical lines or installing pipe lines or for estimating

D

the stocks of mineral and organic elements (mainly C, Fe and Al) contained in this

PT E

deep horizon.

Further work on the effect of the road-side sampling design should also be conducted

AC

calibration.

CE

in order to better take into account the effect of this kind of sampling on the model

5. Conclusion

The probability of iron pan presence in a region of South-West France was mapped using a large number of field data, most of them having been collected using a roadside survey sampling strategy. For this purpose, we tested two digital soil mapping methods. One is using only point data, and the other one uses the same point data plus a set of ancillary variates. Adding the ancillary variates resulted in an

15

ACCEPTED MANUSCRIPT improvement of the predictions, despite the rather high spatial density of our point data. The modelled iron pan presence in combination with a risk-assessment is considered very practical for various end-users. The accuracy of the prediction is currently mainly hampered by the resolution of the employed DEM and by the spatial distribution of the observations points. We believe that a higher resolution DEM (both

PT

in vertical and horizontal resolution) will help to improve the prediction accuracy in the

AC

CE

PT E

D

MA

NU

SC

RI

near future.

16

ACCEPTED MANUSCRIPT References Arrouays, D., Grundy, M.G., Hartemink, A.E., Hempel, J.W., Heuvelink, G.B.M., Hong, S.Y., Lagacherie, P., Lelyk, G., McBratney, A.B., McKenzie, N.J., Mendonça-Santos, M.D.L, Minasny, B., Montanarella, L., Odeh, I.O.A., Sanchez,

PT

P.A., Thompson, J.A., Zhang, G.-L., 2014. GlobalSoilMap: towards a fine-

RI

resolution global grid of soil properties. Advances in Agronomy 125, 93-134. Bierkens M.F.P., Burrough PA. 1993. The indicator approach to categorical data. 1.

SC

Theory. Journal of Soil Science 44(2), 361-368.

NU

Brus, D.J., Kempen, B., Heuvelink, G.B.M., 2011. Sampling for validation of digital

MA

soil maps. Eur J Soil Sci. 62(3), 394-407.

Casenave, A., 1970. Contribution à l’étude de l’alios. Université de Bordeaux. 196 p.

D

(PhD Thesis).

PT E

Collard F., Kempen B., Heuvelink G. B. M., Saby N. P. A., Richer-de-Forges A. C., Lehman S., Nehlig P. and Arrouays D. 2014. Refining a reconnaissance soil map

CE

by calibrating regression models with data from the same map (Normandy,

AC

France). Geoderma Regional 1, 21-30. de Gruijter, J., Brus, D.J., Bierkens, M.F.P., Knotters, M. 2006. Sampling for natural resources monitoring. Springer-Verlag, Berlin Heidelberg, 334 p. Goovaerts P. 2009. AUTO-IK: a 2D indicator kriging program for the automated nonparametric modeling of local uncertainty in earth sciences. Comput Geosci. 35(6), 1255-1270.

17

ACCEPTED MANUSCRIPT Grinand, C., Arrouays, D., Laroche, B., Marin, M.P., 2008. Extrapolating regional soil landscapes from an existing soil map : sampling intensity, validation procedures, and integration of spatial context. Geoderma 143, 180-190. Hengl, T., 2009. A Practical Guide to Geostatistical Mapping, 2nd Edt. University of

PT

Amsterdam, www.lulu.com, 291 p.

SC

2/r142. http://R-Forge.R-project.org/projects/gsif/

RI

Hengl, T., 2014. GSIF: Global Soil Information Facilities. R package version 0.4-

IUSS Working Group WRB., 2014. World Reference Base for Soil Resources 2014.

NU

International soil classification system for naming soil and creating legends for soil

MA

maps. World Soil Resources Reports N° 106. FAO, Rome. Jolivet, C., Arrouays, D., Andreux, F., Lévèque, J., 1997. Soil organic carbon

D

dynamics in forested spodosols converted to maize cropping. Plant and Soil 191,

PT E

225-231.

Jolivet, C., Augusto, L., Trichet, P., Arrouays, D., 2007. Les sols du massif forestier

CE

des Landes de Gascogne : formation, histoire, propriétés et variabilité spatiale. Revue Forestière Française 59(1), 7-30.

AC

Kempen, B., Brus, D.J., Heuvelink, G.B.M., Stoorvogel, J.J., 2009. Updating the 1:50,000 Dutch soil map using legacy soil data: a multinomial logistic regression approach. Geoderma 151(3-4), 311-326. Kempen, B., Brus, D.J., Stoorvogel, J.J., Heuvelink, G.B.M., De Vries, F., 2012. Efficiency comparison of conventional and digital soil mapping for updationg soil maps. Soil Sci. Am. J. 76(6), 2097-2115.

18

ACCEPTED MANUSCRIPT Legigan, P., 1979. L’élaboration de la formation du sable des Landes, dépôt résiduel de l’environnement sédimentaire pliocène-pléistocène centre aquitain. PhD thesis, University of Bordeaux, France. Lemercier, B., Lacoste, M., Loum, M., Walter, C., 2012. Extrapolation at regional

PT

scale of local soil knowledge using boosted classification trees: a two-step approach. Geoderma 171-172, 75-84.

RI

Marsman, B.A., de Gruijter, J.J., 1986. Quality of soil maps. A comparison of survey

SC

methods in a sandy area. Soil Survey Papers, 15. Netherlands Soil Survey

NU

Institute, Wageningen.

Naesset, E., 1996. Conditional tau coefficient for assessment of producer's accuracy

D

Remote Sensing, 51(2), 91-98.

MA

of classified remotely sensed data. ISPRS Journal of Photogrammetry and

PT E

Oreskes, N., Shrader-Frechette, K., Belitz, K., 1994. Verification, validation, and confirmation of numerical models in the earth sciences. Science, 263, 641-646.

CE

http://doi.org/10.1126/science.263.5147.641. Oreskes, N., 1998. Evaluation (not validation) of quantitative models. Environmental

AC

Health Perspectives, 106(Suppl 6), 1453-1460. Penin, F., 1980. Le prisme littoral aquitain: histoire holocène et évolution récente des environnements morphosédimentaires. PhD Thesis, University of Bordeaux I, France. n°1577, 129 p. R Development core team, 2013. R: A language for statistical computing. R Foundation for Statistical Computing, Vienna,
19

ACCEPTED MANUSCRIPT Righi, D., Wilbert, J., 1984. Les sols sableux podzolisés des Landes de Gascogne

AC

CE

PT E

D

MA

NU

SC

RI

PT

(France): répartition et caractères principaux. Science du Sol 4, 253–264.

20

AC

CE

PT E

D

MA

NU

SC

RI

PT

ACCEPTED MANUSCRIPT

Figure 1

21

AC

CE

PT E

D

MA

NU

SC

RI

PT

ACCEPTED MANUSCRIPT

Figure 2 22

NU

SC

RI

PT

ACCEPTED MANUSCRIPT

AC

CE

PT E

D

MA

Figure 3

23

AC

Figure 4

CE

PT E

D

MA

NU

SC

RI

PT

ACCEPTED MANUSCRIPT

24

AC

Figure 5

CE

PT E

D

MA

NU

SC

RI

PT

ACCEPTED MANUSCRIPT

25

ACCEPTED MANUSCRIPT

AC

CE

PT E

D

MA

NU

SC

RI

PT

Table 1: list of co-variates used for regression-kriging of indicators Co-variate Classical relief attributes, as derived from a 25m resolution Digital Elevation Model: alt Elevation (m) slope Local slope angle (%) curvm Global curvature curv-long Horizontal curvature curv-trans Profile curvature dppr Relative hydrological distance to the nearest river (m) dppr3_3 mean values of dppr on 3 × 3 pixel windows dppr7_7 mean values of dppr on 7 × 7 pixel windows dppr21_21 mean values of dppr on 21 × 21 pixel windows dppr41_41 mean values of dppr on 41 × 41 pixel windows hppr Relative height to the nearest river (m) Land cover clc Corine Land Cover class

26