Estimation and mapping of surface soil properties in the Caucasus Mountains, Azerbaijan using high-resolution remote sensing data

Highlights

• Terrain attributes gain more spatial information than spectral indices.

• High resolution remote sensing data are promising to map mountain soils.

• Assumptions made for uncertainty modelling is significantly important.

Abstract

Soil surveys and mapping with traditional methods are time-consuming and expensive especially in mountainous areas while demand for detailed soil information is steadily increasing. This study tested two spatial hybrid approaches to predict and map basic soil properties using high resolution digital elevation model (DEM) and multispectral satellite imagery in a study area located in the Caucasus Mountains, Azerbaijan. Terrain attributes and spectral indices extracted from DEM with 12.5 m spatial resolution and Pléiades-1 data were used as auxiliary variables. A total of 115 soil samples were collected from the surface layer of 423 ha area and tested for soil organic carbon, soil reaction (pH in H₂O and KCl solutions), calcium carbonate (CaCO₃), sand, silt, clay and hygroscopic water content. The predictive capability of Universal Kriging (UK) and Random Forest Kriging (RFK) was evaluated using spatial cross-validation technique. To model and quantify the associated uncertainty of these models a probabilistic framework, kriging variance approach was applied. The uncertainty models were validated using independent and randomly selected control points (20% of the reference samples). For this, the actual fraction of true values falling within symmetric prediction intervals was calculated and visualized known as accuracy plot. Although the performances of the tested models were similar, RFK was superior in view of both accuracy and computed biases. The models were capable of delineating spatial pattern, mostly elevation dependent as well as the local patterns attributed by e.g., variations in vegetation, land use and soil erosion. UK model produced a few local erratic spatial patterns (e.g., in the case of pH) corresponding to the artifacts such as roads and houses in the image that should be considered in future applications. When comparing the uncertainties, both the models produced considerable underestimations and overestimations depending on soil property. RFK provided better uncertainty estimation for the most of soil properties than UK, the latter technique was more appropriate for the clay and pH_KCl prediction. This case study confirmed the importance of assumptions made in uncertainty modelling and quantification. Those soil properties were therefore reliably predicted that their residuals were compatible with the normality assumption and showed apparent spatial correlation, e.g., both the models severely overestimated uncertainty of CaCO₃ due to lack of normality assumption and low spatial correlation. This study showed that high resolution remote sensing data are promising, and the procedure presented in this study can be reliably used to map the studied soil properties and extended to partially larger adjacent areas characterized by similar environmental conditions in the Caucasus Mountains. However, with respect to future digital soil mapping, we assume that it is important to consider sampling design, testing other modelling approaches their uncertainties and multi-scale digital terrain analysis as well.

Introduction

The soil, particularly in the mountainous regions is an important controlling factor of many environmental processes such as forest and crop growth, hydro-geochemical cycle, nutrient leaching and greenhouse gas emission (Lal, 2004; Wilding, 1985). The quantitative estimation of spatial variability of soil is important in a better understanding of complex relationships between soil properties and environmental factors. This complexity includes intrinsic and extrinsic factors (Heuvelink and Webster, 2001) varying greatly depending on topography, climate, vegetation and anthropogenic activity all of which significantly affect spatial variability of soil properties (Shi et al., 2009).

Availability of numerous digital covariates as well new types of space-born remote sensors and advances in computers and information technologies have created a quantitative approach in soil mapping, termed digital soil mapping (DSM) (McBratney et al., 2003). By making use of large diversity of auxiliary data sources and modelling approaches from simple linear models to complex machine learning techniques, DSM has shifted from a research phase into operational use (Minasny and McBratney, 2016) and become comparable, more accurate and an indispensable approach compared to conventional mapping (Bazaglia Filho et al., 2013; Collard et al., 2014).

Overall, DSM especially in mountainous areas, is a complex task due to heterogeneity, the nature of soil-forming factors, and the close relationships that occurred between them. This is especially complex in the Caucasus Mountains where the local communities host traditional agriculture, and typical land ownership is characterized by small patches (Salukvadze and Medvedkov, 2011) when high resolution soil mapping is considered. On the other hand, soil surveys in mountainous areas can be costly, and DSM usually faces a general lack of available observations that often limits the quantity and quality of the samples on which the mapping procedure is based. In this sense, high resolution DSM in mountainous areas can contribute to the proper evaluation and sustainable management of soil resources as well different environmental modelling applications.

In mountainous areas, the topography plays a crucial role and influences significantly most of the other soil-forming factors such as e.g., climate and vegetation cover (Pendleton and Jenny, 1945). The character of such close relations between soil and environmental factors encourages considerably linear regression technique (McKenzie and Ryan, 1999), whereas non-linear relationships between soil properties and environmental factors usually occur and reduce the applicability of the linear regression to limited areas (Lian et al., 2009). These problems can be partially solved by non-parametric machine learning techniques and hybrid spatial models as well auxiliary data used to compensate for the data shortage.

A variety of approaches, such as traditional statistical (Lian et al., 2009), geostatistical and machine learning techniques have been extensively applied in DSM (Keskin and Grunwald, 2018; Veronesi and Schillaci, 2019; Wiesmeier et al., 2011). The geostatistical approach is one of the best spatial interpolation techniques (Webster and Oliver, 2008) and its reliable performance owing to e.g., sampling design, the nature of the study area was presented in many studies in particular with main focus on soil organic carbon mapping (Cambule et al., 2014; Pouladi et al., 2019; Veronesi and Schillaci, 2019). However, it was stressed that kriging becomes uncertain with a small number of samples and unable to model local variations (stationarity assumption) of soil properties caused by differences in e.g., topography, climate and vegetation (Elbasiouny et al., 2014; Eldeiry and Garcia, 2010).

The application of machine learning techniques and hybrid spatial models is steadily increasing in environmental and soil science (Keskin and Grunwald, 2018; Li et al., 2011; Veronesi and Schillaci, 2019). Random Forest (RF) (Breiman, 2001) is one of the most applied machine learning technique quickly gained popularity due to its ability to model non-linear relationships on high dimensional data and the possibility to provide variable importance (Liaw and Wiener, 2002; Myles et al., 2004). Therefore, its superiority over other geostatistical and machine learning techniques was well presented in several studies (Subburayalu and Slater, 2013; Wiesmeier et al., 2011) as well as in trend estimation when hybrid approach was considered (Hengl et al., 2015). Hybrid spatial models employ simultaneously both kriging interpolation and auxiliary variables in order to increase prediction accuracy (Odeh et al., 1995; Phachomphon et al., 2010). Recent studies have emphasized outperformances of the hybrid models over other approaches in combination with different scales, usage of varying resolution auxiliary and response variables (Chen et al., 2019; Tziachris et al., 2019).

As regards auxiliary data, digital elevation models (DEM) and multispectral satellite images (MI) have been widely used in DSM (McBratney et al., 2003). Majority of DSM discussions preferably used freely available DEMs and MIs. Together with the modelling approach and auxiliary data, the scale is an essential consideration in DSM. As the soil-forming and environmental factors vary and respond at different scales, selecting DEM and MI resolution is a noteworthy phase. The terrain features can influence the prediction performance of soil properties, the scale at which pedogenic processes occur in the landscape, and whether the DEMs can represent the features in the terrain (McBratney et al., 2003). More recently, there has been growing interest in evaluating how the characteristics of auxiliary variables contribute to the success of DSM. It is commonly accepted that more spatial information the auxiliary variables gain the more accurately they describe environmental conditions (Hengl et al., 2013). Accordingly, several studies focusing on the effect of multiscale terrain analysis (Behrens et al., 2010; Cavazzi et al., 2013), spatial resolution of DEM (Smith et al., 2006) and MI (Xu et al., 2018) confirmed that more detailed covariates lead to more accurate predictions. The overall accuracy of finer resolution DEM is better than that of medium resolution such as SRTM and ASTER DEM, and the rougher and steeper topography is satisfyingly depicted with high resolution DEM (Kramm and Ho, 2019). However, Kim and Zheng (2011) found inverse results that fine-scale topographic information is not always optimal for understanding soil spatial variability. Likewise, the effect of fine and high-resolution MIs delivered by commercial satellites such as SPOT 6, Rapid Eye, Pléiades-1, World View-2/3 have commonly less studied. Sumfleth and Duttmann, 2008 and Xu et al. (2018) confirmed advantages of using higher spatial resolution images for predicting different soil properties in terms of accuracy and error assessment, while others (Samuel-Rosa et al., 2015; Xu et al., 2018) pointed out that the auxiliary variables extracted from higher spatial resolution images may not always produce the most accurate soil prediction.

In this context, discussions for addressing high resolution auxiliary data are limitedly available in mountainous regions. Especially, no study was concerned with DSM using high resolution data in the Caucasus Mountains. Therefore, this study specifically focused on evaluating the potential of high-resolution DEM and MI (Pléiades-1) in predicting soil properties in a local test area in the Caucasus Mountains. The DEM used in this study was a product of the Advanced Land Observing Satellite (ALOS) that has been made publicly available recent years (Alaska Satellite Facility). To produce a model possibly higher accuracy, high-resolution auxiliary data was combined with modelling approaches. Therefore, we tested two spatial hybrid models, Universal Kriging (UK) and Random Forest Kriging (RFK) as well as their associated uncertainty to evaluate and map basic soil properties, soil organic carbon (SOC), soil reaction (pH), particle sizes (sand, silt, and clay), calcium carbonate (CaCO₃) and hygroscopic water content (WC).