Spatiotemporal and spectral analysis of sand encroachment dynamics in southern Tunisia

ABSTRACT

Aeolian processes in drylands often transcend into sand encroachment, a common form of land degradation. Highly reflective desert features, hence sandy areas, often cause spectral confusion, and mapping through remote sensing techniques can be challenging. This work aims at designing an efficient classification method that minimises spectral confusion of desert features, hence two types of sandy areas. Moreover, we employ land cover (LC) change detection over the last 30 years. The extraction and spatiotemporal variations of LC and sand encroachment areas in the Dahar-Jeffara Medenine site (southeastern Tunisia) are assessed by employing Landsat imagery (1984 and 2014), a 30 m digital elevation model of Shuttle Radar Topography Mission (SRTMGL 1 arc second), field data and X-ray diffraction analyses of sand samples. Five new spectral indices were designed and employed in a Decision Tree (DT) classifier for the extraction of 11 LC classes, including two different types of sandy areas. The DT map yielded an overall accuracy of around 89%. Change detection results showed substantial change in several landscape components and an increase of sand units by 29% within the Jeffara-Medenine plain over the last three decades. Geomorphological observations and multi-temporal, spectral and mineral analyses indicate a main, possible in-situ source area of sand.

KEYWORDS:

• Landsat series
• sand encroachment
• post-classification change detection
• spectral indices
• decision tree classification
• X-ray diffraction

Introduction

Interaction of wind and water over time is arguably the key influence upon the contemporary geomorphology and longer-term landscape development of arid regions of the world (Nash, 2000). Aeolian processes in arid and semi-arid areas need constant monitoring, as they often transcend into sandy desertification or sand encroachment. It is a common form of land degradation resulting from a combination of climate change and intensive or unsustainable human activities in ecologically fragile environments (Mihi, Tarai, & Chenchouni, 2017; Ren, Dong, Hu, Zhang, & Li, 2016), being strongly related to regional wind regime (Ge et al., 2016). Therefore, it is a very dynamic phenomenon. Its evolution depends not only on the biophysical conditions of such arid environments but also on the management and exploitation methods of natural resources. Indicators of sand encroachment are represented by depositional landforms, such as sand dunes and erosional landforms due to deflation, with the production of wind depressions or yarding. Sandy desertification determines degradation of land resources; slow recovery of vegetation; and decrease of ecosystem productivity, services and livelihood (Ahmady-Birgani, McQueen, Moeinaddini, & Naseri, 2017; Duan, Wang, Xue, Liu, & Guo, 2014).

Sand encroachment has become a critical issue in the last decades, especially in susceptible drylands, and is one of the most serious environmental problems in central and south Tunisia. In the Jeffara-Medenine-Dahar study area (southeastern Tunisia), the main desertification issues are related to in-situ sand production, deflation and accumulation. In these environments, which are a priori fragile, it was proven that any change in the ecosystem could have immediate and strong impact. In this case, soils are often sandy or developed on fine-grained deposits, and natural vegetation, even if poorly developed, plays and important role in sand control and stabilization. In the last decades, modifications of socioeconomic policies in the past decades led to rapid and extensive changes in land use and management (Talbi, 1997), including:

• an intensive use of natural resources

• accelerated expansion of rain fed and irrigated agriculture (olive trees and annual crops, respectively)

• significant change of agrarian system

• change of vegetation cover and/or land use (from grazing land to arable land, intensive ploughing)

• development of multi-sector activities for income support (urbanisation, services, migration, etc.)

These variations represent the anthropogenic direct and indirect contribution to intensification of land degradation processes, particularly sand encroachment (M Ouessar & Zerrim, 2014; Schiettecatte et al., 2005). In the Jeffara-Medenine-Dahar study area, in particular, previous research showed that the anthropic component contributes to sand encroachment intensification through: (1) the depletion of natural vegetation and change in land use, that is, transforming vast natural areas into arable land, (2) reworking of extended areas, and (3) ploughing method (Msadki, Bouzaida, Taamallah, & Ouessar, 2014; Talbi, 1997).

Remote sensing and GIS are powerful tools for monitoring sand encroachment dynamics, which, integrated with ground observations, provide valid mapping techniques (Brownett & Mills, 2017; Ge, Li, Luloff, Dong, & Xiao, 2015). Several authors have employed these techniques in similar studies in Tunisia, focusing on geomorphological, socio-economic or ecological aspects (Essifi, Ouessar, & Rabia, 2009; Hanafi & Jauffret, 2008; Lorenz, Gasmi, Radebaugh, Barnes, & Ori, 2013; Ouerchefani, Dhaou, Delaitre, Callot, & Abdeljaoued, 2013). However, studies about sand encroachment in the southeastern Tunisia are limited to either local experiences or regional-scale studies, and deeper remarks on driving forces and trends are still required.

In this work, we aim at (1) evaluating and designing an optimal automated mapping technique specific to local context (highly reflective desert features), (2) assessing and defining appropriate spectral indices for supporting decision tree (DT) classifier, and (3) generating land cover (LC) map of 11 classes, including two types of sandy areas. Other specific objectives are also to (1) qualitatively and quantitatively estimate changes of LC and sandy areas, (2) argue source areas of sand in the Jeffara-Medenine-Dahar study area, by analysing and discussing the origin of sand (allochthon vs autochthon) and its increase in the cultivated areas, and (3) assess the accuracy of results and possible limits of the proposed approach. For this purpose, we employed ground truth data, remote sensing and X-ray diffraction analysis data (XRD).

Starting from a multi-temporal analysis of land cover over the last 30 years (1984–2014), this research highlights the increase of sandy areas, focusing on spectral discrimination of different types of sand (desert vs. inner-plain). Five new spectral indices and a DT classifier have been studied to optimise and “force” extraction of spectral information from bands. Thus, two distinct spectral signatures of sand were recognised for the first time in this study. Change detection analysis drove us to quantitatively and qualitatively estimate change rates and, consequently, discuss possible driving factors (Allbed & Kumar, 2013; Allbed, Kumar, & Aldakheel, 2014; Allbed, Kumar, & Sinha, 2014; Fichera, 2012; Vacca et al., 2014; Vogiatzakis & Melis, 2015).

The employed DT classification method is flexible, simple, computationally efficient, allowing the proof and management of rules (Elnaggar & Noller, 2010; Friedl & Brodley, 1997; Ghose, Pradhan, & Ghose, 2010; Melis et al., 2013).

To characterise different source areas, several sand samples were analysed from the mineralogical point of view using XRD (Afrasinei, 2016; Afrasinei et al., 2015; Sadiq & Howari, 2009). These data, supported by a database of field collected information, provided valuable support for the identification of possible source areas (Coudé-Gaussen, 1991; Coudé-Gaussen & Rognon, 2016).

Research background

In the Jeffara-Medenine plain, sand encroachment and gypsum crusts represent the particularities of most surfaces, while the coastline is characterised by the presence of sabkha. In these environments, wind plays an important role in the transportation of both sand and gypsum. This can eventually lead to further land degradation, given the climatic and anthropic factors that intensify these processes.

Socio-economic and political factors heavily condition the dynamics and degree of land degradation phenomena in Tunisia, especially in its southeastern part. This relates to the constantly changing political context since the end of the 19^th century onward (Talbi, 1997). This context refers to the transition from the nomadic or pre-colonial phase to the colonial phase after the end of the 19^th century and subsequent changes to present. The pre-colonial phase was based on nomadic pastoralism, agro-pastoral activities and animal husbandry, lasting until the end of the 19th century, when the semi-nomadic/colonial phase started. This period was based on sedentarisation and urbanisation of settlements and population, new agricultural production systems, rural exodus, privatisation of land and introduction of monetary exchange. This meant the destruction of nomadic social structures and production systems, such as abandonment of lands and pastoral activities in the favour of modern urbanised areas (including services) and sedentarisation, among others. This phase was followed by a period of wrongly implemented “collectivist/colonial” policy after the independence of Tunisia in 1956 (Talbi, 1997). Subsequently, the abandonment of the collectivism in 1969 triggered two main activities that eventually were the most important causes of desertification: a) the return to livestock breeding and implicitly over-pasturing and b) massive land clearing using heavy mechanical means. Therefore, a massive change in land use took place, characterised by a generalised mechanisation, especially for land labour/plough, translated into over-exploitation of resources and mismanagement of land. This series of transitions, changes of land management and use of resources finally compromised the ecological equilibrium of the country. The settling process promoted by the colonial authorities and subsequently by the following independent ones resulted in sparse but concentrated rural settlements, entailing a constant human presence, thus a considerable continuous pressure on the environment.

Therefore, past decades rapid and intense changes in agricultural practices and in land use had direct consequences on land degradation and sediment production within the whole Jeffara plain, which is likely to be easily engaged in wind transportation and deposition in close adjacent areas within the plain (Marini, Melis, Pitzalis, Talbi, & Gasmi, 2008; Talbi, 1997).

Study area

The study site is represented by the Dahar-Jeffara of Medenine area, in the northern part of the governorate of Medenine, southeastern Tunisia (approx. 4000 km², centre area UTM WGS84 zone 32 N, coordinates 613,694.87 E, 3,681,555.11 N). It stretches from the border of the Grand Oriental Erg up to the Mediterranean Sea, crossing over the Dahar plateau and the northern plains of Jeffara (Figure 1). The area is located in a transition climatic zone, extending from (1) the desert/pre-desert Saharan climate located to the south-west of the Dahar plateau, to the (2) arid and semi-arid zone of the Jeffara of Medenine plain, located east of the Dahar escarpment, to a (3) Mediterranean relatively wet and moderate climate in the coastal part. The mean annual maximum temperature is over 36°C. The cumulative annual precipitation is generally below 190 mm (Trabelsi, Abid, Zouari, & Yahyaoui, 2012).

Figure 1. a) Location of Jeffara-Medenine-Dahar study area and the main geomorphological units: A – Jeffara-Medenine plain; B – Jeffara-Medenine escarpment glacis; C – Dahar plateau; D – Grand Oriental Erg/Dahar transition zone; E – Grand Oriental Erg. b) Collected field data, consisting of two entities: linear (GPS track log, orange colour) and GPS point data with description (black triangles). The study area limit is represented in green colour (b). Map projection: UTM – WGS 1984, Zone 32N, EPSG 32,632 (metric coordinates).

Analyses showed that the active winds (>3m/s, capable of displacing sand grains, whereas non-active winds are below this speed) blowing from N, NE, SE are more frequent than those from S, W, and SW (M Ouessar, 2007). Active winds are important and represent 44% of all winds in the Sidi Maklouf area (north of the city of Médenine) and 40.7% in Médenine (Khatelli & Gabriels, 2000). Spring is considered the windiest season, followed by winter and, finally, autumn. During summer time, hot winds blowing from Sahara (sirocco), locally known as “chili”, are dominating.

Geological and geomorphological setting

From the geological point of view, the study site is characterised by the presence of sedimentary sequences ranging from Palaeozoic up to present. According to (Marini et al., 2008), the geomorphology of Dahar-Medenine area can be divided into: a) morpho-structural units and b) fluvial and slope morphological units, subsequently organised in forms and processes. The macroforms correspondent to the morphological, litho-structural and slope units, as identified in the area, are: 1) Fluvial processes and forms; 2) Structural and slope forms; and 3) Aeolian forms, processes and deposits.

The following geological and morpho-tectonic domains can be distinguished (Figure 1(a)):

1. The Jeffara-Medenine Plain domain represents a low-flat-tabular topography covered by Tertiary and Quaternary deposits. It reflects the evolution of the parallel retreat of slopes with intense phenomena of sheet and gully erosion and, locally, a slope replacement. Recent Quaternary is represented by the deposits of silts and the “loess of Matmata”, which are very fine detrital particles. These deposits could represent the main source of fine wind-blown material that can accumulate in deep valleys and wadi alluviums as semi-consolidated deposits (Arras et al., 2014, 2016) or over vast areas in the Jeffara-Medenine plain as non-consolidated deposits. Since these deposits are often found in arable lands, they are constantly reworked (mouldboard ploughing), thus contributing to a further in situ sand production.

Morphological fluvial and slope units are related to runoff processes operated by the ephemeral rivers (wadi) and active in the Plain of Jeffara-Medenine. The Jeffara plain – Medenine sector represents a vast gentle dipping surface towards east, constituted by the superposition of several orders of accumulation glacis, whose genesis and evolution are related to the progressive retreat of the Dahar monoclinal ridge.

• (B)The plain sector closest to the eastern Dahar ridge is represented by a vast system of cones of dejection, mostly carved and inactive, whose coalescence determines, on the whole, the formation of gently sloping surfaces, extensively cut by a hydrological sub-parallel network. In these surfaces, the runoff dynamics such as rill and gully erosion are intense, but slightly significant aeolian forms and processes appear. The barrier of the Dahar ridge limits the winds coming from West, thus hindering these aeolian forms. The result is a minor aridity, allowing a limited development of soil and vegetation and a modest but important agricultural activity for the local people. The glacis area is mainly covered by gypsum and calcareous crust; aeolian dynamics and lack of soil limit agricultural use and vegetation growth, thus restricting grazing, the main local economical resource (Marini et al., 2008).

• (C)The Dahar Plateau and the Jeffara – Medenine Plain (Figure 1). The Dahar Plateau consists of a Lower Triassic – Upper Cretaceous marine-sedimentary succession (overall thickness of about 700 m), characterised by the alternation of sandstone, clay, marl and limestone.

• (D)Grand Oriental Erg/Dahar transition zone

• (E)Grand Oriental Erg

In the central and eastern part of the plain, a wide range of detrital-sedimentary cover can be recognised. These relate to ancient continental erosive phases of late Cenozoic and Pliocene times, characterised by the presence of limestone and gypsum crusts, referring to previous arid and desert climate phases.

Aeolian morphologies are evolving on the current surface, corresponding to sand/silt/dust being trapped by bushes to form phreatophyte mounds and nebkas, or to sand shadows, sand drifts and lee dunes fed by materials coming from large surfaces exposed to the wind. An important aspect is given by the influence of agricultural practices on the soil and the formation of new sediments that can be subsequently carried away by winds. The reworking of land, ploughing methods and cultivation of olive trees create other transportable sediments and even bring to surface crumbled material from gypsum and carbonate sub-surficial crusts (Afrasinei, 2016).

Aeolian deposits represent an invaluable indicator for possible source areas, as their tails point in the predominant wind direction. These surfaces are adjacent to the deposition areas inside the Jeffara-Medenine plain. Since predominant active winds blow from N and NE, these determine an overall sediment transport towards the Grand Oriental Erg. This fact is in accordance with various works dealing with this topic at large/small spatial and time scales (Drake, 1997; Khatelli & Gabriels, 1998, 2000; Swezey, 2001).

Materials and data

Large ancillary datasets, georeferenced field data, XRD and remote sensing data were employed in this study. The ancillary data included geological, geomorphological and topographical maps, agricultural calendars/statistics, and pedological reports mainly provided by local entities that were collaborating on the WADIS-MAR European project (www.wadismar.eu).

The field survey was carried out from April to July 2014, covering an area of about 4000 km². Approximately 400 points were surveyed, and 10 sand samples for XRD analysis were collected from both the inner-plain area and the western Dahar area (including the Grand Oriental Erg limitrophe area). Mobile data collection was done using ArcPad (ESRI) installed on Magellan MobileMapper 6 GPS receiver instrument. Physiographical, geological, geomorphological, soil, vegetation and land use observations were registered for each surveyed point. The collected sand samples were prepared according to standard protocol and analysed using the PANalytical X’Pert Pro X-ray diffractometer for powders (X-ray diffraction Laboratory of the Chemical and Geological Department, University of Cagliari).

The Landsat scenes of June 1984 and June 2014 with less than 10% cloud coverage were downloaded from online USGS platforms (Table 1). A summary of sensors’ characteristics is provided in Table 2. The scenes of June 2014 were employed in generating the “final state” map, being timely consistent with the ground truth campaign and the most adequate in terms of quality among recent years’ images of June. The Landsat scenes of June 1984 (used for the “initial state” map) were the oldest available ones with an acceptable degree of cloud coverage and quality. Two adjacent scenes were necessary for the full coverage of the study area, and thus the closest dates as possible were chosen for each year.

Table 1. Landsat scenes employed in the analysis.

Landsat	WRS path	WRS row	Year	Date	Classification type employment
LT5	190	037	1984	18 Jun	DT classification (change detection – initial state)
LT5	191	037	1984	25 Jun
LC8	190	037	2014	21 Jun	DT classification (change detection – final state)
LC8	191	037	2014	28 Jun

Table 2. Landsat 5 TM and Landsat 8 OLI and TIRS characteristics.

Landsat 4–5	Bands	Wavelength (micrometers)	Resolution (meters)
Thematic Mapper (TM)	Band 1 – Blue	0.45–0.52	30
	Band 2 – Green	0.52–0.60	30
	Band 3 – Red	0.63–0.69	30
	Band 4 – Near Infrared (NIR)	0.76–0.90	30
	Band 5 – Shortwave Infrared (SWIR) 1	1.55–1.75	30
	Band 6 – Thermal	10.40–12.50	120
	Band 7 – Shortwave Infrared (SWIR) 2	2.08–2.35	30
Landsat 8 Operational Land Imager (OLI) and Thermal Infrared Sensor (TIRS)	Band 1 – Ultra Blue (coastal/aerosol)	0.435–0.451	30
	Band 2 – Blue	0.452–0.512	30
	Band 3 – Green	0.533–0.590	30
	Band 4 – Red	0.636–0.673	30
	Band 5 – Near Infrared (NIR)	0.851–0.879	30
	Band 6 – Shortwave Infrared (SWIR) 1	1.566–1.651	30
	Band 7 – Shortwave Infrared (SWIR) 2	2.107–2.294	30
	Band 8 – Panchromatic	0.503–0.676	15
	Band 9 – Cirrus	1.363–1.384	30
	Band 10 – Thermal Infrared (TIRS) 1	10.60–11.19	100
	Band 11 – Thermal Infrared (TIRS) 2	11.50–12.51	100

The latest 30 m-resolution Shuttle Radar Topography Mission digital elevation model SRTM Version 3.0 Global 1 arc second (SRTMGL1, released starting from September 2014) was also employed. Two tiles of the coordinates N33, E010 and N33, E009 were obtained to cover the entire study area. The SRTMGL1 was chosen as it is characterised by a 16 m accuracy (at 90% confidence), and it presents several advantages compared to optical sources-derived DEMs (Arras et al., 2017).

Methodology

In the following paragraphs, the methodological workflow of the study is briefly presented (Figure 2).

Figure 2. Methodological flowchart.

The 1984 and 2014 Landsat scenes were pre-processed (radiometric calibration, atmospheric correction, co-registration check), mosaicked and subsequently spectrally analysed and processed.

The following aspects were examined in the spectral analysis and processing steps: a) spectral behaviour of land features, b) selection of existing spectral indices and design of new ones and c) spatial distribution of sand accumulations related to the study area morphology. Ancillary data was extensively used for the visual analysis of the images and land features of interest. A part of the field data collected in 2014 (213 points out of 410 total) were employed as training data (hereinafter referred to as training set) for spectral analysis and image classification of the time-correspondent Landsat image of June 2014, while the remaining part was used for the validation of results of the same image (hereinafter referred to as validation set). Once trained and constructed, the DT was applied to the 1984 image, and correspondent thresholds were calculated following the same math rules as for the 2014 image.

Successively, spectral indices were selected and tested from literature, and new ones were designed when existing ones failed to extract the expected land cover classes. Finally, the customised DT classifier and post-classification processes (validation through confusion matrix and change detection) were carried out.

Pre-processing

The post-classification comparison is the most widely used and simplest quantitative technique for change detection. The main advantage of post-classification comparison method (unlike image-to-image or pixel-by-pixel change detection (Singh, 1989)) lies in the fact that multitemporal images are classified independently, thereby minimising the problem of radiometric calibration. (Bovolo & Bruzzone, 2015). The land cover variation is obtained in an explicit way for each change. Moreover, some authors (Aleksandrowicz, Turlej, Lewiński, & Bochenek, 2014; Jensen, 2005) state that it does not even require radiometric normalisation or atmospheric correction. Nonetheless, accuracy depends on pre- and post-classifications processing, and every error in the individual maps will be present in the final change detection product.

For both years, 1984 and 2014, level 1 precision terrain corrected products were obtained as courtesy of the United States Geological Survey (USGS) and downloaded from USGS Earth Explorer online platform. In 2016, the USGS started reorganising the Landsat archive into a formal tiered data Collection structure. This data Collection structure ensures that Landsat Level-1 products provide a consistent archive of known data quality to support time-series analyses and data “stacking” while controlling continuous improvement of the archive and access to all data as they are acquired (https://landsat.usgs.gov/landsat-collections). The calibration of the Landsat sensors is supported by preflight, post launch-onboard, and ground reference data, undertaken by USGS (https://landsat.usgs.gov/calibration). The absolute radiometric calibration of the Landsat 8 Operational Land Imager (OLI) was updated in February 2017 to reflect changes seen in the instrument’s response to the on-board calibrators indicating a change in the instrument (Micijevic, Haque, & Mishra, 2016).

Even Landsat 4–5 TM scenes are processed through the Level 1 Product Generation System (LPGS) (https://lta.cr.usgs.gov/TM). The downloaded scenes, including the TM ones, belong to TIERS 1 category that meet formal geometric and radiometric quality criteria. These contain the highest quality Level-1 Precision Terrain (L1TP) data considered suitable for time-series analysis. It consists of corrected data that have well-characterised radiometry and are inter-calibrated across the different Landsat instruments. Calibration and validation are important aspects of any remote sensing system. The georegistration of Tier 1 scenes is consistent and within prescribed image-to-image tolerances of ≦ 12-meter radial root mean square error (RMSE).

The images used in this study were subset to the study area extent, radiometrically normalised to obtain top of atmosphere reflectance (Chander, Markham, & Helder, 2009), and, subsequently, dark object subtraction was applied for atmospheric correction, thus obtaining surface reflectance. Landsat 8 OLI scenes were converted from digital numbers (DN) to surface reflectance through rescaling according to the instructions of the provider (https://landsat.usgs.gov/using-usgs-landsat-8-product). As each product employed ground control and relief models (as delivered by the provider, namely USGS), geometric correction was not performed. Given the characteristics of our study area and the fact that the provider already performed a thorough correction, no topographical correction was applied since L1TP (Collection 1,) and L1T (Pre-Collection) was also orthorectified using ground control points and a digital elevation model (DEM) to correct for relief displacement (https://landsat.usgs.gov/geometry). All images were manually verified for coregistration using the Image-to-Image registration module of ENVI image processing software; the accuracy of the co-registration was verified using 20 checkpoints, giving a positional accuracy within an RMSE inferior to 0.5 pixel (Byun, Han, & Chae, 2015).

Classification method

Spectral analysis

Since both field data and Landsat 8 OLI imagery were acquired in June 2014, these two datasets were integrated in order to understand and assess the spectral discrimination between the desert sand and the inner-plain sand and to identify suitable spectral techniques for their delineation: band math and band transformation. These analyses were also supported by the geomorphological considerations on the aeolian forms and processes of the area previously studied by (Marini et al., 2008).

The spectral responses of land features retrieved from the Landsat mosaic of two Landsat 8 scenes acquired on the 12^th and the 21^st of June 2014 were analysed. Two different spectral responses of two types of sand (the inner-plain and the Grand Oriental Erg or desert) could be visually observed (Figure 3 (a–c)). However, the reflectance values of these two types of sand were very similar when inspected pixel by pixel. Two pairs of regions of interest of the two sand types, inner-plain and desert, respectively, were engaged for spectral separability computation (Jeffries-Matusita). The resulting value (1.03) indicated low separation. Further to checking their extraction, a series of first tests of various appropriate supervised classification methods were conducted and misclassification issues emerged. In addition, the two sands regions of interest (ROIs) sets were also compared to other land features ROIs by computing spectral separability. Several values were very close to 1. The reason for these similarities is argued in the discussion section. The average spectra of the 213 training ROIs are presented in Figure 4. The number of pixels and standard deviation of each ROI are presented in Figure 5. ROIs standard deviation statistics are presented in Table 3. The number of pixels of each ROI represents the number of point data collected in the field, belonging to the training set.

Table 3. Standard deviation statistics of training ROIs.

	Hydrosaline	Bare land	Carbonate-rich	Water	Aeolian	Desert	Saline	Steppe	Green veg	Herbaceous	Polygenic
Band 1	0.0089	0.0059	0.0057	0.0109	0.0076	0.0083	0.0198	0.0088	0.0111	0.0052	0.0024
Band 2	0.0132	0.0067	0.0073	0.0129	0.0091	0.0091	0.0243	0.0116	0.0151	0.0061	0.0025
Band 3	0.0241	0.0089	0.0125	0.0182	0.0128	0.0102	0.0312	0.0207	0.0257	0.0088	0.0058
Band 4	0.0394	0.0177	0.0206	0.0291	0.0202	0.0198	0.0346	0.0328	0.0503	0.0159	0.0158
Band 5	0.0602	0.0240	0.0224	0.0276	0.0213	0.0225	0.0351	0.0379	0.0366	0.0205	0.0185
Band 6	0.0861	0.0323	0.0301	0.0369	0.0235	0.0317	0.0352	0.0368	0.0654	0.0208	0.0155
Band 7	0.0850	0.0352	0.0282	0.0373	0.0261	0.0392	0.0324	0.0357	0.0687	0.0258	0.0261

Figure 3. Spectral analysis example. (a) False colour composite (FCC) of Landsat 8 OLI image of June 2014, bands 742 RGB. Field-collected points are represented with green circles (213 points used to generate training ROIs,) and red circles (validation ROIs, 228 points). Map projection: UTM – WGS 1984, Zone 32N, EPSG 32,632. (b) Sand encroachment example in the Jeffara plain area emphasised by the multivariate analysis technique based on principal components analysis (PCA) applied to Landsat 8 OLI of June 2014. Inner plain sandy areas appear ones yellowish with green and orange hues in the RGB 123 (1,2,3 are the first 3 principal components extracted from all seven Landsat 8 OLI bands). (c) Sand encroachment example in the Dahar area emphasised by the multivariate analysis technique based on principal components analysis (PCA) applied to Landsat 8 OLI of June 2014. Desert sandy areas appear with more enhanced yellow hues that the inner-plain in the RGB 123 (1,2,3 are the first 3 principal components extracted from all seven Landsat 8 OLI bands).

Figure 4. Average spectra of training ROIs.

Figure 5. Standard deviation of training ROIs.

Consequently, several image analysis methods were employed in order to deepen the possibility of their separation and to identify optimum bands for their extraction. Horizontal, vertical and Z profiles and 2D scatter plots analysis gave insights on the optimal bands which could be used for the separation of the two different types of sand (Elnaggar & Noller, 2010; Immordino & Melis, 2008; Melis & Pilloni, 2011; Pilloni, Melis, & Marini, 2010). Band transformation and spectral enhancement techniques (among which principal component analysis, Tasseled Cap, and Knepper ratios performed best (Mia & Fujimitsu, 2012)) were used in order to better distinguish the two sand surfaces between them and also to delineate them from the other land features.

Decision tree classifier

Traditionally, classification tasks are based on statistical methodologies such as Minimum Distance-to-Mean, Maximum Likelihood and Mahalanobis Distance Classification (Simard, Saatchi, & DeGrandi, 2000). These classifiers are generally characterised by an explicit underlying probability model and are ideally suited for data in which its distribution within each of the classes can be assumed to follow a normal distribution in a multispectral space (Srimani & Prasad, 2012). A decision tree, as a predictive model, is a type of multistage classifier that can be applied to a single image or a stack of images. Decision trees that are used commonly in the remote sensing community include the automatic machine learning univariate algorithms such as C4.5 and C5.0/See5 (Pandya & Pandya, 2015; Quinlan, 1993) or multivariate algorithms such as Quick Unbiased Efficient Statistical Trees (QUEST) or classification and regression trees (CART) (Loh, 2011; Wang, Liu, Jiang, Zhang, & Zhao, 2012). However, in this study, we implemented the decision tree manually to assure maximum user control and structure manipulation. ENVI® Classic’s (ITT Vis Exelis) interactive decision tree tool is a useful way to explore data, or to find areas in data that fulfil certain criteria. The user does not need to start by having a tree that has been generated outside of ENVI Classic (thus it is free) in order to learn useful information about the dataset using the decision tree tool (Exelis Visual Information Solutions, 2013). Decision trees do not make any a priori assumption about data distribution, and hence they are essentially non-parametric in nature. The DT approach was applied due to its proven high accuracy in this these particular arid environments of very reflective, spectrally similar features, being a supervised, non-parametric, multi-stage classification method (Afrasinei et al., 2015; Melis et al., 2013; Otukei & Blaschke, 2010; Pal, 2012; Pal & Mather, 2003; Rogan et al., 2008). It is very effective and useful for the remote sensing community, especially for LULC classification. This classifier is intuitive, simple, flexible, and efficient in computing non-normal, non-homogenous and noisy data as well as non-linear relation between features and classes, missing value, and both numeric and categorical inputs (Otukei & Blaschke, 2010; Pal & Mather, 2003; Rogan, Franklin, & Roberts, 2002).

Classification assessment

Conventional classification methods were assessed through literature review and also by running some tests. The example of the traditional pixel-based image classification method, Maximum Likelihood (hereinafter MaxLike) classification, is described below. It was applied to the Landsat 8 OLI image of 2014 using training ROIs derived from the training set field data of the same period, as described in the section “Materials and data”.

The MaxLike classification was used due to its acknowledged efficiency among supervised pixel-based algorithms for land cover mapping (Hogland, 2013; Julien, Sobrino, & Jiménez-Muñoz, 2011; Simard et al., 2000). We managed to extract 15 LC classes. However, confusion matrix using the validation set (see the section “Materials and data”) showed an overall accuracy of 31.91% and a Kappa Coefficient of 0.29. The reason for these low values is probably related to high spectral similarity between classes, mixed pixel issues, the high number of classes set in the training areas and the heterogeneity within each class. Furthermore, parametric and single stage classifier often present this shortcoming, not being able to deal with non-linear data (Hogland, 2013; Perea Moreno & Meroño De Larriva, 2012; Rogan et al., 2008; Srimani & Prasad, 2012).

Designing new spectral indices

To extract the 11 LC classes, including sandy areas, we used band ratios and various vegetation or mineral indices as decision nodes, either extracted from literature or proposed in this study. Threshold values of decision nodes were retrieved through a histogram-based thresholding technique according to (Singh, 1989; Wu et al., 2014; Yang & Lo, 2000), and standard deviation and mean statistics were jointly used with skill-based visual observation of each index image. A threshold was tentatively applied according to the mean ± standard deviation (García, Contreras, Domingo, & Puigdefábregas, 2009) to examine whether it can divide the two areas of interest for each node (e.g. vegetated from non-vegetated areas) while compared with the natural colour composite or adequate false colour composite of the original images. The threshold values are presented in Table 4.

Table 4. Threshold values of the Landsat images of both 1984 and 2014.

Index image	Year	Mean	Standard deviation	Threshold value
SMI	1984	0.636	0.042	0.720
	2014	0.596	0.044	0.684
NDVI	1984	0.145	0.014	0.152
	2014	0.163	0.021	0.174
NDWI	1984	−0.081	0.025	−0.106
	2014	−0.085	0.020	−0.105
NDWIusgs	1984	−0.145	0.014	−0.159
	2014	−0.163	0.021	−0.184
IRIswir1	1984	0.784	0.058	0.842
	2014	0.858	0.054	0.912
S2	1984	−0.399	0.037	−0.437
	2014	−0.529	0.038	−0.567
WR	1984	0.747	0.022	0.790
	2014	0.720	0.030	0.780
MI	1984	0.018	0.004	0.022
	2014	0.014	0.004	0.018
1/5	1984	0.273	0.029	0.000
	2014	0.188	0.024	0.000
HSI	1984	1.741	0.149	1.890
	2014	1.435	0.137	1.571
IRInir	1984	1.466	0.058	1.495
	2014	1.528	0.063	1.560

Vegetation, water, mineral indices and band ratios applied in similar studies and bio-geographical areas were constructed or retrieved from literature, choosing the relevant ones for the purpose of this study (Bachir, Essifi, Zerrim, & Ouessar, 2013; Boulghobra, 2016; Hadeel, Jabbar, & Chen, 2011; Ouerchefani et al., 2013; Pandey, Rani, Srivastava, Sharma, & Nathawat, 2013; Sadiq & Howari, 2009). Once vegetation cover and water bodies were extracted, a vast area of bare soil/land remained unclassified. Since our study area is characterised by highly reflective surfaces, with a main mineral component, our objective was to correctly separate the various types of these “mineral surfaces”. From the land cover/land use point of view, these surfaces are classified as “bare lands”. However, from the lithological/mineral point of view, these classes can be separated on the basis of their main mineral component: carbonate-rich areas, halite saline areas, high quartz – content sandy areas, clay minerals – rich sandy areas and polygenic (usually alluvial deposits, slope deposits, mineral soils). Thus, the priority was to identify or define specific mineral spectral indices in order to discriminate these surfaces correctly, implying that the indices proposed in this study were particularly designed to meet these needs.

Five new indices were constructed based on choosing optimal band pairs/groups that presented the highest covariance in terms of spectral information of each feature of interest. To identify optimum bands for indices construction, the existing literature also proposes methods like Optimum Index Factor (OIF) or correlation analysis between measured spectral reflectance and satellite data (Sidike, Zhao, & Wen, 2014; Yu et al., 2010). Based on the principle of absorption bands of each land feature in specific wavelengths, optimal bands were identified. The training set ROIs (corresponding to Landsat 2014 image) were used to visually evaluate the performance of each index and in verifying the type of extracted feature (by overlapping ROIs to each index image).

Band ratios and indices were constructed through mathematical operations of bands (sum, subtraction, multiplication, division, square root, exponent) of highly decorrelated bands in order to discriminate features of interest with the highest possible separability. For Salt Minerals Index (SMI), Mineral Index (MI), InfraRed Index – Short Wave InfraRed 1 (IRI_SWIR1) and InfraRed Index – Near InfraRed (IRI_NIR) indices (Table 5), exponential and square root functions were used to force the emphasis of extreme values, helping to better discriminate areas of similar mineral content (Afrasinei, Melis, Buttau, Arras, et al., 2017; Allbed & Kumar, 2013; Allbed et al., 2014).

Table 5. Indices used in the DT classifier applied to the June 2014 Landsat image, Jeffara-Medenine-Dahar area.

Parent nodes-decision	Band math	Feature extraction	Index	Reference
NDVI		Vegetation/No vegetation	Normalised Difference Vegetation Index	(Gitelson, Peng, Arkebauer, & Schepers, 2014)
NDWI		Steppe/green and herbaceous vegetation	Normalised Difference Water Index	(Gao, 1996; Kross, McNairn, Lapen, Sunohara, & Champagne, 2015)
NDWI USGS		Green vegetation/herbaceous vegetation	Normalised Difference Water Index – USGS	(McFeeters, 1996), also known as Normalised Difference Salinity Index (NDSI) (Khan, Rastoskuev, Sato, & Shiozawa, 2005)
WR		Water bodies	Water Index	derived from (Kross et al., 2015; van der Meer et al., 2012)
SMI		Highly saline areas	Salt Minerals Index	proposed for this study
MI		Bare land – carbonate rich	Mineral Index	proposed for this study
IRI_SWIR1		Desert sand – high quartz content	InfraRed Index – Short Wave InfraRed 1 (TM band 5)	proposed in this study
IRI_NIR		Bare land/bare land (polygenic high mineral content material, usually slope deposits, skeletal mineral soils)	InfraRed Index – Near InfraRed (TM band 4)	proposed in this study
S2		Aeolian sand (silty sand, clay minerals, loess component)	Salinity Index 2	(Allbed & Kumar, 2013; Allbed et al., 2014)
HSI		Hydrosaline areas (sabkha)	Hue Salinity Index	proposed in this study
Ratio Blue/SWIR1		Mask (no data)	Ratio Blue/SWIR1	(Melis et al., 2013)

X-ray diffraction analysis

Main minerals identification and quantity were achieved by XRD analysis based on the Reference Intensity Ratio (RIR) method. Spectral behaviour of samples, as retrieved from remote sensing data, is probably associated with the different proportions in quartz (main mineral component) and its association with clay minerals, calcite and gypsum (in the inner plain areas) or mainly with dolomite and calcite (in the desert area).

Change detection

Change detection analysis was performed through the post-classification approach, defined as the comparative analysis of individually produced classifications of different dates (El-Hattab, 2016; Latifovic & Pouliot, 2006; Nielsen, 2007; Woo & Do, 2015). Usually, this method consents the use of complex legends and, consequently, a more detailed analysis of land cover changes (Serra, Pons, & Saurí, 2003). High accuracy of each classification is required by the employment of this method. However, this is not always achieved, especially when a legend containing several agricultural categories (probably similar in terms of spectral signature) is required. In this case, considering fragmentation of the landscape, planimetric accuracies and pixel size, we considered that this change detection approach is proper, but other methods will be assessed in future research (Canty & Nielsen, 2012; Singh, 1989). The DT map of 1984 was used as initial-state image and the DT map of 2014 as final-state image. Post-classification change detection tool of ENVI software (ITT VIS Exelis) was applied in order to obtain the statistical input for quantitative change estimation of each LC class.

Results and discussion

DT classification and XRD analysis

Pre-existing indices selected from literature and the new ones used for the DT classifier are reported in Table 5. The DT nodes hierarchy is illustrated in Figure 6. The description of the resulting 11 LC classes is provided in Table 6. The proposed indices SMI, MI, IRI_SWIR1, IRI_NIR and Hue Salinity Index (HSI) provided an accurate extraction of mineral surfaces. Through the employment of these new indices and the pre-existing ones retrieved from literature, we managed to extract a higher number of land cover classes (11) than the average low number of classes (commonly between 4 and 7 classes) usually obtained through traditional classification schemes at local-regional scales (Avelar & Tokarczyk, 2014; Ceccarelli, 2013; Li et al., 2014; Marconcini, Fernandez-Prieto, & Buchholz, 2014; Nutini, Boschetti, Brivio, Bocchi, & Antoninetti, 2013). Furthermore, the discrimination of various mineral surfaces was consistent with the purpose of separating two different types of sandy areas.

Table 6. Description of the resulting classes of the DT classification.

Green vegetation	Oasis vegetation, mainly palm groves, fruit trees plantation, types of small trees and tall shrubs rich in biomass or chlorophyll
Herbaceous vegetation	Annual crops, small natural herbaceous vegetation, small shrubs
Steppe	Typical dry shrub vegetation, woody correspondent to mountainous and piedmont areas
Carbonate-rich areas	Areas with a high carbonate component, limestone crust, outcropping limestone
Polygenic deposits	Deposits correspondent to alluvial fans, recent alluvial deposits, piedmont and glacis accumulations, with a strong clay, sand and coarse materials component
Bare land	Land with no vegetation cover and of no land use, bare rocks
Highly saline areas	Areas rich in salt minerals components
Hydrosaline areas	Areas that present moderate salinity and high moisture content
Water bodies	Water areas of natural or artificial origin
Inner-plain sand	Areas of sand accumulations (dunes) or veils/sheets found within the Jeffara-Medenine plain, covered by sparse or very low vegetation; it can present typical semi-arid phytogenic aeolian coppice dunes (nebkas)
Desert sand	Areas of sand accumulations/dunes in pre-desert areas and/or at the border of the Grand Oriental Erg, usually without vegetation cover

Figure 6. Decision Tree classification nodes.

The DT maps of 1984 and 2014 are presented in Figures 7 and 8, respectively. Mapping results show a clear visual difference between sandy areas mainly found westwards of the Dahar domain (represented with yellow colour in Figures 8 and 9), denoted as “desert sand”, and those within the Jeffara-Medenine plain, denoted as “aeolian sand” (represented with orange colour in Figures 7 and 8). Examples of field observations on aeolian sand, hence sand encroachment areas in the Jeffara plain, and possible on-site indicators are given in Figures 9 and 10. The same points are overlain on the false colour composite (FCC) of Landsat 8 OLI June 2014 image bands 742 RGB (typical combination for geological/lithological observations) (Figure 11) and on the DT map of June 2014 (Figure 12).

Figure 7. Decision Tree classification applied to June 1984 images, Jeffara-Medenine-Dahar area. Map projection: UTM – WGS 1984, Zone 32N, EPSG 32632 (metric coordinates).

Figure 8. Decision Tree classification applied to June 2014 images, Jeffara-Medenine-Dahar area. Map projection: UTM – WGS 1984, Zone 32N, EPSG 32632.

Figure 9. Example of sand encroachment area in the Jeffara-Medenine plain (360-degree-view panoramic photograph).

Figure 10. Example of aeolian erosion indicators: The roots of Olive trees indicate previous level of soil. Field observations were carried out in collaboration with the technicians from the Eremology Laboratory of the Institute des Regions Arides (IRA, Medenine, Tunisia), one of the WADIS-MAR project partner.

Figure 11. Figures 9 and 10 location represented on the 742 RGB FCC of Landsat 8 OLI image of June 2014. Map projection: UTM – WGS 1984, Zone 32N, EPSG 32632.

Figure 12. Figures 9 and 10 location represented on the DT map of June 2014. Map projection: UTM – WGS 1984, Zone 32N, EPSG 32632.

The spectral separation of the two sand types may have been possible mainly due to the different mineralogical composition of the two sets of samples, as shown by the XRD analysis. The XRD results are presented in Figures 13 and 14, where the corresponding locations of the points are highlighted in the same colours. The mineralogical composition seem to be strongly connected to the spectral separability of the two types of sandy surfaces. Despite that, the XRD data are not conclusive, and further insights are needed, as well as more samples. An important aspect observed during the analyses is the contribution of granulometry and petrographic characteristics in separating the inner-plain samples from the desert ones. This aspect will be deepened in a future analysis as well as collecting more sand samples from both sectors (inner-plain and desert).

Figure 13. XRD analysis and percentages (semi-quantitative R.I.R. method): The inner-plain sand samples are grouped in the red rectangle and the Grand Oriental Erg sand samples are clustered in the yellow rectangle, in correspondence with the same separation in Figure 14.

Figure 14. Location of XRD sampling points within the study area, displayed on the FCC 742 RGB, where the inner-plain sand samples are differentiated from the Grand Oriental Erg sand samples through the red and yellow rectangles, respectively, in correspondence with the same separation in Figure 13. Map projection: UTM – WGS 1984, Zone 32N, EPSG 32632.

Accuracy assessment

The validation of the DTC maps showed an overall accuracy of around 89% and a kappa coefficient of 0.87. The confusion matrix statistics are presented in Table 7 and accuracy percentages, in Table 8. The accuracy assessment of the DT map and field data confirm the validity of the methodological approach and the need of dedicated spectral techniques for challenging extraction of classes known to be prone to spectral confusion (Fares & Philip, 2008). In fact, the only classes that presented an accuracy below 80% are carbonate-rich areas (68.18% producer accuracy) and polygenic areas (76.92% user accuracy). Indeed, the confusion matrix show that the most important spectral confusion issues, though slight, occurred in these class pairs (a) polygenic areas and bare land, (b) highly saline areas and carbonate-rich areas class, (c) carbonate-rich areas and desert sand and (d) steppe and herbaceous. These results are explained and confirmed by existing literature, among which we mention:

• The mineralogy of these features, in particular for pair (b), whose fundamental anion groups (HOH, OH⁻, CO₃²⁻, SO₄²) (Afrasinei et al., 2017; Clark, 1999) the internal vibration modes and thus spectral response.

• In the case of pair (a), (b) and (c), (G. Metternicht & Alfred Zinck, 2008; G. I. Metternicht & Zinck, 2003; Xu, 2010) explained the spectral confusion between light-coloured bedrock or high-mineral content areas (salt crusts, sebkhas, rocky desert “reg”, limestone outcrops, sand dunes, sand sheets) and light soil, and between salty crusts and bright silt loam structural crusts.

• The pair (d) can be explained by the fact that even though steppe contains a high number of woody and typical dry shrub vegetation, and the herbaceous class refers to annual crops, natural annual herbaceous vegetation, small shrubs (see Table 5), areas where a cereal crop had been recently harvested and the stubble and straw were left covering the field, can contribute to this confusion, as confirmed by (Pellikka & Clark, 2009).

• The case of pair (c) is also reposted in (Fares & Philip, 2008) mentioning that separating quartz minerals from carbonate minerals is challenging. However, this can be also explained by the fact that in the western part of the area, there are vast areas where thin sheets or layer of blown desert sand overlay irregularly the reg (rocky desert often consisting of local bedrock, in our case limestone, calcareous concretions), thus often leaving small windows of the underneath feature to be detected by the sensor.

However, other factors must be taken in consideration in order to explain our results. One regards the salt and pepper effect that a substantial percentage of the area has, since we did not apply any post-classification refinement (Latifovic & Pouliot, 2006), thus maintaining, as far as possible, the original information contained in the images. These areas are numerous where complex vegetation occurs, where natural vegetation is mixed with tabias and jessour crops (both in the mountain and lain areas), with interspersed sand sheets and nebkas and patches of sebkhas (in the coastal areas).

Table 7. DT confusion matrix.

	Ground points ROIs (Pixels)
DT map classes	Class 10	Class 11	Class 12	Class 1	Class 2	Class 3	Class 4	Class 5	Class 6	Class 7	Class 8	Class 9	Total
Class 10 Herbaceous vegetation	17	0	0	0	0	0	0	0	0	0	0	2	19
Class 11 Mask	0	10	0	0	0	0	0	0	0	0	0	0	10
Class 12 Polygenic	0	0	15	0	1	0	0	1	0	0	0	0	17
Class 1 Hydrosaline	0	0	1	9	0	0	1	0	0	0	0	0	11
Class 2 Bare land	1	0	3	1	48	1	0	1	0	0	0	0	55
Class 3 Carbonate-rich areas	0	0	2	0	1	19	0	0	2	2	0	0	26
Class 4 Water bodies	0	0	0	0	0	0	4	0	0	0	0	0	4
Class 5 Aeolian sand	0	0	0	0	0	0	0	19	2	0	1	0	22
Class 6 Desert sand	0	0	0	0	0	1	0	1	13	0	0	0	15
Class 7 Highly saline	0	0	0	0	0	2	0	0	0	11	0	0	13
Class 8 Steppe	0	0	1	0	0	0	0	0	0	0	12	0	13
Class 9 Green vegetation	2	0	0	0	0	0	0	0	0	0	0	21	23
Total	20	10	22	10	50	23	5	22	17	13	13	23	228

Table 8. Accuracy assessment.

	Commission	Omission	Commission	Omission	Prod. Acc.	User Acc.	Prod. Acc.	User Acc.
DT map classes	Percent	Percent	Pixels	Pixels	Percent	Percent	Pixels	Pixels
Class 10 Herbaceous vegetation	10.53	15	2/19	3/20	85	89.47	17/20	17/19
Class 11 Mask	0	0	0/10	0/10	100	100	10/10	10/10
Class 12 Polygenic	11.76	31.82	2/17	7/22	68.18	88.24	15/22	15/17
Class 1 Hydrosaline	18.18	10	2/11	1/10	90	81.82	9/10	9/11
Class 2 Bare land	12.73	4	7/55	2/50	96	87.27	48/50	48/55
Class 3 Carbonate-rich areas	23.08	13.04	6/26	3/23	86.96	76.92	20/23	20/26
Class 4 Water bodies	0	20	0/4	1/5	80	100	4/5	4/4
Class 5 Aeolian sand	9.09	9.09	2/22	2/22	90.91	90.91	20/22	20/22
Class 6 Desert sand	6.67	17.65	1/15	3/17	82.35	93.33	14/17	14/15
Class 7 Highly saline	7.69	7.69	1/13	1/13	92.31	92.31	12/13	12/13
Class 8 Steppe	7.69	7.69	1/13	1/13	92.31	92.31	12/13	12/13
Class 9 Green vegetation	8.7	8.7	2/23	2/23	91.3	91.3	21/23	21/23
Overall Accuracy	(202/228) 88.5965%
Kappa Coefficient	0.8713

Another aspect concerns possible interference by vegetation cover (olive trees cultivations that overlap all the inner-plain sandy areas, distributed following a regular pattern). Furthermore, fine aeolian gypsum deposits are often blown, deposited and mixed with the sandy component, and thus the sand encroachment areas have a substantial gypsum component as shown by the XRD analyses.

Another aspect to be considered is that ROIs were not grown or buffered, being derived from point data, thus making validation challenging for us. This regards also the issue of mixed-pixels and the fact that, given a pixel of 30m by 30m, a GPS receiver may present low accuracies in remote areas such as the present one, especially in the mountain parts.

Finally, this procedure is a priori subjected to some errors due to the in-field recognition (Gomarasca, 2009), the field of view of the observer, and the evaluation subjectivity of the user when estimating the predominant land cover facet when collecting a field point information and eventually determining the correspondent land cover class.

Change detection

The results and statistics of the change detection analysis indicated an increase of the inner-plain sand by 29% and a decrease in desert sand by 20%, as shown in Figure 15. This analysis demonstrated substantial change in several landscape components in the last decades, particularly concerning sand encroachment expansion within the Jeffara-Medenine plain by 29%. Spectral and mineral analyses indicate that the Jeffara-Medenine plain may have a main, possible in-situ source area of sand. A decrease of 20% of the “desert sand” class areal extent would be consistent with these findings. In fact, transport of sand from western directions to eastern parts is not significant enough to create vast areas of sand encroachment in the Jeffara-Medenine plain (Bachir et al., 2013; Khatelli & Gabriels, 2000; M. Ouessar, 2010).

Figure 15. Change detection statistics obtained from the DT maps of 1984 and 2014, Jeffara-Medenine-Dahar area.

The geomorphological and soil/vegetation parameters define the current potential, although limited, of the agricultural and livestock practices, despite the high incidence of human activity. The strong predisposition to the onset of arid conditions suggests a thorough planning of the agricultural and livestock husbandry activities. The geomorphological interpretation of the territory (100 km in length) shows the presence of three distinct morpho-climatic systems: semi-arid, arid and desert. Transition zones between adjacent systems are very limited. The rapid passing from a morpho-climatic system to another can be related to physiographic causes, for example, orography control over the distribution of rainfall and influence of local wind regime that participate to generate aeolian forms, sand encroachment and dune morphologies. This aspect, apart from having a key role over surface morphologies and processes, highlights the importance of the assessment of the desertification problem, highly contextualised to the morphological and orographic setting.

The methodological approach and the use of Landsat 8 data has provided results that show a clear distinction between the spectral reflectance of the desert sand and the inner-plain sandy surfaces. Their spectral signatures extracted from remotely sensed data are separable, confirming XRD analysis of the sand samples, which show slightly different mineralogical composition. We can infer that by means of the remote sensing techniques, integrated with ground truth information and the X-ray diffraction, it is possible to discriminate the two types of sandy surfaces and to assess spatial and temporal variations of each class. The XRD analysis showed that the inner plain samples are characterised by the presence of quartz (70–80%) mainly associated with clay minerals as well as gypsum and calcite, whereas desert sand revealed a greater proportion of quartz (usually over 90%) mainly associated with calcite and dolomite.

The results gave a reading of the dynamics of forms and processes of aeolian morphologies and their connection to either natural or anthropic factors. However, a clear distinction should be made between factors which influence a process and driving factors that trigger that process. The indicated changes are likely to be related to increased anthropic pressure and settlement, agricultural policies and national development strategies, in accordance with other works conclusions. Considering both factors, it would be accurate to say that the wind is a vector, not a trigger. The anthropic factor can be a trigger, related to the type and the degree of agricultural practices, change in land use, overexploitation of resources. We may argue that the sand encroachment phenomenon present within the Jeffara-Medenine plain has an autochthone, internal main source area. This inner-plain source area is also shown by predominant wind circulation and macro-scale studies of current aeolian processes. During the three-months internship at the “Institute des Regions Arides” headquarters at El Fje, Medenine, the issue of in-situ sand production has been extensively discussed with local researchers, field technicians and local farmers. Furthermore, field observations, past agricultural reports on pastures extent and land use change helped understand that around 30 years ago, the current sandy areas of the Jeffara-Medenine plain and costal area were mainly covered by natural vegetation. Sandy areas increased their expansion in the last 30 years concomitant with the transformation of those lands in vast Olive trees plantations, which have been constantly ploughed ever since (mouldboard).

Conclusion

The customised and replicable methodological approach proposed in this study, involving remote sensing techniques and field data analysis, proved to be efficient for monitoring spatio-temporal dynamics of sand encroachment phenomena. Using a Decision Tree classifier, both new and already existing spectral indices were employed for the extraction of eleven land cover classes. The new indices were specifically designed for the extraction of high-mineral content surfaces, among which two different sand types: desert (Grand Oriental Erg) and inner plain (Jeffara-Medenine). The classification results presented an overall accuracy of approx. 90%. The Decision Tree classifier has proven to be more flexible and adequate for the extraction of sandy areas and major land cover types, as it allows multi-source information and higher user control.

The results of 1984–2014 change detection analysis show an increase of inner-plain sandy areas and a decrease of desert sand. This is consistent with a possible in-situ production of sandy material and would be related to extensive change in land use over the past decades and intensive agricultural practices. Therefore, the assessment of land cover change through remote sensing and GIS is a valuable and efficient approach for delivering comprehensive findings to stakeholders and decision makers. A better evaluation of the sustainable strategies against desertification problem and taking effective measures of soil degradation reduction in such arid and semi-arid areas can be achieved through the proposed approach.

In the Jeffara-Medenine plain and costal area, land use change, overgrazing, expansion of cultivation and use of tillage or mouldboard ploughing instead of deep mole plough may have a great contribution in fine particles/sand production, potentially transportable by active winds (>3m/s). Mouldboard or turn plough is commonly applied in less fertile areas. To grow crops regularly, the soil must be turned to bring nutrients to the surface. In the Jeffara-Medenine plain area, mouldboard ploughing is practiced on small farming patches for crops and trees (private family cultivated perimeters, subsistence agriculture) but also on extended surfaces for Olive trees plantations (in the whole Jeffara coastal plain and inner, east of Dahar plain, mainly in areas already affected by sand encroachment). Therefore, the intensive and wide-ranging agricultural activities in the Jeffara-Medenine plain give a major contribution to the additional formation of dust, silt and sandy sediments through ploughing, re-working and intense mechanical labour of land in the whole Jeffara-Medenine plain. Sandy areas expansion in the Jeffara plain, massive land use change and decrease in natural vegetation cover infer intense anthropic pressure and agricultural practices.

The importance of this research is given by the likelihood that driving factors not related to climate change or the environment have strongly marked the landscape functionalities in the last decades, such as the anthropic ones. The necessity of this study is given by the need of constant monitoring through repeatable mapping of ongoing land degradation processes, which is indispensable for a proper land management designation and the implementation of adequate agricultural practices. In its turn, these operations require a good understanding of the dynamics and processes of these phenomena. Therefore, the methodologies and tools employed in this study contribute to a knowledge of desertification and land degradation processes, and hence a better territorial planning.

Acknowledgements

This study was developed within the frame of the doctoral research (2013–2016) of Gabriela Mihaela Afrasinei, financed by the Italian Ministry of Education, Universities and Research (MIUR), in “Soil preservation, environmental vulnerability, and hydrogeological protection” of the Doctoral school in Environmental and Land Sciences and Engineering, University of Cagliari, Italy. We wish to acknowledge the NRD team of the University of Sassari, the WADIS-MAR team for the scientific and financial support and the possibility of conducting the current research. We also thank the TeleGIS Laboratory team of the University of Cagliari (Italy) for their comprehensive support and collaboration. We also acknowledge the contribution given by the fellow researchers and professors of the Chemical and Geological Department of the University of Cagliari, and particularly Prof. Franco Frau for the support provided for the XRD analyses. We also wish to thank the WADIS-MAR partners for the provided data and support. Among the Tunisian partners, we owe distinct gratitude to the Institut des Regions Arides of Medenine, Tunisia (IRA), in particular to Dr. Mohamed Ouessar, Amar Zerrim, Messaoud Guied, Kamel Nagaz, Mohamed Moussa, Bouajila Essifi, Mongi Ben Zaied, Hanen Jarray, Hanen Dhaou, Dalel Ouerchefani and Amor Jlali for their invaluable support and collaboration during the field campaign, and the undertaken three-months internship.

Disclosure statement

No potential conflict of interest was reported by the authors.

Additional information

Funding

This work was supported by the Ministero dell’Istruzione, dell’Università e della Ricerca; WADIS-MAR Project SWIM programme [ENPI/2011/280-008].