Normalized Sand Index for Identification of Bare Sand Areas in Temperate Climates Using Landsat Images, Application to the South of Romania

Abstract

The expansion of bare sand surfaces indicates a tendency towards desertfication in certain periods as a result of the improper agricultural use of sand soils and of the significant changes in the climate in the past 30 years. The Normalised Sand Index (NSI) is a new index used to identify bare sand areas and their spatio-temporal evolution in SW Romania. Landsat scenes (1988, 2001, 2019), spectral and soil texture analysis (36 samples), covariates (e.g., soil map), and field observations allowed for the validation of the results. The performance of the NSI was compared with indices from the sand index family (e.g., Normalized Differential Sand Areas Index) and supervised classifications (e.g., Maximum Likelihood Classification) based on 47 random control square areas for which the soil texture is known. A statistical analysis of the NSI showed 23.6% (27,310.14 hectares) of bare sands in 1988, followed by an accelerated increase to 47.2% (54,737.73 hectares) in 2001 because of economic and land-use changes, and a lower increase by 2019, which reached 52.5% (60,852.42 hectares) due to reforestation programs. Compared to the NSI, the bare sand areas obtained with the tested indicator were almost 20% higher. The traditional classification shows smaller areas of bare sands but uses a higher complexity of land use classes, while the producer accuracy values are lower than those of the NSI. The new index has achieved a correct spatial delimitation of soils in the interdune-dune and major riverbed-interfluvial areas, but it is limited to the transition Arenosols-Chernozems by humus content and agrotechnical works. The new spectral index favours bare sand monitoring and is a fast and inexpensive method of observing the desertification trend of temperate sandy agroecosystems in the context of climate change.

1. Introduction

The issue of climate change and its effects on natural and man-made environments is of increasing concern to decision-makers around the world. In 2022, several countries of the world adopted the Glasgow Climate Pact (COP26, United Nations Framework Convention on Climate Change) to try to respond to the climate emergency by switching away from fossil fuels [1].

In this context, it should be mentioned that desertification is a topic of utmost importance at the global level; however, traditionally, desertification refers to the extension of desert ecosystems to neighbouring areas due to the reduction of rainfall or due to the displacement of dunes by wind.

In the case of Romania, the mere appearance of sand is interpreted by environmental activists as the result of desertification, but things are much more complex than they seem. In this article, we intend to trace the changes that have taken place over 30 years in the Oltenia area of southern Romania, where there were traditionally shifting sands that had been stabilized by land improvement works during the communist period.

After the fall of the communist system, the changes that occurred in the area had different effects on the local ecosystems.

In Europe, sands of aeolian origin are the most commonly found and widespread, forming the so-called European Sand Belt, while continental sands have multiple sources [2]. The sands of the Oltenia Plain are alluvium transported by the Danube or its tributaries from the Carpathians and the Balkan mountains [3].

Over the last half century, landforms, soils, and land use have undergone several stages of transformation.

In 1970, the development of the Sadova–Corabia irrigation system (>74,000 ha) began, which included the southern part of the study area, and one of the subprojects was to establish modern plantations (e.g., peaches, apricots, cherries, grapes) [4].

Part of the dunes were levelled and excavated for irrigation canals, and the sandy material was redistributed to the interdune areas. The sands were covered with organic material, and the agricultural plots were delimited by forest buffer strips at larger distances for permanent crops and smaller distances for arable land.

The post-communist legislative framework was not very favourable for the maintenance of the irrigation system, plot sizes, and crop types. The irrigation system was an excessive consumer of energy as water was pumped upstream, but during the communist period, little attention was paid to the cost of energy [5]. The administrative transfer of the irrigation system from one entity to another (from national authorities to private associations) due to high energy and maintenance costs, coupled with the lack of subsidies in the following years, led to the deterioration of the irrigation system, especially after 1991 [6]. Following the retrocession of agricultural land [7], the farming areas decreased (5–10 ha), which is also one of the reasons for the abandonment of the irrigation system [8].

The Psamosols of the Oltenia Plain have undergone several stages of deforestation and reforestation. At the end of the 19th century, the sand dunes were stabilized by acacia plantations (about 9000 ha), which were cleared during the development of the irrigation system, thus leading to the reactivation of the dunes [9]. After 1970, however, protective forestry barriers (8–10 m wide, at a distance of 288 m) were set up perpendicular to the prevailing wind direction, some of which were illegally cleared after 1989 [10]. The National Rural Development Programme (2014–2020), funded by the Romanian Government, is a prerequisite for increasing the reforested areas in the Oltenia Plain.

The legislative and economic context has led to changes in the crop structure. Thus, the permanent crops grown on sandy soils (vines, peaches, mulberries) during the period of the irrigation system were replaced by traditional crops (cereals).

The present study aimed at observing how these changes were reflected in the agricultural spatial planning that was shown in the satellite images by using the classical sandy surface mapping indices and proposing a new index.

The problem of desertification in the Oltenia Plain was assessed using remote sensing indices for vegetation (NDVI, MSAVI2) [11] and via the strong correlation between the ground surface temperature obtained from satellite images and the surface temperature from meteorological stations in the area [12]. The intensification of climate aridization, which has been emphasized by future climate scenarios (e.g., RCP4.5), indicates that in the plain and plateau areas outside of the Carpathian Mountains, the desertification process can be significant and long-lasting [13].

Desertification is a dynamic process; therefore, it is necessary to monitor the dynamics at certain time intervals [14]. Spectral indices are used to detect, map, and monitor land degradation. Many studies use sand indices [15,16,17], such as the bare sand index [18,19] or its derivatives [20], which belong to the spectral index family that uses covariates and field surveys to capture the dynamics of sandy areas from medium-resolution satellite imagery.

The mapping and inventory of bare sand provide the baseline data needed in order to prioritize areas and prepare action plans for combating desertification.

The objectives of the study are as follows:

a.

To propose a new index (Normalized Sand Index—NSI) that is able to capture bare sand areas and scattered vegetation covering the sandy soils;

b.

To evaluate of sandy surfaces based on literature indicators and supervised classifications (Maximum Likelihood Classification—MLK, Support Vector Machine—SVM);

c.

To perform a comparative analysis of the accuracy of the NSI with indicators in the literature and with the supervised classifications (MLK, SVM);

d.

To assess the spatio-temporal dynamics of sands (1988–2019) based on satellite images (Landsat sensor).

In remote sensing, most spectral indices can detect dunes and mobile sands in arid and semi-arid climates using spectral images such as visible (VIS), near infrared (NIR), and shortwave infrared (SWIR) by applying normalization [15,17,21,22] (

Table 1. Indices used for mapping bare sand areas.

Index	Bands Math	Feature Extraction, Sand Value/Climate	Satellite	References
Normalized Differential Sand Dune Index (NDSDI)	$NDSDI=\frac{R-SWIR2}{R+SWIR2}$	Sand, ˂0/dry	Landsat 5 TM, Landsat 7 ETM	[15]
Normalized Differential Sand AreasIndex (NDSAI)	$NDSAI=\frac{SWIR1-R}{SWIR1+R}$	Sand, ˂0/dry or humid	Landsat 5 TM, Landsat 7 ETM, Landsat 8 OLI	[17]
Normalized Difference Enhanced Sand Index (NDESI)	$NDESI=\left(\frac{b4-b2}{b4+b2}\right)+\left(\frac{b7-b6}{b7+b6}\right)$	−2, +2/arid	Sentinel 2 and Landsat 8 OLI	[21]
Sand differential emissivity index (SDEI)	$SDEI=\frac{b13-b12}{b13+b12}$	1 to 0.28/extremely arid	Aster	[22]

). Other indices use band 1 (Coastal aerosol) and partial normalization [19] or employ a surface reflectivity analysis using red and near-infrared bands and correlation with fractional vegetation cover [23,24]. Different from single index modelling, desertification is quantified by combining an active surface property (e.g., albedo) with vegetation and soil indices (e.g., soil-adjusted vegetation index) in point-to-point and point-to-line models [25].

The indices for sand use simple mathematical equations [21] that are easy to apply, and if the validation is supported by the analysis of properties with low variability over time (e.g., soil texture versus humus content), the results are good. There is no best method to map bare sands. In many papers, data normalization is used [15,17,26], while in others, partial normalization is preferred [19,27] as a way to differentiate sand from soils. The disadvantage of indices that detect sand is their exclusive use in arid climates [21,26], thus making the results obtained in temperate climates not adequate [27].

2. Study Area

The study area is bounded to the west and south by the slope that transitions to the floodplains of the Jiu and Danube rivers (Figure 1). The eastern boundary is characterized by a change in soil type and texture, i.e., the transition from Arenosols to Chernozems. Towards the east, a buffer of 2000 m was generated based on the soil sampling points so that the sandy boundary could be correctly identified.

Sand deposits in the study area are located on the left side of the Jiu River, taking the form of an elongated strip in the north–south direction (about 70 km) and with a narrower width in the north (4–6 km), but reaching 28–30 km in the southern sector (Dăbuleni). These are closely related to the paleogeographic evolution of the Jiu and Danube rivers. In the Pleistocene and Holocene, the course of the Danube moved southwards, leaving a well-developed system of up to six river terraces on the left side. The terrace deposits consist of alternating sand, clay, and loess deposits.

Morphologically, two types of dunes can be distinguished. There are old, stabilized, or semi-stabilized dunes that have a symmetrical profile. The interval between two rows of dunes varies between 100 and 500 m, and their height does not exceed 8–15 m. Towards the north, the gap between the dunes increases with distance from the Danube while the height decreases. There are recent, mobile, or semi-stabilized dunes that have an asymmetrical morphology, with a less steep slope towards the wind and steeper slopes on the opposite side. They may occur in the form of barrens that are very sparsely covered by vegetation, and they consist of quartz sands and additional minerals such as mica [3,28].

The climate of the Oltenia Plain is temperate-continental, but spring time is the most critical period of deflation because the soil has not yet been covered by vegetation, and the most severe droughts affect agricultural crops [29], therefore increasing the sandy areas that are not covered with vegetation for a longer period of the year.

The dominant wind direction is NW-SE, with an average speed ranging from 3.2 m·s⁻¹ from the NW to 3.8 m·s⁻¹ from the W, but with speeds that can exceed 15–17 m·s⁻¹ becoming more frequent in recent years [30]. The most active and harmful form of wind erosion occurs in the steppe and partly wooded-steppe zones, where wind speeds exceed 4 m·s⁻¹ [31].

Soil distribution and texture were taken from the Map of Romanian Soils (1:200,000), and the names were adapted according to the International System for Soil Classification (WRB) [32].

Soils are represented by Arenosols, which vary spatially according to the configuration of the land. Sandy deposits are favourable for the formation of Arenosols [33], but mollic arenosols occupy the interdune areas, while Eutric Arenosols are dominated on dune tops [34].

In the sector related to the irrigation system (Figure 2), the sandy soils have been anthropogenically modified by levelling, resulting in Arenosols being displaced and transported. A particular characteristic of the sandy soils of the Oltenia Plain is the very high proportion of coarse sand (60–80%) and the very low proportion of clay (2–19%) [35].

Towards the east, the sand deposits thin out, intercalating into Cambic Chernozems, while in the northeast, the transition is better expressed towards Chromic Luvisols (Figure 2).

The texture on the surface horizon is sandy in the west and changes to sandy loam towards the east and northeast. In the south-eastern half, the sandy-loam texture sometimes forms bands to the south of the sandy-textured soils (Arenosols) as a consequence of the dominant wind direction (VNV-ESE) (Figure 1).

Land use for the first hierarchical level in Corine Land Cover (2018) [19] includes: agricultural land, forests and semi-natural areas, water bodies, and artificial surfaces. Agricultural land nomenclature has been adapted to local characteristics, and the following classes have been identified: arable land and pastures on sandy soils (APSS), arable land-pastures on soils with different texture (APSDT), autumn crops (AC), and permanent crops (PC). Forests and semi-natural areas include: compact forests (CF), scattered forests (SF), and sands (S).

Arable land on sandy soils (Arenosols) is partially covered by vegetation, and the plots are rectangular in shape. Arable land and grassland on soils with other textures have low reflectance values due to humus content on the surface horizon, and the plot size is larger than in the previous case. Autumn crops are found on both soil types. The moisture of the Arenosols influences plant density so that the identification of sandy textured soils can be blurred in rainy years. Permanent crops (PC) are represented by vineyards and orchards, which are traditional crops on sandy soils.

Compact forests (CF) are represented by even-aged acacia plantations, and scattered forests (SF) are older acacia plantations where sparse trees form mixtures with grassy vegetation and scrub. Detecting and analysing the spatial dynamics of scattered forests can be an indicator of areas at risk of waste. Bare sand areas (S) are frequently encountered in the SE of the territory, within modified relief. Water bodies include two categories: small lakes (L) resulting from the damming of streams and lakes of semi-permanent character, located between dunes.

3. Materials and Methods

The study methodology involves both the use of remote sensing data such as satellite images and the use of existing GIS databases such as soil maps, topographic maps, Corine Land Cover, etc., in addition to field data processed in the laboratory. A brief methodological workflow is presented in Figure 3.

Three methods were used for mapping areas with uncovered or poorly covered sand: (a) the introduction of a new index, called the Normalized Sands Index (NSI), (b) the use of spectral indices from the literature that are able to map sandy and bare soils (Table 1), and (c) the use of a traditional pixel-based image classification method (e.g., Maximum Likelihood image) to facilitate the interpretation of sand bars in relation to land use. Finally, the performance of NSI is validated by statistical analysis.

3.1. Data Types

Due to the fact that the study of land use evolution on sandy soils in Oltenia spans a period of 30 years, we decided to use Landsat satellite images as these are characterized by a high temporal consistency; images from the communist period (1988), the first agricultural developments, and recent images (2019) are available with comparable spatial and radiometric characteristics. In this study, we used the highest quality data available for Landsat scenes (Tier 1), which are suitable data for time-series processing analysis [36] and are consistently georegistered within a ≤12 m radial root mean square error, making them suitable for time-series pixel-level analysis [37]. The study area is entirely included in path 184, row 29 for all Landsat scenes, and therefore spatial enhancement is not necessary.

Multispectral Landsat TM5, Landsat 7 ETM, and Landsat 8 OLI data were used for the remote sensing of surface soil attributes, and each row/path date is shown in

Table 2. Landsat scenes used in the analysis.

Satellite	Row	Path	Date of Acquisition
Landsat 5 TM	29	184	27 January 1988
Landsat 7 ETM	29	184	7 February 2001
Landsat 8 OLI	29	184	17 February 2019

. The scenes were obtained from the United States Geological Survey (USGS) and downloaded from the USGS Earth Explorer online platform; all scenes are included in the Level 1 precision terrain (L1TP).

The selection of spectral data was based on several factors that take into account the characteristics of the objects on the images and the variation of the climate during the year. Firstly, it was ensured that the entire study area was free of cloud, haze, and fog. The selection of scenes and time intervals took into account the economic changes outlined in the introduction and the criteria for processing Landsat images. The selection of spectral bands in sandy soil areas for the multi-temporal analysis of land degradation dynamics took into account annual climatic conditions in semi-arid climates, with images from both high rainfall months and months of maximum vegetation growth [38] and low rainfall [39]. Unlike in the drylands, in the temperate zone, the use of images from the same season and the successive years was not possible due to cloud cover.

The choice of images from the winter season took into account the following aspects:

a.

Identification of the maximum uncovered sand areas during the year;

b.

Avoiding confusion between the light (yellow) colour of the sand and the straw resulting from grain harvesting;

c.

In autumn, plant residues are burnt and dark coloured surfaces can be confused with the soft horizon of Chernozems, while smoke influences the quality of the images;

d.

In winter, evapotranspiration is lower than in summer, therefore image quality will be less affected by the volume of vapour in the atmosphere. Winter evapotranspiration dynamics, analysed over a long period of time, was more stable compared to other seasons [40].

Soil reflectance is affected by the presence of vegetation, water bodies, irrigation, and other land features. In this study, an attempt was made to obtain the highest possible image quality (with the maximum number of pixels) of the bare soil in the research area. To get the least vegetation cover, the satellite images were acquired in winter (January–February), when the soil was not covered by snow.

The field study looked at two aspects. A total of 11 pairs of soil samples (1a,b to 11a,b) (Figure 1) were taken from the eastern extremity of sandy soils (a) to show the maximum distribution of Arenosols and the transition mode to other soils, and (b) to identify land use types in the vicinity of the sampling points that validate the Landsat 8 OLI (2019) scene classification and NSI accuracy.

Field observations have allowed for a differentiation of the soils, based on the texture of the surface horizon, into two classes: sandy (sandy and sandy loam) and soils with a different texture (loam). For each point (Figure 1), soil texture was estimated by testing the sandy texture between fingers, according to the FAO flowchart [41], with the relative deviation from laboratory analysis being lower (3.8%) for sand [42]. The coordinates of the sampling points were acquired with a RTK UNISTRONG G975 geodetic GPS.

Soil samples weighing 500 g were collected from the surface soil horizon. After drying in the laboratory for 14 days, the colour of each soil sample was determined with the Munsell Soil Colour Chart (2000) [43]. After sample preparation [44], soil texture was determined using Analysette 22NanoTec (Fritsch GmbH, Idar-Oberstein, Germany). Then, the sand (2–0.02 mm), silt (0.02–0.002 mm), and clay (˂0.002 mm) were classified in accordance with the Romanian soil Taxonomy [45], which is compatible with the texture used by the International Society of Soil Science. This approach allowed for a comparison of the texture results of this study with those in the literature. In order to validate the NSI results, the weight of sand on the surface horizon from 14 soil profiles (12 to 25) specified in the literature was used (Figure 1) [34,35].

3.2. Data Processing

The Landsat scenes were pre-processed (radiometric calibration, atmospheric correction) and subsequently spectrally analysed and processed.

The Landsat scenes were radiometrically corrected [46], thus obtaining surface reflectance, and a Quick atmospheric correction was then applied; all images were reprojected in Stereographic 1970 (EPSG 31700) using the ENVI (Exelis Visual Information Solutions) software.

In many studies, soil properties were investigated based on the visible (VIS), near infrared (NIR), and shortwave infrared (SWIR) spectral range [47]. Landsat bands are known for particular applications: band 3 (discriminant of soil from vegetation), band 5 (soil–rocks), band 7 (geology) [48]; the red band (Landsat 5 TM) was used to discriminate bare soil from other land [9].

The multispectral digital classification of remote sensing data requires the selection of only the bands containing the most information. The high value of the Optimal Index Factor (OIF) indicates that the bands contain a lot of information (e.g., high standard deviation) with little duplication (e.g., low correlation between the bands) [49].

ILWIS 3.8 Academic software and the Optimal Index Factor (OIF) algorithm were used to select the bands. For each scene, two combination possibilities were identified based on the highest Index Highest Ranking (IHR) values: bands 4,5,7 (1527.26) and bands 3–5 (1282.48) for Landsat 5 TM; bands 4,5,7 (1597.26) and bands 3,5,7 (1406.53) for Landsat 7 ETM; bands 5–7 (2982.32) and bands 4,5,7 (2982.32) with bands 3,5,7 (2449.33) for Landsat 8 OLI.

Prior to soil sampling, the eastern boundary of the sandy textured soils for scene Landsat 5 TM (27 January 1988) was estimated by transforming the synthetic colour image (Transformation module in Envi 5.3), which enabled us to locate the soil sampling points, noted on the map from the 1st to the 11th (Figure 1). The pair of each point was identified in the field based on texture and noted from 1b to 11b (Figure 2). Using ArcGIS 10.2 (Environmental Systems Research Institute, Redlands, CA, USA), the twenty-two points were generated in the eastern part of the sand area, which were used for sampling and land-use validation. In the second step, the position of the points was maintained or modified according to the texture estimation in the field.

Using ArcGIS 10.2, the representative area for each land use type was digitized, resulting in a different number of ROIs. When the landscape of a study area is complex-heterogeneous, selecting sufficient samples becomes difficult [50]. Auxiliary data were used to extract training areas associated with each land cover type. Thus, certain land use types were extracted from the topographic map (1984) and used for the classification of the Landsat 5 TM scene (1988), and the CLC (2000) [51] was used as a covariate to identify certain land use types (e.g., vineyards and orchards) for scene Landsat 7 ETM (2001). Only a few land cover types could be identified on the basis of topographic map information (lakes, vineyards and orchards, forests), while others could be inferred from plot shape (e.g., agricultural land on sandy soils is rectangular in shape and separated by tree belts).

Within localities, different pixels’ reflectance values are determined by the diversity of material types and agricultural practices (e.g., addition of natural fertilizers, burning of vegetation). Under these conditions, the small areas within the localities induce a high level of noise in the identification of some classes. In the study area, there were only rural settlements whose area had not changed considerably during the study period. Based on the above considerations, the perimeter of the localities was extracted from Landsat 8 OLI (2019), which represents the unique mask for the whole set of analysed images.

The Jeffries–Matusita (JM) separability criterion is a common statistical index that allows for the selection of a suitable feature subset associated with high class separability [52]. A spectral curve for each training area associated with one land use was directly determined by the software ENVI 5.3. For increased sand separability, we used a set of non-overlapping areas for each Landsat scene. This allowed us to retrieve the average spectral patterns of bare sand and compare them with other land-use types (Figure 4). In such a comparison, the spectral curve could be used to indicate the spatial cover of sands with the other land components (e.g., lakes). In the case of sandy soils, some studies show that spectral reflectance in Landsat bands increase rapidly to B4, with a peak in B5, especially in the absence of opaque minerals [19,53].

Based on the JM-transformed divergence algorithm, which was applied to eight classes of land use, twenty-eight pairs with different values were obtained. The separability of the selected training sites for all classes was examined by computing their spectral separability in the ENVI 5.3. software.

3.3. Normalized Sand Index (NSI)

To date, the most commonly used indices for assessing sand surfaces, based on pixel information from satellite images, are the bare sand indices (

Table 3. Indices used to extract sandy surfaces.

Index	Band Math	Feature Extraction, Sand Value/Climate	References
Normalized Differential Sand Dune Index (NDSDI)	$NDSDI=\frac{R-SWIR2}{R+SWIR2}$	Sand, ˂0/dry	[15]
Normalized Differential Sand Areas Index (NDSAI)	$NDSAI=\frac{SWIR1-R}{SWIR1+R}$	Sand, ˂0/dry or humid	[17]

).

Two traditional supervised classifications were tested: Maximum Likelihood Classification (MLK) and Support Vector Machine (SVM); however, MLK was used as it is one of the most popular and widely used conventional classifications, thanks to its robustness [54] and its high performance for regional scale land cover mapping [55].

In order to get better results for our geographical area, we proposed a new index—Normalized Sand Index (NSI), which is calculated based on the Red (R), Green (G), and Short-wave Infrared (SWIR1) bands (Equation (1)). In order to minimize the differences induced by atmospheric and moisture conditions between images collected in different time periods, the sand index values were max–min normalized, resulting in the linear transformation of the data in the range 0–1, which preserves the relationships among the original data values. The application of normalization is a widely used method that restricts a large range of values to a small range (e.g., 0–1). Under these conditions, we can compare images from different years or results from different indices, applied to the same type of surface (e.g., sand) [56,57,58].

NSI = (G + R)/(log (SWIR1))

(1)

NSI is an index that can be used in the temperate zone, where bare sands and soils have different moisture values. For example, in a major riverbed, soil moisture is higher than on the slopes, and the soil is covered with vegetation. NSI applied to satellite imagery from 2019 maps these details better than NDSAI. Under these conditions, NSI is able to map sandy surfaces without being influenced by moisture. NSI does not generate confusion between light-coloured concrete irrigation channels or other surfaces (e.g., roads) and sand bars.

In remote sensing, some software (e.g., ENVI 5.3) uses traditional transformations (Principal Component Analyses, Tassaled Cap) in order to highlight certain properties of land surface objects, but also the dynamics of some phenomena (e.g., desertification) [59,60].

The mapping of some soil properties (e.g., salinization, sodicity) can be improved by transforming variables (1/log), and the spectral reflectance increases with the value of electro-conductivity (e.g., EC) [61]. Bare sand as well as salinized soils show high reflectance values. The distribution of DN values for SWIR 1 (2001) analysed in SPSS 10 does not clearly show the presence of sandy textured soils. Following the SWIR1 log transformation, an increase in reflectance values was observed for sandy textured soils, with a clear clustering in the range of 7.5–8.5. Cartographically, sandy soils are better highlighted and will have higher values than soils with other textures.

The classification accuracy check followed two components:

(a)

Assessment of the accuracy of the NSI classification against the selected spectral indicator (NDSAI);

(b)

Assessment of the ability of traditional classification (TL) to map bare sand surfaces.

The priority was to identify the threshold at which the NSI value and the tested indicators show the presence of sand. For this purpose, 47 polygons of 3 × 3 pixel sides were constructed in Arc GIS 10.8; these will be referred to as control square areas hereafter. By observing each control square area, it was found that it may also contain pixels indicating sand with sparse vegetation cover. Using only pure pixels to assess accuracy would result in unduly high values for overall accuracy [62]; therefore, decreasing the size of the control square areas would not be justified.

Since NSI tests for the existence of uncovered sands, each control square area was assigned to sandy textured soils (Arenosols) based on the soil map (Figure 1). For some of the 1–3 control square areas, there was the uncertainty of not belonging to Arenosols as they were placed at the transition to soils with other textures.

The minimum threshold value at which pixels were assigned to sand, as estimated by the NSI index based on satellite images, was extracted by the statistical calculation of the 9 pixels integrated in the averaging. The threshold value was identified from the histogram of the 47 control square areas by summing the mean ± standard deviation [63]. The minimum threshold value for NSI varied from year to year as follows: 0.38 in 1988, 0.24 in 2001, and 0.19 in 2019.

3.4. Accuracy of Traditional Classifications

The confusion matrix from the ENVI 5.3. software was the basis for the evaluation of the MLK classification performance, and the correlation matrix for the NSI and NDSAI indicators was generated in Arc GIS 10.8. Due to the large time difference between the 1988 and 2001 scenes, a confusion matrix using ground truth image was preferred. For the scene of 17 February 2019, ground truth ROIs were generated based on observations in the vicinity of the sampling points (1a,b–11a,b) and texture on the surface horizon [34,44].

The classification accuracy for multi-temporal studies and soil properties was evaluated based on overall accuracy and the Kappa coefficient [64], or the user’s and producer’s accuracy per class [65].

4. Results and Discussion

The assessment of bare sands or soils with high sand content is really difficult in temperate areas because the precipitation excess—compared to the average value in desert areas—creates favourable conditions for the growth of grassy vegetation. Natural and cultivated vegetation is closely related to the consumption habits of the community.

The different configurations of the spectral curves in SWIR 1 and SWIR 2 in 2019 as compared to 1988 and 2001 may be a consequence of the change in land use and thus the number of pixels allocated to each training area (

Table 4. Pixels count for each training areas.

Land Use	Pixels Count for Each Land Use
	27 January 1988	7 February 2001	17 February 2019
Autumn crops (AC)	7702	7730	4274
Arable lands-pastures on sandy soils (APSS)	3724	3213	2167
Arable lands-pastures on different soil textures (APSDT)	8651	7688	6104
Permanent crops (PC)	2995	1298	778
Compact forests (CF)	933	809	1619
Scattered forests (SF)	2123	4013	2140
Sands (S)	464	627	635
Lakes (L)	120	112	146

).

For the distance, Jeffries–Matusita (JM) values higher than 1.9 were interpreted as very good, and separability values between 1.5 and 1.9 indicated good separability [66].

Low separability values may have been the result of the combined types of vegetation (such as cultivated, herbaceous, bushes) and soil surface characteristics (e.g., colour). In the case of low separability, two classes indicating the same object may be merged into one. The weak values of the JM results obtained for the CF-SF pair are not a major problem because the two forest types can form a single class, but this segregation was preferred for observing the spatial dynamics of forest plantations (Figure 5).

The band selection was based on the spectral behaviour of sand surfaces in the range of visible bands (Green and Red) and for other land uses in the Short-wave Infrared band (SWIR 1). In the studied area, the sand showed a continuous increase in the reflectance of the Green (b2, Landsat 5 TM, Landsat 7 ETM 0.53–0.59 µm, b3 Landsat 8 OLI = 0.53–0.59 µm) and Red (b3, Landsat 5 TM, Landsat 7 ETM 0.63–0.69 µm, b4 Landsat 8 OLI = 0.64–0.67 µm) bands, but with a variation for SWIR1 (Figure 4). The reflectance of the SWIR1 bands (b5, Landsat 5 TM, Landsat 7 ETM 1.55–1.75 µm, b6 Landsat 8 OLI = 1.57–1.65 µm) accentuated the difference between soil moisture content and vegetation.

From the category of indices used for bare sand and sand dunes, those that generated normalized values were preferred (Table 4), such that the results are comparable with those obtained by the indicator proposed in this study (NSI). Two indicators were tested to identify sandy surfaces: NDSAI [17] and NDSDI [16]. However, NDSAI was selected because its spatial distribution was closer to NSI, and because it could capture sandy areas in wet climates [17]. The dynamics of non-vegetated sandy areas are associated with the desertification process. One of the methods used to monitor bare sands is the use of spectral indices (e.g., Bare soil index, Normalized Differential Sand Dune Index, Normalized Differential Sand Areas) [15,17,20]. Some of these monitor sandy surfaces in arid climates, while the Normalized Differential Sand Areas (NDSAI) uses Red and SWIR1 bins that are sensitive to soil moisture [17]. As the study area is located in the temperate zone, the snow sands received variable amounts of water. Of the tested indicators, NDSAI best estimated the areas with sand that lacked vegetation, including wet areas where sand rapidly loses moisture, and was compared with the indicator proposed in this study (NSI).

Compared to the uncovered sandy areas identified by NSI and shown above, NDSAI captures larger areas at about 20% of the total number of pixels. Thus, in 1988, the uncovered sand area was 45.5% (586,408 pixels), it reached 57.4% (739,699 pixels) in 2001, and it continued to increase to 71.6% (922,397 pixels) in 2019. Climatically, February 2019 was considered the warmest month, with an average of 2.7 °C, and it was excessively rainy; the last decade, however, has seen cooling air masses and wind intensities reaching 80 km h⁻¹ [67]. Thus, as the wind velocity increases, plant roots are exposed, and plants are damaged by grain bombardment or by being covered with sand [68]. These rapid changes in weather can cause the unvegetated sand surfaces to manifest in a pulsating manner from year to year, and even within the same season.

The comparative analysis of the user/producer accuracy matrix shows the close values of the two indicators, but NSI stands out with higher values for producer accuracy in Arenosols only in 2001 (

Table 5. NSI and NDSAI calibration accuracy analysis for Arenosols.

Index	NSI 1988				NDSAI 1988
	NSI (1)	NSI (2)	Total (User)	User Accuracy (%)	NDSAI (1)	NDSAI (2)	Total (User)	User Accuracy (%)
Arenosol (1)	42	8	50	84	44	9	53	83
Other soil classes (2)	3	5	8	62.5	1	4	5	80
Total (Producer)	45	13	58	0	45	13	58	0
Producer accuracy (%)	93.3	38.5	0	81.0	97.8	30.8	0	82.8
Overall accuracy for arenosols (%)				81.4				82
	NSI 2001				NDSAI 2001
Arenosol (1)	44	8	52	84.6	43	9	52	82.7
Other soil classes (2)	1	5	6	83.3	2	4	6	66.7
Total (Producer)	45	13	58	0	45	13	58	0
Producer accuracy (%)	97.8	38.5	0	84.5	95.6	30.8	0	81
Overall accuracy for arenosols (%)				84.5				97
	NSI 2019				NDSAI 2019
Arenosol (1)	43	8	51	84.3	44	10	54	81.5
Other soil classes (2)	2	5	7	74.4	1	3	4	75
Total (Producer)	45	13	58	0	45	13	58	0
Producer accuracy (%)	95.6	38.5	0	85.8	97.8	23.1	0	81
Overall accuracy for arenosols (%)				82.4				81

). For the other soil types (e.g., Cernozems), which also have a sandy texture, the user accuracy results obtained with NSI are better, as a spatially correct delimitation of soils in the interdune-dune and major riverbed-interfluvial areas was obtained.

The overall accuracy values obtained with NDSAI for the study territory are comparable to that determined by Sahar et al. [17], but NSI has a lower deviation from the three-year average (82.7%) as compared to NDSAI. If NDSAI is the best sandy index for mapping most of the sandy surfaces [17], NSI highlights the bare sand surfaces in a temperate zone territory characterized by a higher complexity of soil types and cultivated species, but also by anthropogenic influence on the environment.

The statistical processing of the NSI data shows that the area of bare sand continuously increased at an accelerated rate from 1988 (23.6% representing 303,446 pixels) to 2001 (47.2%, representing 608,197 pixels), and then moderately until 2019 (52.5%, representing 676,138 pixels).

The higher frequency of dry years in the 1983–2002 period [44] may be one of the causes of the expansion of uncovered sand areas. This is compounded by intervals of extreme temperatures, such as those in June and July 2000, which exceeded 43 °C in the southern part of the territory and had heavy rainfall recorded on a single day [69]; these increased the area of uncovered sand due to erosion. It is possible that the NSI performance shown by producer accuracy was also influenced by the climatic characteristics of the year 2000.

The classification maps produced from the implementation of the MLK classifier is illustrated in Figure 6, and the accuracy statistics are summarized in

Table 6. User’s and producer’s accuracy by classes.

MLK 1988					MLK 2001				MLK 2019
Land Use	Prod. Acc.	User Acc.	Prod. Acc.	User Acc.	Pro. Acc.	User Acc.	Prod. Acc.	User Acc.	Prod. Acc.	User Acc.	Prod. Acc.	User Acc.
	(%)	(%)	Pixels	Pixels	(%)	(%)	Pixels	Pixels	(%)	(%)	Pixels	Pixels
AC	98.87	95.64	263/266	263/275	93.24	92.34	193/207	193/209	100	99.38	641/641	641/645
APSS	95.68	86.46	332/347	332/384	92.27	85.42	334/362	334/391	62.99	93.27	291/462	291/312
APSDT	92.35	88.95	169/183	169/190	94.01	90.23	157/167	157/174	100	87.77	165/165	165/188
PC	79.82	77.39	178/223	178/230	81.72	80	228/279	228/285	62.07	31.03	36/58	36/116
CF	80.88	85.94	55/68	55/64	75.7	76.42	81/107	81/106	79.29	91.79	425/536	425/463
SF	58.08	80.83	97/167	97/120	61.54	79.28	88/143	88/111	85.96	57.83	251/292	251/434
S	66.67	94.74	18/27	18/19	52.17	92.31	12/23	12/13	92.65	94.19	227/245	227/241
L	75	100	3/4	3/3	0	0	0/1	0/0	100	100	56/56	56/56

.

Although climatic conditions are not very favourable for revegetation during dry years, which have also been very frequent after 2001 [30,70], natural revegetation through the abandonment of agricultural land [71] and national revegetation programs slowed the pace of sand expansion until 2019.

Overall accuracy (Oa) and the kappa coefficient (k) obtained for MLK in the three years were close (1988-L5, Oa = 86.77%, k = 0.83%; 2001-L7, Oa = 84.79%, k = 0.81%; 2019-L8, Oa = 85.21%, k = 0.82), but it is recommended that the kappa coefficient be discarded, which is not an adequate index for describing classification accuracy [63]; alternatively, user’s and producer’s accuracy per classes was used.

For the sand class, user accuracy was considerably higher for the three years (Table 6), and this could be explained by the constant shape of the spectral signature in the green and red bands. The homogeneous spectral signature associated with high building density surfaces resulted in higher values for user accuracy than producer accuracy [34], and sand indicated higher homogeneity than the other land use classes in the studied case.

Some Chernozems do not meet the diagnostic criterion for colour to be classified as a mollic horizon (e.g., samples 4a,4b,5a) (

Table 7. Soil properties and land use.

Soil Sample Code	Latitude	Longitude	Fine Sand a	Coarse Sand b	Sand	Land Use/Vegetation	Soils	Soil Colour c
1a	44°19′91″	23°91′59″	24.9	62.1	87	Forest	Arenosols	10YR8/8
1b	44°20′16″	23°91′82″	71.2	22.1	93.3	Arable		10YR6/4
2a	44°12′31″	23°93′63″	49,8	36.2	93	Herbaceous plants	Arenosols	10YR6/4
2b	44°12′16″	23°93′17″	12.6	87.1	86	Herbaceous plants		10YR7/6
3a	44°86′72″	23°59′34″	29.8	59.4	89.2	Herbaceous plants (grassland)	Luvisols	10YR6/3
3b	44°83′93″	23°59′39″	49.5	41	90.5	Herbaceous plants (grassland)		10YR6/3
4a	44°24′64″	24°71′31″	19.2	70.5	89.7	Arable	Chernozems	10YR6/3 c
4b	44°24′10″	24°65′21″	13.8	79.6	93.4	Arable		10YR6/4 c
5a	44°01′30″	24°10′49″	17.3	75.6	92.9	Arable	Chernozems	10YR6/3 c
5b	44°01′38″	24°10′49″	14.9	79	93.9	Arable		10YR5/3
6a	43°59′14″	24°12″	17.6	75.8	93.4	Arable	Chernozems	10YR5/3
6b	43°98′41″	24°20′91″	8.3	87	95.3	Arable		10YR4/3
7a	43°55′49″	24°15′21″	14.8	77.7	92.5	Vineyard	Chernozems	10YR4/2
7b	43°92′85″	24°26′41″	19.9	72.4	92.3	Arable		10YR5/3
8a	43°84′46″	24°25′56″	20.1	72.9	93.0	Arable	Arenosols	10YR4/3
8b	43°84′49″	24°25′66″	19.8	71.9	91.7	Arable		10YR5/4
9a	43°79′66″	24°2′46″	52.4	33.4	85.8	Herbaceous plants	Arenosols	10YR5/2
9b	43°79′64″	24°24′78″	25.4	62.9	88.3	Arable	Cherozems	10YR4/2
10a	43°76′56″	24°25′77″	0.7	85	85.7	Arable	Arenosols	10YR4/4
10b	43°76′34″	24°26′4″	57.6	38	95.6	Arable	Cherozems	10YR4/3
11a	43°75′95″	24°23′46″	23.2	72.4	95.6	Abandoned vineyard	Arenosols	10YR4/4
11b	43°75′88″	24°23′79″	25.1	62.4	87.5	Abandoned vineyard		10YR4/3

Fine sand ^a = 0.02–0.2 mm, coarse sand ^b = 0.2–2 mm, 10YR6/3 ^c = colour criteria not met for mollic horizon.

), being less than five within dry material [32]. This feature suggests the idea of the sands extending eastward from the active areas (e.g., east of Stefan cel Mare locality).

There are no significant differences between the texture values in the sampling points (Table 7) and those on the surface horizon of the soils in the literature (average sand in samples = 91.5%, in profiles = 89.14%; StdDev in samples = 3.52, in profiles = 4.35) (Figure 7). A detailed study of the soils of the Oltenia Plain revealed the existence of several soil types (e.g., Cambisols) included in the Arenosols area [35], but also some particular properties of the latter [34]. The maximum and minimum texture values express both the pedological diversity and anthropic influence on soils (e.g., stripping/covering).

The statistical analysis of land use classes obtained from MLK did not reflect the increase in uncovered sand areas in 2019 as compared to 2001, and APSS was observed to be declining. The lower exposure of soils to deflation through the decrease in the area of autumn crops and the areal extension of compact forests in the period of 2001–2019 explains the slower evolution of uncovered sandy areas (

Table 8. Pixels count and percent for each land use (according to MLK).

Land Use	27 January 1988		7 February 2001		17 February 2019
	MLK		MLK		MLK
	Pixel Count	%	Pixel Count	%	Pixel Count	%
AC	266,683	8.14	207,204	6.33	132,014	4.02
APSS	346,531	10.58	362,030	11.06	265,433	8.1
APSDT	184,225	5.62	168,077	5.13	265,382	8.1
PC	224,362	6.85	278,978	8.52	208,548	6.36
CF	68,577	2.09	106,872	3.26	141,542	4.32
SF	167,130	5.10	142,871	4.36	256,403	7.82
S	26,718	0.81	22,530	0.68	20,179	0.61
L	3920	0.11	1278	0.03	1214	0.03

). In addition to these, there was also the abandonment of some agricultural land due to the low productivity of sandy soils (e.g., Sadova) [71]. The values for S and APSS resulting from the MLK classification (Table 8) indirectly validated the NSI, which reflects a moderate increase in the last interval.

NSI emphasizes certain peculiarities in the distribution of unvegetated sands according to land use. The rectangular shape of the unvegetated areas in the north-west corresponds to large agricultural land from the period of planned agriculture (1988) (Figure 6a). The high NSI values from the same area in 2001 are due to excessive land parcelling (Figure 6b) and poor soil cover by autumn crops (Figure 6h), which are factors favouring the reactivation of sands. In response to land degradation, the area of compact forests was extended by acacia plantations in the north-west of the territory (e.g., near Marsani) [72] (Figure 6c,i).

After 1989, the surfaces occupied by vineyards or by orchards decreased abruptly, as opposed to the surfaces used as pastures and agricultural land; in the latter case, an increase could be observed [73,74,75]. The deforestation of vine plots to make way for new plantations damaged the fragile equilibrium of sandy soils. In the Oltenia region during the period of 2007–2018, through the conversion/restructuring programs, there was a significant increase in the areas planted with vines destined to produce quality wines [76]. The disappearance of some traditional crops on Arenosols (Figure 8a) in the southern part of the territory is evidenced by NSI in the rectangular shape of the unvegetated plots, (Figure 8d) but with smaller dimensions as compared to the agricultural lands in the north-west.

With the levelling of the relief, the lakes between the dunes (Figure 8a) disappeared, leading to the expansion of the uncovered sandy areas (Figure 8b). The mitigation of hot areas was achieved by the expansion of compact forests (Figure 8d, Table 8), with the rate of reforestation in SW Oltenia being higher than the rate of deforestation [77], at least in some periods.

In the central and southern part, linear dunes are more clearly expressed in the relief than in the north, their characteristics being parallelism, regular spacing, and partial vegetation [78]. The dune peaks, in the shape of an arrow, intersect with the rectangular surface of the farmland, penetrating the Chernozems area (Figure 9a). The change in colour of the A horizon over a distance of several meters marks the transition from Arenosols to Chernozems (Table 7). The agricultural use of dune soils keeps the peak active, but agricultural work and wildlife in Chernozems limit the correct estimate of peak dynamics (Figure 9d).

In the areas where the dunes were modelled, the unaltered sand brought to the surface had high reflectivity, as was the NSI value for 1988 (Figure 10a). Over time, the alteration and accumulation of organic matter attenuated the NSI value (2001), and through land cover, bare sands became less identifiable (Figure 10c).

Most often, the boundary between two different soil types is gradual, but while one property (texture) may be similar for both types, another differentiates them significantly. For example, Arenosols and the other soils (Chernozems and Luvisols) have a high proportion of sand, while the humus content is low in Arenosols (<1–2%) and high for Chernozems (>5%) [34].

Soil is a dynamic system, and profile development through the aeolian input of material is a component of progressive pedogenesis in which current processes and land use are able to assimilate the addition of material from the surface. The aeolian input of material can be assimilated with regressive pedogenesis if the material prevents horizon differentiation or profile development at depth [79]. The transport and accumulation of sand on the surface horizon of Chernozems (CZ) near Arenosols influences the spectral reflectance and therefore the NSI value at different time intervals (Figure 10d).

5. Conclusions

In this study, the evolution of bare sand surfaces was assessed using a new Normalized Sand Index (NSI) indicator based on a Landsat medium-resolution image series from 1988, 2001, and 2019.

The existence of bare sandy areas is conditioned by anthropogenic factors, in particular land use and the variability of climatic conditions. The development of the irrigation system in the south-eastern part of the territory and the large plots in the north, which lack erosion protection from the period of planned agriculture before 1989, have kept the bare sand areas active. After this date, the decrease in irrigated areas, the abandonment of traditional crops on sandy soil, and the decline in forested areas led to an increase in bare sand surfaces. Sand areas increased by 2019, but at a slower pace due to national reforestation programmes and the natural revegetation of abandoned agricultural land.

Statistical analysis of the NSI showed the presence of bare sands at 23.6% (27,310.14 hectares) in 1988, followed by an accelerated increase to 47.2% (54,737.73 hectares) in 2001, and a smaller increase to 52.5% (60,852.42 hectares) by 2019.

Compared to NSI, the bare sand areas detected with the tested indicator were almost 20% higher. The traditional classification (MLK) showed smaller areas of bare sand but used a higher complexity of land use classes; the producer accuracy values were also lower.

Bare sand surfaces in temperate climates can increase or decrease from one year to the next, depending on climate conditions and even during the same season. The uncovered sandy area in 2001 may have been a consequence of the previous year′s drought, but the use of images from the succeeding years was limited by cloud cover in the temperate climates.

The spatial distribution of bare sands is the result of a complex of economic, social, and political factors. The deforestation of forest vegetation and the collapse of the Sadova–Corabia irrigation system (1970), followed by its abandonment, the clearing of forestry fences, and excessive land subdivision following the change in legislation (1989) influenced the distribution of bare sands. Two areas with bare sands persisted during the analysis period: one in the NW, where sand protection measures were lacking, and another in the SE, in the area of levelled dunes. The national reforestation programmes of recent years (2019) have decreased the areas with bare sands, especially in the south.

NSI accurately captures sandy areas, but it is limited to the transition to Chernozems by humus content and agrotechnical works.

This study will continue in the future because the area studied is vulnerable to climate change.