3D mapping of soil texture in Scotland

    3D mapping of soil texture in Scotland. Laura Poggio, Alessandro Gimona PII: DOI: Reference:

S2352-0094(16)30128-6 doi: 10.1016/j.geodrs.2016.11.003 GEODRS 104

To appear in: Received date: Revised date: Accepted date:

1 June 2016 2 November 2016 21 November 2016

Please cite this article as: Poggio, Laura, Gimona, Alessandro, 3D mapping of soil texture in Scotland., (2016), doi: 10.1016/j.geodrs.2016.11.003

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.

T

ACCEPTED MANUSCRIPT

RI P

3D mapping of soil texture in Scotland. Laura Poggio∗ 1 , Alessandro Gimona1

The James Hutton Institute - Craigiebuckler, AB158QH, Aberdeen, Scotland (UK)

NU

1

[email protected]

SC

∗

Abstract

MA

Reliable spatially explicit information about soil is important for global environmental challenges. Soil texture is one of the soil most important characteristics as it drives several physical, chemical, biological, and hydrological properties and processes. Despite the importance, there is scarcity

ED

of information on soil texture, especially at the resolution required for environmental modelling. Many recent efforts modelled soil texture with different approaches focussing on the spatial

PT

relationships with environmental covariates. This study aimed at i) modelling and map soil particle classes for Scotland at medium resolution (250m), for topsoil and the whole profile,

AC CE

using an operational DSM approach following specifications from the GlobalSoilMap project; ii) assessing the spatial uncertainty of the modelling approach, and iii) evaluating the impact of spatial and modelling uncertainty on soil texture classification of the topsoil. An extension of the scorpan-kriging approach, i.e. hybrid geostatistical Generalized Additive Models (GAMs), combining GAM with Gaussian simulations was used on Additive-Log-Transformed soil particle classes. The R2 calculated with the validation dataset was between 0.55 and 0.60 and the RMSE values were below 13%. The set of covariates used in this study explained about 40% of the variance of the data. The significant covariates included morphological features, vegetation index and information about the phenological season. The results also showed a large percentage of the variability to be spatially structured. The assessment of the uncertainty on the soil texture classification showed variability and class shift. The resulting datasets can be used as input for further modelling in a number of areas, and they are also important for soil functions modelling

1

ACCEPTED MANUSCRIPT

in the context of provision of Ecosystem services. The accompanying uncertainty can be used to

T

provide supporting information for land management choices.

RI P

Key words:

soil properties, geostatistics, compositional modelling, MODIS, uncertainty. Cambisols, Gleysols,

AC CE

PT

ED

MA

NU

SC

Leptosols, Podzols, Stagnisols, Umbrisols.

2

ACCEPTED MANUSCRIPT

1

Introduction

T

Reliable spatially explicit information about soil is important for global environmental challenges

RI P

such as climate change, food and water shortage, land degradation, and loss of biodiversity (Hartemink and McBratney, 2008).

Soil texture is one of the soil most important characteristics as it drives several physical,

SC

chemical, biological, and hydrological properties and processes, e.g. water and nutrient retention, infiltration, drainage, aeration, SOC content, pH buffering, and porosity. It is therefore

NU

fundamental to understand and model soil functions and provision of ecosystem services. Soil texture is used in soil classification systems. In soil taxonomy it distinguishes soil orders and it is used in the diagnosis of some key epipedons (Bockheim and Hartemink, 2013). Soil texture

MA

also determines the suitability of the soil for a particular use and management and for water management (Thompson et al., 2012). The capacity of soils to maintain organic carbon is influenced by its clay and silt content (Hassink, 1997; Bationo et al., 2007). Particle-size fractions

ED

are inputs in most hydrological, ecological, climatic, and environmental risk assessment models (Liess et al., 2012). The proportions of soil particles have been used to create pedotransfer func-

PT

tions to estimate difficult-to-measure soil properties such as bulk density, hydraulic conductivity, and water holding capacity (Wosten et al., 2001; Minasny and Hartemink, 2011; Toth et al.,

AC CE

2015), important for many plant and soil-water modelling approaches. Soil texture and derived properties are used for soil quality, soil productivity, nutrient dynamics, sensitivity of soils to compaction and soil resilience (Dane et al., 2002; McLauchlan, 2006). The identification of fine and coarse soils allows policy makers to develop soil management techniques. The soil erodibility (K-factor in RUSLE models) largely depends on soil texture (Panagos et al., 2014). Soil texture is also important for the provision of ecosystem services (Adhikari and Hartemink, 2016; Calzolari et al., 2016). Despite the importance, there is scarcity of information on soil texture, especially at the resolution required for environmental modelling. In modelling, quantitative and continuous soil properties rather than classes are required (Arrouays et al., 2014). Soil texture, has conventionally been mapped with polygons where each polygon illustrates a texture class. However, due to the presence of intra-polygon texture variability, there can be a large degree of uncertainty in

3

ACCEPTED MANUSCRIPT

the textural composition within the area labelled by a polygon (Heuvelink and Huisman, 2000). Thus, an alternative way to cope with this problem is to map different textural fractions numer-

T

ically (van Meirvenne and van Cleemput, 2005) so that intra-polygon textural variability can be

RI P

better assessed. Soil properties in soil profiles vary continuously with depth and recent efficient approaches were proposed to map soil properties in 3D space (Malone et al., 2009; Poggio and Gimona, 2014; Orton et al., 2016). Numerous recent methods exist to map soil texture using a

SC

scorpan-kriging approach (McBratney et al., 2003) considering both lateral and vertical variability, e.g. (Adhikari et al., 2013; Akpa et al., 2014; Ballabio et al., 2016; Chagas et al., 2016;

NU

Niang et al., 2014; Gooley et al., 2014; Buchanan et al., 2012). These approaches used a variety of data mining techniques (such as linear models, regression trees, multivariate adaptive regression

MA

splines) and of covariates to describe soil forming factors and management inputs. Most of used covariates were derived from digital elevation model and remote sensing. Some attempts were made with multispectral remote sensing (Ben-Dor, 2002; Mulder et al., 2011), hyper-spectral

ED

remote sensing (Selige et al., 2006; Lagacherie et al., 2008; Ouerghemmi et al., 2016) and radar remote sensing (Niang et al., 2014).

The main aims of this study were:

PT

1. model and map soil particle classes for Scotland at medium resolution (250m) for topsoil and the whole profile using an operational Digital Soil Mapping (DSM) approach following

AC CE

specification from GlobalSoilMap project (Arrouays et al., 2014); 2. assess the spatial uncertainty of the modelling approach, and 3. evaluate the impact of spatial and modelling uncertainty on soil texture classification of the topsoil.

2

Data and test area [Figure 1 about here.] The test area considered was the whole of Scotland (' 78, 000km2 ), excluding the Shetland

islands, providing a wide range of morphological features and soils.

4

ACCEPTED MANUSCRIPT

2.1

Soil data

The Scottish Soils Database contains information and data on soils from locations throughout

T

Scotland. It contains the National Soil Inventory of Scotland (NSIS) profile samples collected

RI P

on a regular 10 km grid of sampled locations (Lilly et al., 2010) and physical and chemical data from a large number of soil profiles taken to characterize the soil mapping units. The properties

SC

necessary were measured to a maximum depth of 1 metre. The positional uncertainty of profile locations and the uncertainty of analytical procedures were not taken into account. About 9000 sampled locations were available for the topsoil model. Of those, 5600 had

NU

recorded granulometry information, while the remaining were classified as organic (fig. 1). In total 26,000 horizons were available for the 3D modelling with 20,255 horizons with recorded

MA

granulometry information. Most of the horizons without granulometry information were classified as organic. In order to provide validation of the models the data available were split in two sets, training and validation, with a ratio of approximately 3:1. The validation set was randomly

ED

sampled, taking into account the geographical coverage and the range of the variables of interest. The two data-sets have similar spatial and values distributions.

PT

The percentage of sand (> 0.05 mm), silt (0.002-0.05 mm) and clay (< 0.002 mm) were modelled using the available data in the Scottish database. The values were defined only for

2.2

AC CE

mineral soils. Therefore pixels with high probability of containing organic soils were masked.

Covariates

The covariates included are freely and globally available and were selected to describe, directly or indirectly, the most important scorpan factors, namely topography, vegetation, climate and geographical position. 2.2.1

Morphology

The Digital Elevation Model (DEM) used as a covariate in the fitted models was SRTM (Shuttle Radar Topography Mission), further processed to fill in no-data voids (Jarvis et al., 2006; Rodriguez et al., 2006). SRTM has a spatial resolution of 90m with global coverage. The measures used were elevation, topographic wetness index and slope as the steepest slope angle, calculated

5

ACCEPTED MANUSCRIPT

using the D8 method (O’Callaghan and Mark, 1984). The DEM was resampled to the resolution used in the analysis (250m). The medians in each

RI P

2.2.2

T

grid cell were used. Remote sensing: MODIS

A set of indices was derived from the Terra Moderate Resolution Imaging Spectro-radiometer

SC

(MODIS) 8 and 16 day composite products. The data were acquired from the NASA ftp website (ftp://e4ftl01u.ecs.nasa.gov/MOLT/) for 12 years between 2000 and 2012. The single images

NU

were restored to fill the cloud gaps (Poggio et al., 2012). The indices were selected for their capability to differentiate spectral responses from different bare soils, vegetation cover and mixed

MA

situations:

1) Enhanced Vegetation Index (EVI; Huete et al., 2002).

ED

2) the Normalised Difference Water Index (NDWI; Gao, 1996).

N DW I =

N IR − SW IR N IR + SW IR

(1)

PT

NDWI was calculated with NIR (Near InfraRed) and SWIR (Short Wave InfraRed) band: SWIR = 2130 (Gu et al., 2008).

AC CE

3) Land Surface Temperature (LST; Wan, 1999) during the day. 4) primary productivity (Running et al., 1999, 2004). 5) Phenology information: Phenology is the response of the vegetation to seasonal climatic cycles in irradiance, temperature and rainfall. Therefore phenology constitutes an essential land surface parameter in atmospheric and climate models. Using fitted functions reduces the uncertainties and leads to more stable measures. Timesat (J¨ onsson and Eklundh, 2002, 2004; Eklundh and Jonsson, 2011) is primarily designed to process time-series of vegetation index derived from satellite spectral measurements. The processing method used was based on least-squares fits to the upper envelope of the vegetation index data with local polynomial functions in the fitting, and the method can be classified as an adaptive Savitzky-Golay filter. The parameters calculated are: (a) length of the season (LOS); time from the start to the end of the season.

6

ACCEPTED MANUSCRIPT

(b) seasonal amplitude; difference between the maximum value and the base level, given as the average of the left and right minimum values.

T

The medians over 12 years (2000-2011) were used as covariates. The medians were downscaled

RI P

to 250 m resolution using the approach described in Poggio and Gimona (2015) with the land cover map (Morton et al., 2011) as only covariate.

Methods Soil texture transformation

NU

3.1

SC

3

Soil texture was expressed as the relative percentage of sand, silt and clay, with the sum of the

MA

components always equals to the unit (one or 100). They can be treated as compositional variables. They were transformed using Additive Log-Ratio (ALR) transformation. ALR is an extension of the logit transformation when the response variable has more than two components. The

ED

back-transformation is a numerical approximation following GaussHermite quadrature (Aitchison, 1986). ALR was previously applied to soil texture data (Lark and Bishop, 2007; Akpa et al.,

PT

2014; Ballabio et al., 2016) and it was showed (Lark and Bishop, 2007) that ALR transformed variables preserve information on the spatial correlation and maintain the compositional aspect

AC CE

of the variables.

Given a composition of D elements z = [z1 , . . . , zD ], such as zi > 0∀i = 1, . . . , D and PD i=1 zi = 1, ALR is defined as: x = ALR(z) =



ln

z1 zD−1 , . . . , ln zD zD



(2)

The inverse ALR transformation can be described as:

z=

w) exp(w j T exp(w)

(3)

w ) denotes the vector [exp(w w1 ), exp(w w2 ), exp(w wD−1 )] and and j is a vector of length where exp(w D with all elements equal to one (Lark and Bishop, 2007). In case of soil texture, D is equal to 3. In this study, clay was used as the denominator

7

ACCEPTED MANUSCRIPT

variable after testing the possible combination and comparing validation results and model fits

sand clay silt ALR2 = ln clay

T

(not shown here). Therefore the two ALR components that were interpolated can be defined as:

(4)

Soil texture interpolation: full profile and topsoil

SC

3.2

RI P

ALR1 = ln

An extension of the scorpan-kriging approach, i.e. hybrid geostatistical Generalized Additive

NU

Models (GAM Wood, 2006), combining GAM with Gaussian simulations (3DGAM+GS Poggio and Gimona, 2014) was used with , in particular:

MA

1) the fitting of a GAM to estimate the trend of the variable with related covariates; and 2) kriging or Gaussian simulations of GAM residuals as spatial component to account for local details.

ED

The prediction grid had a resolution of 250 m for the lateral dimensions. A smoother of the geographical coordinates was included in the modelling. The prediction matrix of the GAM

PT

model was used to obtain multiple realisation of the trend, according to the approach suggested by Wood (2006). See also Poggio et al. (2010) for an implementation example. Trend estimations were derived by simulating 100 replicate parameter sets from the posterior distribution of the

AC CE

vector of model parameters.

The GAM produced residuals that were further modelled for spatial correlation (e.g. Poggio et al., 2010) using a Gaussian simulation approach (Journel, 1996; Goovaerts, 1997). A variogram (Cressie, 1993) was fitted for the residuals. Exponential and spherical models (Deutsch and Journel, 1998; Goovaerts, 1997) were tested and the model providing the lowest AIC (Akaike Information Criterion; Akaike, 1973) was retained. Anisotropy was also taken into account and the variograms were fitted accounting for the principal anisotropy axes (Goovaerts, 1997). The sum of trend and corresponding residual realisations from the Gaussian simulations was calculated and at each node with 500 iterations in total for each ALR component. The realisations thus obtained were summarised with median and upper (quant95) and lower (quant5) percentiles. The variability of the target property values (ΔV) was estimated as the difference between the

8

ACCEPTED MANUSCRIPT

upper and the lower percentiles. For each node, the deviation of the realisations from the median was also expressed as a percentage as a further measure of the spread of credible model outputs.

T

When considering the full profile, the prediction grid had a vertical resolution of 100cm with

RI P

5cm intervals. The trend was modelled with a 3D spatial approach (Augustin et al., 2009). The model accounts for full 3D-space correlation, exploiting the values of the neighbouring pixels in 3D-space and taking into account the autocorrelated error (Wood, 2006). The trend is modelled

SC

using a scale-invariant tensor product smooth of 3D-space dimension, handling the estimation of the 3D space trend and interaction, as described in Augustin et al. (2009). The 3D-space

NU

smoother used allows separate smoothing parameters and penalties for the 3D-space dimensions and thus avoids the need to make user specific choices about the relative scaling of 3D-space.

MA

The approach does not rely on a regular grid and allows the incorporation of a wide range of correlation structures (Augustin et al., 2009). The GAM implementation used rely on an internal

3.3

ED

cross-validation for model fitting (Poggio and Gimona, 2014).

Mask for organic soils

PT

The mask for organic soils was derived from the work described in Poggio et al. (2013), where the information of a soil being organic was modelled and mapped using a scorpan-kriging ap-

AC CE

proach. The analysis was run for 250m resolution cells. The obtained probability distribution was thresholded using an Area Under the curve (AUC) approach. The area under the ROC (Receiver Operating Characteristic) curve was approximated with a Mann-Whitney U statistic (DeLong et al., 1988). The threshold used in this study was 0.69.

3.4

Soil texture classification

The observed and predicted values of the soil particle size classes (clay, silt and sand) were classified according to the USDA soil texture classification (Soil Survey Division Staff, 1993) (table 4 and fig. 2). The USDA soil texture classification was used instead of the UK Soil Survey of England and Wales texture classification (Natural England, 2011) for easier integration and comparison with international initiatives (Hengl et al., 2014; Viscarra et al., 2015; Ballabio et al., 2016). Soil texture can be plotted on a ternary plot (also called triangle plot). In a ternary plot,

9

ACCEPTED MANUSCRIPT

3D coordinates, whose sum is constant, are projected in the 2D space, using simple trigonometry rules. The texture of a soil sample can be plotted inside a texture triangle. The impact of the

T

modelling uncertainty on the classification was also assessed and maps of probability of a pixel

RI P

to belong to a certain soil texture class were produced. The map of the texture classification obtained from the median of the simulations was also produced.

Validation

NU

3.5

SC

[Figure 2 about here.]

As measure of the spatial variability a spatial structured variance ratio (SSVR) (Vaysse and Lagacherie, 2015) was calculated for the soil properties for the topsoil, the whole profile (3D

MA

variogram) and at certain depths. The SSVR is defined as the different from one of the nugget to sill ratio (Kerry and Oliver, 2008). The values range from zero to one. Values closer to one mean a higher proportion of the data explained by the spatial component.

ED

The results were assessed using in-sample and out-of-sample measures and compared for distribution similarity, spatial structure reproduction and uncertainty ranges. The results of the

PT

modelling were summarised for the whole profile and at five depth layers (i.e. 0-5, 5-15, 15-30, 30-60 and 60-100cm) and compared with values from the validation set at corresponding depths.

AC CE

The main statistics used were:

1) Root Mean Square Error (RMSE), 2) Normalised Root Mean Square Error (nRMSE) using the range of each soil property, 3) R2 derived from a linear model between observed and modelled data (R2 LM ), and 4) The misclassification error and the kappa statistics were calculated (Smeeton, 1985; Congalton, 1991) for the soil texture classes. Kappa coefficient is a statistic which measures inter-rater agreement for qualitative (categorical) items.

3.6

Softwares used

The analyses were performed using open source software: 1) GRASS GIS (GRASS Development Team, 2016) for data management, preparation and visualisation; 10

ACCEPTED MANUSCRIPT

2) the R software (R Core Team, 2016) for the statistical modelling. The following packages were used: i) raster for data management, preparation and visualisation (Hijmans and van

T

Etten, 2013); ii) mgcv for GAM (Wood, 2006); iii) geoR (Ribeiro and Diggle, 2001) for fitting

RI P

the variograms of the residuals; iv) gstat (Pebesma, 2004) for kriging; v) rgdal for data management (Keitt et al., 2009); vi) soiltexture for the soil texture classification (Moeys,

AC CE

PT

ED

MA

NU

SC

2015).

11

ACCEPTED MANUSCRIPT

Results and discussion

4

Properties distribution and predictors

T

4.1

RI P

The distribution at the sampled locations of the three particle size classes considered are presented in fig. 3. Sand showed a higher median and wider range, while clay and silt generally showed lower values. The spatial structure of the data was explored with the auto- and cross-variograms

SC

(fig. 4). Silt and clay showed covariance decreasing with distance, while sand showed covariance increasing with distance (repulsion) with the other particle sizes. The SSVR values (table 2) show

NU

that a large proportion of the soil particle classes were spatially structured. The relative nugget effect i.e., the complement to 1 of the spatially structured variance, relates to measurement errors

MA

and spatial sources of variation at distances smaller than the shortest sampling interval. Such large values for SSVR could indicate that despite the large number of sampled locations, there are not sufficient information to describe the short range variability of the values.

ED

[Table 1 about here.]

PT

[Figure 3 about here.] [Figure 4 about here.]

AC CE

Table 1 showed the significant covariates for the topsoil and full profile models for both the ALR components eq. 4. The explained deviance was similar for both components and higher for the topsoil-only model. This behaviour was already found in previous studies (Vaysse and Lagacherie, 2015; Adhikari et al., 2013) and can be explained with the higher number of profiles with available data in the topsoil layer. The 3D model had a lower explained deviance despite having more significant variables. However the AIC of the models increased when dropping variables. The lower explained deviance can be explained with the weaker relationship between deeper layers and the selected covariates. The significant predictors were similar for both components of the ALR transformation, with morphological features and vegetation indices describing most of the variability. The water index was only selected for the 3D models. The land surface temperature was always significant acting as a proxy for climate conditions. The phenological information was used in the models with 12

ACCEPTED MANUSCRIPT

the length of the season always significant, while amplitude was never selected. The selected significant covariates described the different soil forming factors and in particular the differences

T

in vegetation and climate information. Some of the RS covariates are mainly describing the

RI P

vegetation cover of the soils. This is particularly useful for temperate regions where the soil is vegetated for most of the year with different ratios of evergreen and broadleaved and therefore different seasonality. Furthermore in this study the variables were choose to maximise the

SC

predictive power of the models not their explanatory capabilities. Covariates were only derived from morphological features and optical remote sensing. Further improvements in the models

NU

could be reached using radar remote sensing (e.g. Niang et al., 2014). Radar remote sensing is an attractive alternative as an environmental covariate for soil surface texture mapping be-

MA

cause it is not significantly affected by atmospheric attenuation and clouds, and at a wavelength of 5 cm (C-band), the data can provide estimates of near-surface soil moisture to a depth of

4.2

Topsoil model

ED

approximately 0 to 5 cm (Moran et al., 2004).

PT

The validation assessment of the back-transformation of the ALR components showed good results summarised in table 2 and in fig. 5. Sand is the soil particle class with the least accurate

AC CE

validation results. The RMSE is the highest and almost double the values obtained for silt and clay. The original range of sand is also wider compared to silt and clay. This is confirmed by the results of the normalised RMSE, where the value for sand is only slightly higher. The qq-plot (fig. 5) showed the correspondence between the measured and the predicted values for all three particle classes for most of the range. There is over-prediction of lower values for sand and a corresponding under-prediction for the silt. The distribution of clay values was closer to the one to one line for most of its range. The largest discrepancies can be found in the mountain areas where there is alternance of less developed soils with organic and mode developed soils. The validation results with the high values for SSVR suggest that the covariates used do not capture all the small range variability. The use of further covariates may help to improve the results. The reproduction of the spatial structure was also assessed in the validation, comparing the variograms of predicted and validation sets (fig. 5). As expected, the variograms have a lower

13

ACCEPTED MANUSCRIPT

sill because the predicted values are smoother, but they follow a similar shape of the variograms from the validation set.

T

The final maps are presented in fig. 6 with a mask for regions with organic soils. The maps

RI P

show the higher content of sand in large part of Scottish soils. Higher percentages of silt and clay are modelled in the regions with higher productivity soils, i.e. in the north-east and the Eastern coast in the flatter areas. Sand shows also the higher variability, i.e. the differences

SC

between 95th and 5th quantiles. This is related with the larger range of sand values. Higher uncertainty for clay values are present in the mountain soils where the clay values are normally

NU

reduced. Fig. 6 shows the delta (difference between 95 t h and 5t h percentiles) as percent of the median. Sand and silt had a percentage variability of about 5% of the median values, while

MA

clay presented much higher values with a maximum value of 70% and a median of 20%. The area with higher percentage variability for sand and silt are close to the north coast and in the areas where agricultural soils are. The areas with higher percentage variability for clay are in

ED

the mountains, where little amount of clay can be found. This study used GAMs instead of a regression tree approach (e.g. Adhikari et al., 2013; Akpa et al., 2014; Chagas et al., 2016; Niang

AC CE

(Poggio et al., 2013).

PT

et al., 2014). In previous studies, the use of GAM provided more accurate results in Scotland

[Figure 5 about here.] [Table 2 about here.] [Figure 6 about here.]

14

ACCEPTED MANUSCRIPT

4.3

Whole profile model (3D).

Table 3 shows the validation summary statistics for the 3D model at five depths. The RMSE

T

results are higher for sand and silt both for absolute and normalised values. The validation is less

RI P

accurate with depth showing increased values of RMSE and decreasing values of R2 . SSVR values show a large proportion of the variance explained by the spatial variation, especially for clay.

SC

Values of SSVR increased in deeper layers. This can be explained by the lower number of points available and the weaker relationships with the selected environmental covariates (Adhikari et al., 2013; Vaysse and Lagacherie, 2015).

NU

Fig. 7 shows the maps for the whole profile model at the five considered depths. The spatial pattern is similar to the one described for the topsoil model (fig. 6), with decreasing contents of

MA

silt and clay at lower depths.

[Table 3 about here.]

AC CE

PT

ED

[Figure 7 about here.]

15

ACCEPTED MANUSCRIPT

4.4

Soil texture classes

Fig. 8 shows the texture triangle for both validation and predicted values. The patterns look

T

similar with most results in sand and sandy loam classes. The predicted values have more values

RI P

in loam, silty loam and sandy loam classes. These results showed patterns similar to the qq-plots in fig. 5. The over-prediction of lower values of sand and corresponding under-prediction for high

SC

values of silt can explain the higher number of points in the more loamy texture classes. Fig. 10 shows the distribution of pixels in the soil texture classes across the 500 iterations. For most of the classes, the corresponding values calculated from the validation dataset were included in the

NU

whiskers range of the boxplot. The main exceptions were Sandy loam and, in minor measure, Sandy clay loam and Loamy sand classes. This is the consequence of the overestimation of sand

MA

for the lower values of the distribution.

Fig. 9 shows the map of the soil texture classes obtained from the median of the simulations. The spatial patterns are comparable with the map presented in The James Hutton Institute

ED

(2014). The map produced in this study is less smooth, and the pattern can be explained by the cell resolution and the approach used reproducing the spatial variability. Fig. 11 shows the

PT

probability of a pixel to belong to a soil texture class. The probability was derived over the 500 iterations. The spatial patterns are similar to the ones shows in fig. 6 for the separate soil particle

soils.

AC CE

components. The hotspots of finer texture soils are located in the regions of higher productivity

[Figure 8 about here.] [Figure 9 about here.] [Figure 10 about here.] [Figure 11 about here.]

16

ACCEPTED MANUSCRIPT

5

Conclusions

T

This study used available environmental information to map 3D soil texture with uncertainty

RI P

assessment and evaluation of the uncertainty propagation on soil texture classification. It showed an operational application of the previous developed GAM+GS approach to compositional data in Scotland, a region with high soil spatial variability. The results obtained showed good agree-

SC

ment between values at validation locations with the R2 calculated with the validation dataset between 0.55 and 0.60 and the RMSE values below 13%.

NU

The accuracy of DSM results depends on the statistical relationships between measured soil observations and environmental covariates. It is difficult to identify a set of covariates that fully describe the relationships with soil-forming factors. The set of covariates used in this study

MA

explained about 40% of the variance of the data. Further covariates derived from radar remote sensing at different resolution (e.g. Sentinel1, RadarSAT2, SMOS, SMAP) could be helpful to further describe the landscape patterns of soil texture.

ED

The data modelled are complementary with continental and global products providing more detailed and regionally-specific results. The modelled particle size distribution datasets can be

PT

used as input for further modelling in a number of areas, such as soil and water management, hydrology, agricultural and crop modelling ecosystem services and in assessment of soil erosion

AC CE

risk. The accompanying uncertainty can be used to provide supporting information for land management choices.

17

ACCEPTED MANUSCRIPT

Acknowledgements

T

This work was funded by the Scottish Government’s Rural and Environment Science and Analyt-

RI P

ical Services division. Many thanks to the team that sampled and analysed the soils and to the team that set up the database. MODIS data are distributed by the Land Processes Distributed Active Archive Centre (LP DAAC), located at the U.S. Geological Survey (USGS) Earth Re-

SC

sources Observation and Science (EROS) Center (lpdaac.usgs.gov). Thanks are due to Dr Mark

NU

Brewer (BioSS) for comments on an early version of this manuscript.

References

MA

Adhikari, K., Hartemink, A. E., 2016. Linking soils to ecosystem services - A global review. Geoderma 262, 101–111.

Adhikari, K., Kheir, R. B., Greve, M. B., Bocher, P. K., Malone, B. P., Minasny, B., McBratney,

ED

A. B., Greve, M. H., 2013. High-Resolution 3-D Mapping of Soil Texture in Denmark. Soil Science Society of America Journal 77 (3), 860–876.

PT

Aitchison, J., 1986. The statistical analysis of compositional data. Chapman & Hall, London.

AC CE

Akaike, H., 1973. Information theory and an extension of the maximum likelihood principle. In: Petrov, B., Csaki, F. . (Eds.), Second International Symposium on Information Theory. Akademia Kiado, Budapest, Hungary, pp. 267–281. Akpa, S. I. C., Odeh, I. O. A., Bishop, T. F. A., Hartemink, A. E., 2014. Digital Mapping of Soil Particle-Size Fractions for Nigeria. Soil Science 78 (6), 1953–1966. Arrouays, D., Grundy, M. G., Hartemink, A. E., Hempel, J. W., Heuvelink, G. B., Hong, S. Y., Lagacherie, P., Lelyk, G., McBratney, A. B., McKenzie, N. J., d.L. Mendonca-Santos, M., Minasny, B., Montanarella, L., Odeh, I. O., Sanchez, P. A., Thompson, J. A., Zhang, G.-L., 2014. Chapter three - globalsoilmap: Toward a fine-resolution global grid of soil properties. Vol. 125 of Advances in Agronomy. Academic Press, pp. 93 – 134. Augustin, N., Musio, M., von Wilpert, K., Kublin, E., Wood, S., Schumacher, M., 2009. Modeling

18

ACCEPTED MANUSCRIPT

spatiotemporal forest health monitoring data. Journal of the American Statistical Association 104 (487), 899–911.

RI P

ropean scale using the LUCAS database. Geoderma 261, 110–123.

T

Ballabio, C., Panagos, P., Monatanarella, L., 2016. Mapping topsoil physical properties at Eu-

Bationo, A., Kihara, J., Vanlauwe, B., Waswa, B., Kimetu, J., 2007. Soil organic carbon dy-

SC

namics, functions and management in West African agro-ecosystems. Agricultural Systems 94 (1, SI), 13–25, Scientific Workshop on Land Management for Carbon Sequestration in the

NU

Sudan-Sahel Region, Bamako, Mali, 2004.

Ben-Dor, E., 2002. Quantitative remote sensing of soil properties. In: Sparks, DL (Ed.), Advances

MA

in Agronomy. Vol. 75 of Advances in Agronomy. pp. 173–243. Bockheim, J. G., Hartemink, A. E., 2013. Distribution and classification of soils with clay-

ED

enriched horizons in the USA. Geoderma 209, 153–160. Buchanan, S., Triantafilis, J., Odeh, I. O. A., Subansinghe, R., 2012. Digital soil mapping of compositional particle-size fractions using proximal and remotely sensed ancillary data. Geophysics

PT

77 (4, S), WB201–WB211.

AC CE

Calzolari, C., Ungaro, F., Filippi, N., Guermandi, M., Malucelli, F., Marchi, N., Staffilani, F., Tarocco, P., 2016. A methodological framework to assess the multiple contributions of soils to ecosystem services delivery at regional scale. Geoderma 261, 190–203. Chagas, C. d. S., de Carvalho Junior, W., Bhering, S. B., Calderano Filho, B., 2016. Spatial prediction of soil surface texture in a semiarid region using random forest and multiple linear regressions. Catena 139, 232–240. Congalton, R., 1991. A review of assessing the accuracy of classifications of remotely sensed data. Remote Sensing of Environment, 37, 35–46. Cressie, N., 1993. Statistics for Spatial Data. Wiley, New York. Dane, J., Topp, G., Campbell, G., Horton, R., Jury, W., Nielsen, D., van Es, H., Wierenga, P., Topp, G., 2002. Part 4, physical methods. Methods Soil Analytics. 19

ACCEPTED MANUSCRIPT

DeLong, E., Delong, D., Clarke-Pearson, D., 1988. Comparing areas under two or more correlated Receiver Operating Characteristic curves: a nonparametric approach. Biometrics 44 (3), 837–

T

845.

edition Edition. Oxford University Press, New York.

RI P

Deutsch, C., Journel, A., 1998. GSLIB: Geostatistical Software Library and User’s Guide, second

SC

Eklundh, L., Jonsson, P., 2011. Timesat 3.1 Software Manual. Tech. rep., Lund University, Sweden.

NU

Gao, B., 1996. NDWI - a normalized difference water index for remote sensing of vegetation liquid water from space. Remote Sensing of Environment 58, 257 – 266.

MA

Gooley, L., Huang, J., Pag, D., Triantafilis, J., 2014. Digital soil mapping of available water content using proximal and remotely sensed data. Soil Use and Management 30 (1), 139–151.

ED

URL http://dx.doi.org/10.1111/sum.12094

Goovaerts, P., 1997. Geostatistics for Natural Resources Evaluation. Oxford University Press.

PT

GRASS Development Team, 2016. Geographic Resources Analysis Support System (GRASS GIS) Software, version 7.0.3.

AC CE

URL http://www.grass.osgeo.org Gu, Y., Hunt, E., Wardlow, B., Basara, J. B., Brown, J. F., Verdin, J. P., 2008. Evaluation of MODIS NDVI and NDWI for vegetation drought monitoring using Oklahoma Mesonet soil moisture data. Geophysical Research Letters 35 (22), L22401. Hartemink, A. E., McBratney, A., 2008. A soil science renaissance. Geoderma 148 (2), 123–129. Hassink, J., APR 1997. The capacity of soils to preserve organic C and N by their association with clay and silt particles. Journal of Plant and Soil 191 (1), 77–87. Hengl, T., de Jesus, J. M., MacMillan, R. A., Batjes, N. H., Heuvelink, G. B. M., Ribeiro, E., Samuel-Rosa, A., Kempen, B., Leenaars, J. G. B., Walsh, M. G., Gonzalez, M. R., AUG 29 2014. SoilGrids1km-Global Soil Information Based on Automated Mapping. PLOS ONE 9 (8).

20

ACCEPTED MANUSCRIPT

Heuvelink, G., Huisman, J., 2000. Quantifying spatial uncertainty in natural resources: theory and applications for GIS and remote sensing. Ann Arbor Press, Ch. Choosing between abrupt

T

and gradual spatial variation?, pp. 111–117.

version 2.1-25.

SC

URL http://CRAN.R-project.org/package=raster

RI P

Hijmans, R. J., van Etten, J., 2013. raster: Geographic data analysis and modeling. R package

Huete, A., Didan, K., Miura, T., Rodriguez, E. P., Gao, X., Ferreira, L. G., 2002. Overview of the

NU

radiometric and biophysical performance of the MODIS vegetation indices. Remote Sensing of Environment 83 (1-2), 195 – 213.

MA

Jarvis, A., Reuter, H., Nelson, A., Guevara, E., 2006. Hole-filled seamless SRTM data V3. Tech. rep., International Centre for Tropical Agriculture (CIAT). J¨onsson, P., Eklundh, L., 2002. Seasonality extraction by function fitting to time-series of satellite

ED

sensor data. IEEE Transactions on Geoscience and Remote Sensing 40, 1824–1832. J¨onsson, P., Eklundh, L., 2004. TIMESAT - a program for analyzing time-series of satellite sensor

PT

data. Computer and Geosciences 30, 833–845.

AC CE

Journel, A., 1996. Modelling uncertainty and spatial dependence: Stochastic imaging. International Journal of Geographical Information Systems 10 (5), 517–522. Keitt, T., Bivand, R., Pebesma, E., Rowlingson, B., 2009. rgdal: Bindings for the Geospatial Data Abstraction Library. R package version 0.6-21. URL http://CRAN.R-project.org/package=rgdal Kerry, R., Oliver, M., 2008. Determining nugget:sill ratios of standardized variograms from aerial photographs to krige sparse soil data. Precision Agriculture 9, 33–56. Lagacherie, P., Baret, F., Feret, J.-B., Netto, J. M., Robbez-Masson, J. M., 2008. Estimation of soil clay and calcium carbonate using laboratory, field and airborne hyperspectral measurements. Remote 112 (3), 825–835.

21

ACCEPTED MANUSCRIPT

Lark, R., Bishop, T., 2007. Cokriging particle size fractions of the soil. European Journal of Soil Science 58, 763–774.

T

Liess, M., Glaser, B., Huwe, B., 2012. Uncertainty in the spatial prediction of soil texture.

RI P

Comparison of regression tree and Random Forest models. Geoderma 170, 70–79. Lilly, A., Bell, J., Hudson, G., Nolan, A., Towers, W., 2010. National Soil Inventory of Scotland

SC

1 (NSIS1): site location, sampling and profile description. (1978-1998). Tech. rep., Macaulay Institute.

NU

Malone, B., McBratney, A., Minasny, B., Laslett, G., 2009. Mapping continuous depth functions of soil carbon storage and available water capacity. Geoderma 154 (1-2), 138–152.

MA

McBratney, A., Mendonca-Santos, L., Minasny, B., 2003. On digital soil mapping. Geoderma 117, 3–52.

ED

McLauchlan, K. K., DEC 2006. Effects of soil texture on soil carbon and nitrogen dynamics after cessation of agriculture. Geoderma 136 (1-2), 289–299.

PT

Minasny, B., Hartemink, A. E., 2011. Predicting soil properties in the tropics. Earth-Science Reviews 106 (1-2), 52–62.

AC CE

Moeys, J., 2015. soiltexture: Functions for Soil Texture Plot, Classification and Transformation. R package version 1.3.3.

URL https://CRAN.R-project.org/package=soiltexture Moran, M. S., Peters-Lidard, C. D., Watts, J. M., McElroy, S., 2004. Estimating soil moisture at the watershed scale with satellite-based radar and land surface models. Canadian Journal of Remote Sensing 30 (5), 805–826. Morton, D., Rowland, C., Wood, C., Meek, L., Marston, C., Smith, G., Wadsworth, R., Simpson, I., 2011. Final report for LCm2007 - the new uk land cover map. Tech. Rep. Countryside Survey Technical Report No 11/07 112pp. (CEH Project Number: C03259)., NERC/Centre for Ecology & Hydrology.

22

ACCEPTED MANUSCRIPT

Mulder, V., de Bruin, S., Schaepman, M., Mayr, T., 2011. The use of remote sensing in soil and terrain mapping - A review. Geoderma 62, 1–19.

RI P

http://publications.naturalengland.org.uk/file/83081.

T

Natural England, 2011. Natural England Technical Information Note TIN037 - Soil texture.

Niang, M. A., Nolin, M. C., Jego, G., Perron, I., 2014. Digital Mapping of Soil Texture Using

SC

RADARSAT-2 Polarimetric Synthetic Aperture Radar Data. Soil Science Society of America Journal 78 (2), 673–684.

NU

O’Callaghan, J. F., Mark, D. M., 1984. The extraction of drainage networks from digital elevation data. Computer Vision Graphics Image Processing 28, 323–344.

MA

Orton, T. G., Pringle, M. J., Bishop, T. F. A., 2016. A one-step approach for modelling and mapping soil properties based on profile data sampled over varying depth intervals. Geoderma

ED

262, 174–186.

Ouerghemmi, W., Gomez, C., Naceur, S., Lagacherie, P., 2016. Semi-blind source separation for the estimation of the clay content over semi-vegetated areas using vnir/swir hyperspectral

PT

airborne data. Remote Sensing of Environment 181, 251 – 263.

AC CE

Panagos, P., Meusburger, K., Ballabio, C., Borrelli, P., Alewell, C., 2014. Soil erodibility in Europe: A high-resolution dataset based on LUCAS. Science of the Total Environment 479, 189–200.

Pebesma, E., 2004. Multivariable geostatistics in S: the gstat package. Computers & Geosciences 30, 683–691.

Poggio, L., Gimona, A., 2014. National scale 3D modelling of soil organic carbon stocks with uncertainty propagation - An example from Scotland. Geoderma 232-234, 284–299. Poggio, L., Gimona, A., 2015. Downscaling and correction of regional climate models outputs with a hybrid geostatistical approach. Spatial Statistics 14 (A), 4–21. Poggio, L., Gimona, A., Brewer, M., 2013. Regional scale mapping of soil properties and their uncertainty with a large number of satellite-derived covariates. Geoderma 209-210, 1–14. 23

ACCEPTED MANUSCRIPT

Poggio, L., Gimona, A., Brown, I., 2012. Spatio-temporal MODIS EVI gap filling under cloud cover: an example in Scotland. ISPRS Journal of Photogrammetry and Remote Sensing 72,

T

56–72.

RI P

Poggio, L., Gimona, A., Brown, I., Castellazzi, M., 2010. Soil available water capacity interpolation and spatial uncertainty modelling at multiple geographical extents. Geoderma 160,

SC

175–188.

R Core Team, 2016. R: A Language and Environment for Statistical Computing. R Foundation

NU

for Statistical Computing, Vienna, Austria, ISBN 3-900051-07-0. URL http://www.R-project.org/

MA

Ribeiro, P., Diggle, P., 2001. geoR: a package for geostatistical analysis. R-NEWS 1 (2), 14–18. URL http://CRAN.R-project.org/doc/Rnews/

Rodriguez, E., Morris, C., Belz, J., Chapin, E., Martin, J., Daffer, W., Hensley, S., 2006. An

ED

assessment of the SRTM topographic products. Tech. Rep. JPL D-31639, NASA-Jet Propulsion Laboratory.

PT

Running, S., R. R., Nemani, F. A., Heinsch, M., Zhao, M., Reeves, Hashimoto., H., 2004. A continuous satellite-derived measure of global terrestrial primary production. BioScience

AC CE

54 (6), 47–560.

Running, S., Nemani, R., Glassy, J., Thornton, P., 1999. MODIS daily photosynthesis (PSN) and annual net primary production (NPP) product Algorithm Theoretical Basis Document. Tech. rep.

Selige, T., Boehner, J., Schmidhalter, U., 2006. High resolution topsoil mapping using hyperspectral image and field data in multivariate regression modeling procedures. Geoderma 136 (1-2), 235–244. Smeeton, N., 1985. Early History of the Kappa Statistic. Biometrics 41, 795. Soil Survey Division Staff, 1993. Soil survey manual, volume handbook 18,chapter 3. Tech. rep., Soil Conservation Service. U.S. Department of Agriculture.

24

ACCEPTED MANUSCRIPT

The James Hutton Institute, 2014. Map of topsoil texture in Scotland using British Standards Institute classes. http://www.soils-scotland.gov.uk/.

T

Thompson, J., Roecker, S., Grunwald, S., Owens, P., 2012. Hydropedology. Elsevier, Ch. Digital

RI P

Soil Mapping: Interactions with and Applications for Hydropedology.

Toth, B., Weynants, M., Nemes, A., Mako, A., Bilas, G., Toth, G., 2015. New generation of

SC

hydraulic pedotransfer functions for Europe. European Journal of Soil Science 66 (1, SI), 226–238.

NU

van Meirvenne, M., van Cleemput, I., 2005. Environmental Soil-Landscape Modeling. CRC Press, Ch. Pedometrical Techniques for Soil Texture Mapping at Different Scales, pp. 323–342.

MA

Vaysse, K., Lagacherie, P., 2015. Evaluating digital soil mapping approaches for mapping globalsoilmap soil properties from legacy data in Languedoc-Roussillon (France). Geoderma regional

ED

4, 20–30.

Viscarra, R. A., Chen, C., Grundy, M. J., Searle, R., Clifford, D., Campbell, P. H., 2015. The Australian three-dimensional soil grid: Australia’s contribution to the GlobalSoilMap project.

PT

SOIL RESEARCH 53 (8), 845–864.

AC CE

Wan, Z., 1999. MODIS Land-Surface Temperature Algorithm Theoretical Basis Document (LST ATBD). Tech. rep., NASA. Wood, S., 2006. Generalized Additive Models: An Introduction with R. Chapman and Hall/CRC Press.

Wosten, J., Pachepsky, Y., Rawls, W., 2001. Pedotransfer functions: bridging the gap between available basic soil data and missing soil hydraulic characteristics. Journal of Hydrology 251 (34), 123–150.

25

ACCEPTED MANUSCRIPT

List of Tables

T

RI P

SC NU MA ED

4

PT

2 3

Significance of covariates for models of the ALR components for both the topsoil and the whole profile models. . . . . . . . . . . . . . . . . . . . . . . . . . . . . Validation summary statistics and SVVR values for the topsoil model (2D). . . Validation summary statistics and SVVR values for the whole profile model (3D) at 5 depths. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Texture classes for the system used. . . . . . . . . . . . . . . . . . . . . . . . . .

AC CE

1

27

28 29 30 33

NU

SC

RI P

T

ACCEPTED MANUSCRIPT

MA

2D

Coordinates smoother Elevation Slope Topographic wetness index

PT

ED

EVI NDWI Primary productivity Amplitude EVI LOS EVI LST d

Explained deviance (%)

3D

ALR1

ALR2

ALR1

ALR2

Y (x,y) -

Y (x,y) Y Y

Y(x,y,z) Y Y Y

Y(x,y,z) Y Y Y

Y Y Y Y

Y Y Y

Y Y Y Y Y

Y Y Y Y Y

40.4

42.2

28.5

26.7

AC CE

Table 1: Significance of covariates for models of the ALR components for both the topsoil and the whole profile models.

28

RMSE

nRMSE

R2LM

SSVR

ED

MA

NU

SC

RI P

T

ACCEPTED MANUSCRIPT

0.13 0.10 0.11

0.58 0.56 0.57

0.81 0.74 0.77

Sand Silt Clay

12.21 7.34 6.95

AC CE

PT

Table 2: Validation summary statistics and SVVR values for the topsoil model (2D).

29

RMSE

nRMSE

Sand

12.75 15.51 14.54 16.08 20.44

0.17 0.16 0.16 0.17 0.22

Silt

9.96 11.42 10.36 11.53 13.91

0.61 0.62 0.72 0.62 0.89

0.47 0.40 0.36 0.31 0.26

0-5 5-15 15-30 30-60 60-100

0.27 0.15 0.14 0.15 0.18

0.71 0.72 0.81 0.88 0.90

0.40 0.39 0.33 0.30 0.28

0-5 5-15 15-30 30-60 60-100

0.13 0.14 0.14 0.15 0.16

0.89 0.84 0.89 0.87 0.89

0.45 0.43 0.35 0.32 0.27

0-5 5-15 15-30 30-60 60-100

MA

R2LM

ED

PT

6.24 7.51 7.05 7.15 7.57

Depth (cm)

SVRR

AC CE

Clay

NU

SC

RI P

T

ACCEPTED MANUSCRIPT

Table 3: Validation summary statistics and SVVR values for the whole profile model (3D) at 5 depths.

30

ACCEPTED MANUSCRIPT

List of Figures

8 9 10 11

T

RI P

SC

7

NU

6

MA

5

ED

4

PT

2 3

Maps of test area: sampling locations with texture information available (black circles) and sampling locations with organic soils without texture information (blue triangles). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Texture triangle (TT) with limits for the classes considered. . . . . . . . . . . . . Values distributions for the three considered particle sizes (sand left, silt centre, clay right). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Auto- and Cross- empirical variograms for the three particle size classes (observed values). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Validation results for the topsoil model. QQ-plots between predicted and validation data (left) and comparison of variograms of the predictions (right). . . . . . Maps of the topsoil particle size classes in mainland Scotland (sand left, silt centre, clay right) with mask for organic soils (grey pixels). . . . . . . . . . . . . . . . . Maps of the particle size classes in mainland Scotland at five depths (sand left, silt centre, clay right) with masks for organic soils and depth (only the pixels with depth >= the threshold are shown). The five depths are 0-5,5-15,15-30,30-60,60100cm. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Texture triangle for validation (left) and prediction (right) values. The misclassification error is 12 % and the kappa is 0.64. . . . . . . . . . . . . . . . . . . . . . Map of soil tetxure classes obtained from the median of the simulations. . . . . . Distribution of pixels in the soil texture classes for the 500 iterations. The red square indicates the values of pixels in the validation dataset. . . . . . . . . . . . Maps of probabilities of belonging to a soil texture class, derived from the 500 iterations (topsoil). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

AC CE

1

31

32 33 34 35 36 37

38 39 40 41 42

AC CE

PT

ED

MA

NU

SC

RI P

T

ACCEPTED MANUSCRIPT

Figure 1: Maps of test area: sampling locations with texture information available (black circles) and sampling locations with organic soils without texture information (blue triangles).

32

Full name

Cl SiCl SaCl ClLo SiClLo SaClLo Lo SiLo SaLo Si LoSa Sa

clay silty clay sandy clay clay loam silty clay loam sandy clay loam loam silty loam sandy loam silt loamy sand sand

AC CE

PT

ED

Abbreviation

MA

NU

SC

RI P

T

ACCEPTED MANUSCRIPT

Table 4: Texture classes for the system used.

Figure 2: Texture triangle (TT) with limits for the classes considered.

33

AC CE

PT

ED

MA

NU

SC

RI P

T

ACCEPTED MANUSCRIPT

Figure 3: Values distributions for the three considered particle sizes (sand left, silt centre, clay right).

34

AC CE

PT

ED

MA

NU

SC

RI P

T

ACCEPTED MANUSCRIPT

Figure 4: Auto- and Cross- empirical variograms for the three particle size classes (observed values).

35

NU

SC

RI P

T

ACCEPTED MANUSCRIPT

(b) Silt

AC CE

PT

ED

MA

(a) Sand

(c) Clay

Figure 5: Validation results for the topsoil model. QQ-plots between predicted and validation data (left) and comparison of variograms of the predictions (right).

36

NU

SC

RI P

T

ACCEPTED MANUSCRIPT

AC CE

PT

ED

MA

(a) Median

(b) Delta (difference between 95 t h and 5t h percentiles)

(c) Delta as percent of the median

Figure 6: Maps of the topsoil particle size classes in mainland Scotland (sand left, silt centre, clay right) with mask for organic soils (grey pixels). 37

AC CE

PT

ED

MA

NU

SC

RI P

T

ACCEPTED MANUSCRIPT

Figure 7: Maps of the particle size classes in mainland Scotland at five depths (sand left, silt centre, clay right) with masks for organic soils and depth (only the pixels with depth >= the threshold are shown). The five depths are 0-5,5-15,15-30,30-60,60-100cm.

38

AC CE

PT

ED

MA

NU

SC

RI P

T

ACCEPTED MANUSCRIPT

Figure 8: Texture triangle for validation (left) and prediction (right) values. The misclassification error is 12 % and the kappa is 0.64.

39

AC CE

PT

ED

MA

NU

SC

RI P

T

ACCEPTED MANUSCRIPT

Figure 9: Map of soil tetxure classes obtained from the median of the simulations.

40

AC CE

PT

ED

MA

NU

SC

RI P

T

ACCEPTED MANUSCRIPT

Figure 10: Distribution of pixels in the soil texture classes for the 500 iterations. The red square indicates the values of pixels in the validation dataset.

41

AC CE

PT

ED

MA

NU

SC

RI P

T

ACCEPTED MANUSCRIPT

Figure 11: Maps of probabilities of belonging to a soil texture class, derived from the 500 iterations (topsoil).

42

ACCEPTED MANUSCRIPT

Highlights

T

• 3D modelling and mapping soil particle classes for Scotland at 250m resolution

RI P

• Evaluation of the impact of uncertainty on soil texture classification

AC CE

PT

ED

MA

NU

SC

• Remote sensing and morphology to explain about 40% of the variance of the data

26