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...
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
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
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).
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.
+ (𝑠) 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