Characterization of short-term salinity fluctuations in the Neretva River Delta situated in the southern Adriatic Croatia using Landsat-5 TM

Ecological Indicators 110 (2020) 105924

Contents lists available at ScienceDirect

Ecological Indicators journal homepage: www.elsevier.com/locate/ecolind

Original Articles

Characterization of short-term salinity fluctuations in the Neretva River Delta situated in the southern Adriatic Croatia using Landsat-5 TM

T

⁎

Ivan Racetina, , Andrija Krtalicb, Veljko Srzicc, Monika Zovkod a

Department of Geodesy and Geoinformatics, University of Split, Faculty of Civil Engineering, Architecture and Geodesy, Matice hrvatske 15, 21 000 Split, Croatia Department of Photogrammetry and Remote Sensing, University of Zagreb, Faculty of Geodesy, Kaciceva 26, 10 000 Zagreb, Croatia c Department of Hydrotechnical Engineering, University of Split, Faculty of Civil Engineering, Architecture and Geodesy, Matice hrvatske 15, 21 000 Split, Croatia d Department of Soil Amelioration, University of Zagreb, Faculty of Agriculture, Svetosimunska cesta 25, 10 000 Zagreb, Croatia b

A R T I C LE I N FO

A B S T R A C T

Keywords: Agriculture Irrigation Multitemporal analysis Remote sensing Spectral indices

Elevated soil and water salinity can lead to increased risk of crop yield and quality reduction as well as land degradation. This is especially emphasised in arid and semi-arid regions with high levels of evaporation and irrigation practice. In this study, we assessed the potential of applying Landsat-5 Thematic Mapper (TM) to evaluate low to moderately salt-affected agricultural areas of the Neretva River Delta, situated in the southern Adriatic region of Croatia. Monthly measurements of water electrical conductivity (ECw), sampled from drainage/irrigation canals, were used as in situ salinity indicators. The relationship between agricultural plots reflectance (32 Landsat-5 TM images, distributed over 22 months from May 2009 to October 2011) and water salinity (measured at ten monitoring stations) was analysed using six Landsat-5 TM bands and spectral indices. Correlation analysis, simple (SLR) and multiple linear regression (MLR) models were implemented to predict ECw values using the satellite data. Very strong correlations were achieved in the spatial domain for August 2009, September 2010 and September and October 2011, while correlations in the temporal domain had relatively lower values. The Landsat-5 TM data managed to explain 43% of ECw variance using the SLR model, and 62% of the variance using the MLR model with three regressors. Landsat-5 TM remotely sensed data, paired with the in situ observations, showed a promising potential in enhancing the quality of environmental monitoring of low to moderately salt-affected agricultural fields in the irrigated deltas.

1. Introduction Soil salinization adversely affects ecosystems, crop production, and the economy, particularly in arid and semi-arid regions with high evaporation and irrigation practice levels (Daliakopoulos et al., 2016; Metternicht and Zinck, 2008; Pitman and Läuchli, 2004; Rasool et al., 2013; Szabolcs, 1979). Soil salinization is a consequence of human action, while soil salinity is a natural phenomenon. Therefore, sources of salinity can be categorized as primary, or naturally occurring, and secondary, or human-induced (Metternicht and Zinck, 2003). Primary salinity is caused by weathering of rocks, capillary rise from groundwater, seawater intrusion, salt blown by sea winds, and restrained drainage (Rasool et al., 2013; Szabolcs, 1979). Improper management of natural resources, which results in soil and water salt build-up, relates to secondary salinization and includes the inappropriate management of irrigated land, industrial wastewater, overuse of fertilizers, and land use changes (Rasool et al., 2013). The Food and Agriculture

⁎

Organization of the United Nations (FAO, 2017) states that, of 230 million hectares of irrigated land, almost 20% are salt-affected soils, while other authors report that this proportion is as high as 50% (Pitman and Läuchli, 2004). Salinization affects approximately 3.8 million hectares in Europe, and the cost of soil degradation due to salinization processes is estimated to be up to 321 million euros per year (Stolte et al., 2016). At the field scale, soil salinity is principally assessed by field observation and laboratory analysis of different parameters related to salinity level assessment (e.g. saturated soil-paste electrical conductivity (ECe), pH, cation exchange capacity, sodium adsorption ratio, concentration of different hazard ions) as well as related parameters in surface and ground water. Soil and water salinization monitoring does not only include identifying the places where salts accumulate, but also the detection of temporal and spatial variations (Ghassemi et al., 1995; Metternicht and Zinck, 2008). Achieving this, by using only conventional methods, is time-consuming and expensive. Remote sensing

Corresponding author. E-mail addresses: [email protected] (I. Racetin), [email protected] (A. Krtalic), [email protected] (V. Srzic), [email protected] (M. Zovko).

https://doi.org/10.1016/j.ecolind.2019.105924 Received 29 March 2019; Received in revised form 6 October 2019; Accepted 9 November 2019 1470-160X/ © 2019 Elsevier Ltd. All rights reserved.

Ecological Indicators 110 (2020) 105924

I. Racetin, et al.

(“Procjena rizika od katastrofa za Republiku Hrvatsku [Disaster Risk Assessment in Croatia]” 2015). The NRD is criss-crossed by a dense network of canals used for drainage and irrigation. During irrigation periods, salt which is leaching from the soil and salt from groundwater (driven by capillary rise) accumulates in these canals. The planned irrigation system in the NRD intended to draw water from upstream area and transport it by underground structures to the valley. Unfortunately, the construction work has not been finished completely, resulting in insufficient amount of freshwater pumped into the canals to reduce salt concentrations to desirable levels. Hence, the risk of soil salinization can be high (Zovko et al., 2018) and can jeopardize agriculture production and biodiversity if mitigation measures are not applied in timely manners. According to six-year monitoring data, almost 40% of the agricultural area in the NRD is saline (ECe > 2 dS m−1), while the average water electrical conductivity (ECw) is 9.3 dS m−1 for groundwater and 3.8 dS m−1 for surface waters (“Procjena rizika od katastrofa za Republiku Hrvatsku [Disaster Risk Assessment in Croatia]” 2015). Using conventional field and laboratory methods, water and soil salinity in the NRD has been monitored for more than two decades (Klačić et al., 1998; Romić et al., 1999; Zovko et al., 2018). According to the study of Lobell et al. (2007) it might be adequate to assume that soil salinity remains reasonably stable within a period of 5–7 years. In a single season, however, salts redistribute frequently within the soil profile (Scudiero et al., 2016). In general, salinity in surface waters and soil profiles in the NRD showed similar patterns with overall salinity increasing after the irrigation period (September/October) and restoring to background levels after the rainy period (March/April) (Romić et al., 2014), but with average salinity remaining fairly stable in the period from 2009 to 2013. Therefore, ECw of surface water samples from drainage/irrigation canals could be applied as the in situ benchmark for monitoring short-term soil salinity fluctuations in irrigated deltas. The main goal of this study was to examine the potential of multitemporal Landsat-5 Thematic Mapper (TM) reflectance data to characterize short-term soil salinity variations in the NRD. The relationship between soil salinity and irrigation (drainage) canals surface water salinity in the NRD was analysed using the data from 2010 and 2011. The interdependence of agricultural plots reflectance and surface water salinity was then analysed over a three-year period, using six Landsat-5 TM bands and selected intensity, salinity, soil and vegetation indices.

provides complementary data, rapidly covers broad areas, costs less and offers systematic periodic spatiotemporal sampling. Although remote sensing provides tools to characterize soil salinity, it is limited by several factors: atmospheric effects; spatial (including vertical), spectral and temporal resolution, and the physicochemical properties of soil and vegetation (Ben-Dor, 2002; Farifteh et al., 2006; Metternicht and Zinck, 2003). That is why synthesising in situ measurements with remote sensing can potentially optimize resources needed for mapping and monitoring of salinity (Bastiaanssen et al., 2000; Eldeiry and Garcia, 2008; Scudiero et al., 2017; Vermeulen and van Niekerk, 2016). The presence of salinity using remote sensing methods can be detected directly on bare soil surface, or indirectly by the presence of halophytic plants, changes in crop performance or from moisture condition (Allbed and Kumar, 2013; Ben-Dor, 2002; Mougenot et al., 1993). Direct approach used in the detection of salt accumulation is mainly based on indicators such as salt crust, efflorescence, bare soil patches with or without salt crystals, puffy soil and black stains (McGhie, 2005; Metternicht and Zinck, 2003). On the other hand, the use of different land cover responses to salinity as a proxy for detection and mapping salinity in an area might be considered an indirect approach (Elhag and Bahrawi, 2016; Muller and van Niekerk, 2016a; Scudiero et al., 2014; Van Der Werff et al., 2008; Zhang et al., 2012). The presence of halophytic plants could be an useful indirect indicator in areas with high levels of salinity (Dehaan and Taylor, 2002; Fernández-Buces et al., 2006; Metternicht, 1998; Zhang et al., 2011), but these areas are of less economic significance than irrigated agricultural areas, which can be harmed by slight or moderate quantities of salts (Ayers and Westcot, 1985; Machado and Serralheiro, 2017) that seasonally accumulate in the soil root zone. In irrigated areas, salt accumulation in the soil root zone is commonly a consequence of poor irrigation management (saline water, non-adequate salt leaching and drainage), and capillary rise of saline groundwater. Increase of salt content in the soil reduces its osmotic potential, which impedes the movement of water from soil into the root. In salt affected rhizosphere adverse effects on soil, plant growth, crop quality and yield could occur (Ayers and Westcot, 1985; Rhoades et al., 1992). The rate of these adverse effects depends on crop varieties, microclimatic conditions as well as agricultural management. Reduced water uptake slows the rate of plant growth, and may lead to changes in leaf colour to darker green, bluish-green, yellow or brown colour, rapid senescence and wilting. From a spectral point of view, vegetation stress (including salinity) causes increased reflectance of leaves in the visible spectrum and decreased reflectance in the near-infrared part of electromagnetic spectrum (Mulla, 2013; Zhang et al., 2011). Salinity assessment in low to moderately salt-affected irrigated agricultural areas is still a challenge for remote sensing (Lenney et al., 1996; Muller and van Niekerk, 2016a; Scudiero et al., 2017; Tamas and Lenart, 2006; Zhang et al., 2011). In majority of recent studies (Allbed and Kumar, 2013; El Harti et al., 2016; El-Battay et al., 2017; Hamzeh et al., 2013; Lugassi et al., 2017; Rahmati and Hamzehpour, 2017; Zhang et al., 2017) different satellite indices were implemented to achieve better models or correlations with ground truth data depicting soil salinity (e.g. ECe). At the irrigation scheme level, specific problems for remote sensing can also arise; poor farming practices and land mismanagement can lead to reduced plant vigour, which can be incorrectly attributed to salinity effects (Muller and van Niekerk, 2016b), while using vegetation indices in areas with sparse vegetation canopy (e.g. young orchards) can affect the index values negatively and result in incorrect estimates of salinity (Douaoui et al., 2006; Scudiero et al., 2017). The alluvial Neretva River Delta (NRD) is situated in the southern Adriatic region of Croatia (Fig. 1). From a global point of view, the NRD is small delta which covers just around 12,000 ha, but in the Mediterranean region such small deltas are precious environments for agricultural production as well as biodiversity conservation. In the NRD, agricultural production is of great economic importance, and the value of the economic output is estimated at around 49 million euros

2. Data and methods 2.1. Study area The research presented in this paper was performed in the alluvial NRD along the south-eastern Adriatic coast of Croatia, between 2009 and 2011 (Fig. 1). The observed area extends from 17.47° E to 17.65° E and 42.97° N to 43.05° N and has an elevation range from −1 to 5 m above mean sea level. In its lower course, the Neretva River forms a delta which covers 12,000 ha, of which around 5200 ha are agricultural land. The NRD is a semi-arid area characterised by Mediterranean climate with hot, dry summers and mild, wet winters. The mean annual rainfall (1980–2000) is 1230 mm, mostly in the period from October to April. The average annual air temperature is 15.7 °C and the annual PenmanMonteith reference evapotranspiration is 1196 mm. The highest values of both parameters, 25.2 °C and 191 mm respectively, occur in July, which is the month with the greatest demand for water (Romic et al., 2012). 2.2. In situ data The salinity monitoring system in the NRD is organized as integrated monitoring of surface water, groundwater bodies and soil (Romic et al., 2012; Romić et al., 2014; Zovko et al., 2013). Soil salinity 2

Ecological Indicators 110 (2020) 105924

I. Racetin, et al.

Fig. 1. Geographical extent of the study area and the locations of surface water quality sampling stations.

indicators were extracted at five monitoring stations, labelled S-1 to S-5 in Figs. 1 and 2. Soil monitoring stations established in 2009 were positioned to encompass a variety of representative soil types and different structure of agricultural land use in the NRD. Location S-1 was not tilled until 2010 and in 2011 a corn cultivation started, station S-2 is located at parcel used for agricultural production with alteration of spring or winter cabbage and melons or watermelons. Stations S-3 and S-4 are placed at parcels used for intensive production of mandarins (Citrus reticulata L.). Location S-5 is placed on the parcel which is used for conventional production of vegetables (alteration of brassicas and fruiting vegetables). Soil samples were collected and analysed twice a year, at the end of the winter season (late March or early April), and at the end of the summer/dry season (late September or early October). Soil salinity monitoring at the end of the dry season was used to determine salt accumulation in a soil profile after the irrigation period. Soil salinity monitoring at the end of the rainy season was used to evaluate soil capacity for preservation of equilibrium state of salts in a soil profile. Physical and chemical soil parameters were analysed per four soil layers with the step of 25 cm. ECe was analysed in the soil saturated paste extract, following the methodology proposed by Rhoades (1996), and using soil samples collected in March and August 2010 and April and October 2011. According to the World reference base (WRB) classification (IUSS Working Group WRB, 2015) soils at

stations S-1, S-3 and S-5 were classified as Gleysols and soils at stations S-2 and S-4 as Fluvisols. Monitoring of surface waters was organised at thirteen sampling stations in the NRD (Fig. 1) and determined by chemical laboratory analysis of water pH, ECw, and major ion concentrations (Romić et al., 2014). ECw is the parameter used for quantifying water salinity (the overall measure of dissolved salts) and is the in situ value used in further analyses. Surface water samples were collected once a month to determine seasonal changes and spatial trends of salinity fluctuations. The NRD is divided in six monitoring subregions, which have been previously classified as areas with highest risk of salinization (Romić et al., 2014). Surface water sampling stations were placed in lateral canals, melioration canals, main canals near pumping stations, at the mouth of effluent Mala Neretva river, and in the main Neretva River watercourse (near the town of Metkovic, about 21 km upstream from the Adriatic Sea) (Fig. 1). Two sampling locations in the lateral canals and the one in Mala Neretva were excluded from further data processing because there were no agricultural fields of interest affected by their values, therefore data from ten surface water sampling stations were used in this analysis. In June 2009, no ECw data were available for sampling stations 2 and 3, which is why there are only eight (instead of ten) reference ECw values for this time period.

3

Ecological Indicators 110 (2020) 105924

I. Racetin, et al.

Fig. 2. Close-up view of soil quality sampling stations and corresponding surface water sampling stations displayed on top of digital orthophoto (courtesy of the State Geodetic Administration, Republic of Croatia).

2.3. Satellite image acquisition and pre-processing

Table 2 Acquisition dates and WRS-2 coordinates of Landsat-5 TM images used in the study.

The Landsat-5 TM sensor has seven bands covering the visible, nearinfrared, shortwave infrared, and thermal parts of electromagnetic spectrum. The thermal band has a spatial resolution of 120 m, while all other bands have a spatial resolution of 30 m (Table 1). Only bands with a spatial resolution of 30 m were considered in this study. This study utilized 32 Landsat-5 TM Level-1 satellite images acquired during three consecutive years (2009–2011). All Landsat images are courtesy of the U.S. Geological Survey (USGS) and are publicly available through the USGS Center for Earth Resources Observation and Science (see http://earthexplorer.usgs.gov). Ten images belong to path 187 and row 30 of the Worldwide Reference System-2 (WRS-2) path/ row system and the remaining 22 images are from path 188 and row 30 (Table 2). Although it would have been better to use only one path/row combination for the time-series pixel analysis, the excess cloud cover (especially in July and August 2010) made the 188/30 images

Acquisition date of Landsat-5 TM image Path/Row

Spectral range (micrometers)

Spatial resolution (meters)

Band 1 – Blue Band 2 – Green Band 3 – Red Band 4 – Near Infrared (NIR) Band 5 – Shortwave Infrared (SWIR) 1 Band 6 – Thermal Band 7 - Shortwave Infrared (SWIR) 2

0.45–0.52 0.52–0.60 0.63–0.69 0.76–0.90 1.55–1.75

30 30 30 30 30

10.40–12.50 2.08–2.35

120 30

188/30

188/30

17 March 2010 18 April 2010 23 July 2010 8 August 2010 24 August 2010 21 April 2011 7 May 2011 24 June 2011 10 July 2011 28 September 2011

24 May 2009 9 June 2009 25 June 2009 27 July 2009 12 August 2009 28 August 2009 29 September 2009 15 October 2009 19 January 2010 9 April 2010 16 September 2010

2 October 2010 7 February 2011 11 March 2011 14 May 2011 30 May 2011 15 June 2011 17 July 2011 2 August 2011 18 August 2011 3 September 2011 5 October 2011

unusable, and we opted to also include the 187/30 images. This proved to be somewhat problematic because of the relationship between the projected Landsat pixel area (900 m2) and relatively small field sizes in the NRD. The pixel footprints for images from two different rows are not perfectly aligned, so they do not cover the exact patch of the Earth’s surface. This resulted, along with the spectral mixing of classes, in considerable differences in spectra at any given point for two subsequent images from different rows. In 2016, the USGS introduced Landsat Collection 1 Tiers: Tier 1, Tier 2 and Real-Time. The Tiers are used to sort Landsat data based on the quality and level of processing. Tier 1 offers Landsat images of the highest quality, which are suitable for time-series pixel analysis. As time-series data were explored in this study, we used only Tier 1 products. Tier 1 processing level includes Level-1 Precision and Terrain (L1TP) corrected data. L1TP data are radiometrically calibrated and

Table 1 Spectral range and spatial resolution of Landsat-5 TM bands (Markham and Barker 1985). Band

187/30

4

Ecological Indicators 110 (2020) 105924

I. Racetin, et al.

Fig. 3. Comparison of Landsat-5 TM data and digital orthophoto (courtesy of the State Geodetic Administration, Republic of Croatia). The images in one row represent exact geographical extent, but the digital orthophoto has ground sample distance (GSD) of 0.5 m and Landsat-5 TM has GSD of 30 m. Green dots represent locations of sampling points of Landsat-5 TM and gauge glass icon represents surface water sampling station. In the first row, all of the 9 Landsat-5 TM sampling points used for calculating reference spectra for surface water station number 5 are presented. Grid coordinates are expressed in the UTM zone 33N (EPSG:32633).

the agricultural field, which is connected with the canal where the surface water was sampled. Selection of representative Landsat sampling points was done using a digital orthophoto with the ground sample distance (GSD) of 0.5 m (Fig. 3). This resulted with at least 9 and up to 25 Landsat sampling points for an area. For each image, the reference spectrum of surface water sampling station was calculated as the average value of corresponding Landsat sampling points. In the months with multiple Landsat images, since the ground truth values of ECw were collected once a month, the average value for these images was calculated as the final reference spectrum. Next, data from 32 Landsat images (Table 2) were transformed into 218 reference spectra (22 months for 10 ECw stations, minus two because of unavailable data for stations 2 and 3 in June 2009). These steps were carried out using

orthorectified using GCPs and a DEM to correct for elevation displacement. Satellite data were corrected for radiometric and atmospheric effects using Fast Line-of-sight Atmospheric Analysis of Spectral Hypercubes (FLAASH) atmospheric correction (ENVI, 2018). FLAASH incorporates the Moderate Resolution Transmission (MODTRAN) 4 radiative transfer code (Matthew et al., 2000). Surface reflectance values were calculated using the FLAASH algorithm incorporated in the ENVI (Exelis Visual Information Solutions, Boulder, Colorado) software. A great diversity of fruit and vegetable cultures is grown in the NRD area. To capture this diversity and to acquire reference spectra, Landsat sampling points were defined in the agricultural fields near surface water sampling stations. The criterion for selection of Landsat sampling points was that pixel footprint on the ground should completely cover 5

Ecological Indicators 110 (2020) 105924

I. Racetin, et al.

spectrum.

Table 3 Spectral indices evaluated in this study. B stands for Landsat-5 TM band and the subscript numbers correspond to the band numbers shown in Table 1.

1 2 3

Indices

Equation

Normalized Difference Vegetation Index Enhanced Vegetation Index

NDVI =

B4 − B3 B4 + B3

EVI = 2.5

B4 − B3 B4 + 6B3 − 7.5B1 + 1

4

Soil Adjusted Vegetation Index Ratio Vegetation Index

RVI =

5

Difference Vegetation Index

DVI = B4 − B3

6

GARI =

B4 − B2 + 0.9(B1 − B3 ) B4 + B2 + 0.9(B1 − B3 )

GVMI =

(B4 + 0.1) − (B5 + 0.02) (B4 + 0.1) + (B5 + 0.02)

9

Green Atmospherically Resistant Vegetation Index Global Vegetation Moisture Index* Normalized Difference Salinity Index Salinity Index 1

S1 =

B1 B3

10

Salinity Index 2

S2 =

B1 − B3 B1 + B3

11

Salinity Index 3

S3 =

B2 × B3 B1

12

Salinity Index 4

S4 =

13

Salinity Index 5

S5 =

B1 × B3 B2

14

Salinity Index 6

S6 =

B3 × B4 B2

15

Salinity Index 7*

SI7 =

B5 B7

16

Salinity Index 8*

SI8 =

B4 − B5 B4 + B5

17

Salinity Index

SI =

B3 × B4

18

Salinity Ratio

SR =

19

Vegetation Soil Salinity Index ASTER Salinity Index

VSSI = 2B2 − 5(B3 + B4 )

7 8

20 21 22

Canopy Response Salinity Index Brightness Index

23

(B4 − B3 )(1 + L) B4 + B3 + L

SAVI =

B4 B3

B3 − B4 B3 + B4

NDSI =

B1 × B3

B3 − B4 B2 + B 4

SIASTER = CRSI =

2.5. Correlation and regression analysis

Reference

B5 − B7 B5 + B7 (B4 × B3) − (B2 × B1 ) (B4 × B3) + (B2 × B1 )

BI =

B32 + B42

Intensity Index 1

INT1 =

B2 + B3 2

24

Intensity Index 2

INT2 =

B2 + B3 + B4 2

25

Soil Index 1

SI1 =

B2 × B3

26

Soil Index 2

SI2 =

B22 + B32 + B42

27

Soil Index 3

SI3 =

B22 + B32

Pearson correlation coefficients (r) were calculated between ECw values and reference Landsat-5 TM spectra and spectral indices. Pearson correlation coefficient, or Pearson’s product moment correlation coefficient is used for quantifying linear relationships between two variables. Correlation analysis was performed by averaging data in the time domain and examining the relationship between ECw and satellite data. Next, the data were averaged by surface water sampling stations so that temporal ECw and satellite data relationships could be explored. In addition to the basic correlation analysis, the simple and multiple linear regression models (Freedman, 2009) for prediction of ECw values using the Landsat-5 TM bands and spectral indices were performed. The simple linear regression (SLR) model uses one continuous independent variable (regressor, x) to predict one continuous dependent variable (response, y), which is related by the equation (Freedman, 2009):

Rouse Jr et al. (1974) Huete et al. (1997) Huete (1988) Birth and McVey (1968) Douaoui et al. (2006) Gitelson et al. (1996) Ceccato et al. (2002) Khan et al. (2001) Abbas and Khan (2007) Abbas and Khan (2007) Abbas and Khan (2007) Khan et al. (2001) Abbas and Khan (2007) Abbas and Khan (2007) Bannari et al. (2008) Bannari et al. (2008) Dehni and Lounis (2012) Dehni and Lounis (2012) Dehni and Lounis (2012) Bannari et al. (2008) Scudiero et al. (2014) Khan et al. (2001) Douaoui et al. (2006) Douaoui et al. (2006) Douaoui et al. (2006) Douaoui et al. (2006) Douaoui et al. (2006)

y = β0 + β1 x

(1)

where β0 and β1 are constant values called regression coefficients. The constant β0 is also called the intercept, and β1 is referred to as slope. The multiple linear regression (MLR) model, expands the SLR model with two or more independent continuous variables (xi) (Freedman, 2009):

y = β0 + β1 x1 + β2 x2 + …+βn x n

(2)

where β0, β1, β2, …, βn are regression coefficients. 3. Results The results obtained are presented in the following order. Descriptive statistics of the in situ measurements and remotely sensed data are displayed first. These sections provide an overview of the input data: water salinity expressed by ECw, soil salinity expressed by ECe, the relationship between water and soil salinity, and the spatiotemporal variations of Landsat-5 TM reflectance spectra. Next, the correlation analysis is divided into a spatial and temporal domain and examined separately. In the final section, the irrigation period (July-September) of analysed years (2009–2011) was compared with the whole year data, and the results of spatiotemporal correlation and regression analysis are presented. 3.1. Descriptive statistics of surface water quality measurements ECw values in the NRD during the analysed period ranged from 0.31 to 8.58 dS m−1, while the average ECw per sampling station ranged from 0.50 to 5.78 dS m−1 (Fig. 4). Standard deviation (SD) of ECw per sampling station was from 0.19 dS m−1 (location 1) to 1.48 dS m−1 (location 2). Using the classification of saline waters (Rhoades et al., 1992), average water quality at location 1 can be classified as nonsaline, at locations 7 and 9 as slightly saline, and at the remaining locations as moderately saline. Temporal fluctuation of salinity in the NRD (Fig. 5), revealed that the maximum average ECw value of 4.38 dS m−1 occurred in October 2009, while the lowest average value of 1.96 dS m−1 was in January 2010. The largest SD and range of ECw values happened in June 2010, with ECw values of 2.49 dS m−1 and 8.23 dS m−1, respectively. The minimum range (3.81 dS m−1) was in February 2011 and the lowest SD (1.04 dS m−1) in January 2010.

* Equation adjusted to the Landsat-5 TM bands.

the ESRI (Environmental Systems Research Institute, Redlands, California) ArcGIS software. 2.4. Spectral indices Spectral indices are dimensionless measures that act as indicators of properties of interest, such as water, vegetation, soil, salinity, and intensity. Based on previous research on salinity using remote sensing (Abbas et al., 2013; Bannari et al., 2008; Dehni and Lounis, 2012; Douaoui et al., 2006), we selected a number of differences indices (Table 3): vegetation (1–7), salinity (8–21), intensity (22–24) and soil (25–27). Soil and intensity indices were included in the analysis because the vegetation in the NRD does not have a dense canopy, and soil background scatter was expected to influence the Landsat-5 TM

3.2. Descriptive statistics of soil quality measurements Among the four soil sampling events, average ECe in the soil profiles ranged from 0.62 dS m−1 to 3.90 dS m−1 (Fig. 7). The highest average ECe (3.66 dS m−1) was exhibited at location S-3 and the lowest 6

Ecological Indicators 110 (2020) 105924

I. Racetin, et al.

classified as non saline (ECe: 0–2 dS m−1). Soil salinity data for locations S-2, S-3 and S-5 showed seasonal behaviour: an increase in salinity after the irrigation period (AugustOctober), and its decrease after the rainy period (March-April). Locations S-1 and S-4, however, did not follow that pattern, due to differences in soil characteristics as well as agricultural management. 3.3. Relationship between surface water and soil quality measurements Surface water quality in the NRD was analysed once a month at thirteen locations, while soil quality data was collected twice a year at five locations. Seasonal changes in soil salinity fluctuations are slower comparing to surface water salinity fluctuations, hence the frequency and number of sampling locations between these in situ measurements is considerably different. Precipitation in the days preceding surface water sampling influences the water quality parameters. This can be seen in the sudden changes in ECw values for consecutive months (Fig. 6). To smoothen the variation in the data, the moving average method (Chatfield, 2004) with the period of three months was applied to the ECw data. When the trends in the smoothed ECw data are compared with the ECe data, by relating the measurements of soil salinity to the measurements of surface water salinity at the nearest surface water sampling location (Figs. 1 and 2), it can be seen that they follow similar pattern at locations S-2, S-3 and S-5 (Fig. 7). Four ECw measurements, taken prior to the month of ECe sampling, were used to define the trends for August 2010 and October 2011. For April 2011 the ECw measurements from October 2010 to January 2011 were selected instead because the ECw values sharply increased in March and April 2011 (Fig. 7). The main cause for this increase was the exceptionally small amount of precipitation in that time period. Meteorological station located in the town of Opuzen (Fig. 1) measured rainfall of

Fig. 4. Descriptive statistics of ECw in the NRD for the analysed period averaged per surface water sampling station. Error bars indicate one standard deviation (SD) above and under average values.

(0.85 dS m−1) at location S-1. The lowest average of ECe at five sampling locations (1.46 dS m−1) was detected in April 2011, while the highest average value of 2.44 dS m−1 was measured in August 2010. Using the general classification of soil salinity (Abrol et al., 1988) station S-3 was at upper level of slightly saline class (ECe: 2–4 dS m−1) in all four sampling occasions. Stations S-2 and S-5 showed slight salinity in August 2010. In the remaining observations soil salinity can be

Fig. 5. Descriptive statistics of temporal fluctuations of ECw in the NRD (average monthly ECw at ten surface water monitoring stations). Error bars indicate one standard deviation above and under average values. 7

Ecological Indicators 110 (2020) 105924

I. Racetin, et al.

Fig. 6. Example of smoothening of measured ECw values using the moving average method with the period of 3 months for surface water sampling location 9.

118.3 mm in March and only 31.2 mm in April 2011. Most of the precipitation occurred during a few days: 46.8 mm fell on March 17th and 29.5 mm on March 29th, while 28.0 mm fell on April 14th, 2011.

section, the level of approximation is quantitatively described with SD and coefficient of variation (CV). Mean values and average SDs needed for determination of CV were calculated as average values of the reference spectra and corresponding SD per ECw station. An example of averaged Landsat-5 TM reflectance values (reference spectra) for surface water sampling station number 10 in March 2010, with the corresponding error bar (one SD), is shown in Fig. 8a. Fig. 8b depicts the CV of reference spectra by station. Band 1 showed the lowest SD of averaged reflectance values for all surface

3.4. Descriptive statistics of satellite data The reference spectra were determined by averaging Landsat-5 TM reflectance values. That introduced an approximation in the further analysis of the relationship between ECw and satellite data. In this

Fig. 7. Values of average salinity in a soil profile (ECe, marked by the lines with the diamonds) and the corresponding smoothed surface water salinity (smoothed ECw, measured from the nearest surface water sampling location, marked by the dotted lines with the dots) for the period from June 2009 to December 2011. 8

Ecological Indicators 110 (2020) 105924

I. Racetin, et al.

Fig. 8. (a) Reference spectra of Landsat-5 TM bands for ECw station 10 in March 2010. Error bars indicate one standard deviation above and under reference spectra. (b) Coefficient of variation of reference spectra (Landsat-5 TM Bands) per surface water sampling station.

water sampling stations. The same band had the lowest average reflectance values at all surface water sampling stations, and it resulted in the highest values of CV. The highest values of SDs at ECw stations 2–7 and 10 occurred in Band 4, at stations 1 and 8 in Band 7, and at station 9 in Band 5. The lowest average CV was noted in Band 4 and Band 5, while the lowest CV was at stations 2 and 5. 3.5. Correlation analysis in the spatial domain Correlation analysis in the spatial domain was carried out by relating the 33 variables, six Landsat-5 TM bands (bands 1–5 and band 7, Table 1) and 27 spectral indices (Table 3), to the ECw values for a given timestamp (Table S1 in the Supplementary Material). In other words, at one point in time the relationship among the ten surface water sampling stations was observed (Fig. 9). The SI8 index achieved the highest correlations for each observation year (2009–2011): in August 2009 (−0.88), September 2010 (−0.85), and September 2011 (−0.93) (Fig. 9, Table 3). Regarding the Landsat-5 TM data, the bands from near-infrared and shortwave-infrared spectra (marked B4, B5 and B7 in Fig. 10) showed the highest positive correlations. The SI8 index and Global Vegetation Moisture Index (GVMI) are basically the same indices, the only difference being in GVMI having small correction factors. These two indices combined showed the lowest negative correlation for seven months, all of them distributed in the warm part of a year: August and September of 2009 and 2011 and July, September and October 2010. Each of the Band 5 and S2 index yielded the highest positive correlation on five occasions (Fig. 10). It should be noted that Band 5, which covers shortwave infrared (SWIR) region (1.55–1.75 µm), had the highest correlations in the August and September 2009, as well as September and October 2010.

Fig. 9. Scatterplot of ECw and SI8 index values for September 2011 for ten ECw sampling stations. Numbers near the circles on the plot correspond to the numbers of ECw sampling stations shown on the Fig. 1. Dotted line represents the linear trendline and Pearson correlation coefficient r is equal to −0.93.

The highest correlation coefficient between ECw reference values and Landsat-5 TM bands (B1-B7) and spectral indices was 0.49 (VSSI index) at station 8. The lowest observed correlation value of −0.50 (S3 index) was also at station 8 and the same value at station 6 for S6 index (Fig. 12). Stations 5 and 10 showed the weakest correlation, having maximum correlation of 0.18 (NDSI index) and 0.22 (S2 index), and minimum correlation of −0.20 (GVMI index) and −0.22 (S3 index) respectively (Fig. 12).

3.6. Correlation analysis in the temporal domain For each sampling station, the ECw values were correlated with the 33 variables using the data for 22 months analysed (Fig. 11, Table S2 in the Supplementary Material). Using this approach, it was possible to examine the relations between temporal fluctuations of water salinity and satellite data at specified locations.

3.7. Spatiotemporal correlation and regression analysis The months from July to September (2009–2011) were selected for spatiotemporal analysis because they mark the irrigation period in the study area. When the data for all 22 months were used (218 reference 9

Ecological Indicators 110 (2020) 105924

I. Racetin, et al.

Fig. 10. Maximum and minimum Pearson correlation coefficient between ECw values and 33 variables for each month analysed. Statistical significance is considered at p < 0.05 (*), p < 0.01 (**) and p < 0.001 (***).

Fig. 11. Scatterplot of ECw and S3 index values for 22 months for ECw sampling station number 8. Dotted line represents the linear trendline and Pearson correlation coefficient r is equal to −0.50.

spectra), a minimum negative correlation of −0.32 was achieved for GVMI. S1 and S2 indices scored the maximum positive correlation of 0.26, while the average of absolute Pearson correlation coefficients of ECw with all 33 variables was 0.13 (Fig. 13, Table S3 in the Supplementary Material). During the irrigation period, the average correlation more than doubled, from 0.13 to 0.31 (Table S4 in the Supplementary Material). The maximum positive correlation increased from 0.26 for S2 index to 0.42 for Landsat-5 TM Band 5, while the

Fig. 12. Maximum and minimum Pearson correlation coefficient between ECw values and 33 variables for each sampling station. Statistical significance is considered at p < 0.05 (*), p < 0.01 (**) and p < 0.001 (***).

minimum negative decreased from −0.32 for GVMI to −0.65 for SI8 index. The similar values of ECw correlation with GVMI and SI8 indices were due to their high correlation, and the only difference between SI8 10

Ecological Indicators 110 (2020) 105924

I. Racetin, et al.

Fig. 14. ECw-SI8 index scatterplot for the period July-September (2009–2011), with single linear regression model and coefficient of determination (R2).

checked the multicollinearity using the Variance Inflation Factor (VIF). The VIF scores for the regressors in the proposed MLR model (3) was 1.134 (S2), 1.512 (SI7) and 1.381 (SI8). If VIF score is less than 3.3, we can suggest no collinearity (Cenfetelli and Bassellier, 2009). Fig. 13. Pearson correlation coefficients between ECw values and 33 variables calculated using the data for all dates (22 months, 218 observations), and also using only the data for months from July to September (2009–2011) (90 observations). First two columns represent the average of absolute Pearson correlation coefficients for all 33 variables (Table S3 in the Supplementary Material). Then, four variables with the highest positive and lowest negative correlations are showed. Statistical significance is considered at p < 0.05 (*), p < 0.01 (**) and p < 0.001 (***).

4. Discussion In the spatial domain, the Global Vegetation Moisture Index (GVMI) and SI8 index showed very strong and the best overall correlation with the surface water salinity data (Fig. 10). These two indices are practically the same, so they can be observed as one. They both utilize a Band 4 from Near Infrared (NIR) region and a Band 5 from Shortwave Infrared (SWIR) region, only GVMI has additional constant correction factors (Tables 1, 3). The reason why exactly those indices achieved the best results for monitoring of spatial salinity variations can be speculated as follows: the NIR region is often related to the present amount of biomass amount and/or the soil and canopy contrast, while the SWIR region is related to the canopy (leaf) moisture content. As water strongly absorbs in SWIR region, low water content in leaves increases the leaf reflectance. Although determined correlation coefficient values were lower in temporal compared to spatial domain, some patterns in temporal domain could be observed. The highest positive correlation value, but only moderately strong, has been observed for the station closest to the sea. Regarding the spectral regions, our results demonstrate that the indices with the strongest correlations in the temporal domain were calculated using the data only from Visible region (S2 and S3) and combination of Visible and NIR region (S6 and Vegetation Soil Salinity Index – VSSI) (Fig. 12), unlike SWIR and NIR bands which held the most information about variations in the temporal domain. It can be noticed that some indices had the highest positive correlations for one station and the lowest negative correlations for the other one. For example, VSSI has the highest correlation at stations 2, 8 and 9 and the lowest at station 4 (Fig. 12). This is clear indicator of high influence of local conditions to surface reflectance temporal variations, especially when taking the whole year data. By comparing the spatial and temporal correlations, correlation coefficient values within the spatial domain (Fig. 10) were generally higher than the temporal ones (Fig. 12). This is not surprising since the presence of several additional factors affecting the agricultural plots reflectance spectra through time, such as crop rotation, phenological cycle of vegetation, meteorological conditions and land management practices. Very strong correlations are achieved in the spatial domain

and GVMI indices was in the correction factors of 0.1 and 0.02 for Bands 4 and 5, respectively (Table 3). Due to the significant increase in the correlation values for the irrigation period, SLR and MLR models for ECw were created. To define the models, 90 observations were used (data from ten stations for nine months). The SI8 index showed the greatest overall correlation with ECw, and was used for the SLR analysis (Table S4 in the Supplementary Material). The ECw-SI8 index scatterplot and SLR model are shown in Fig. 14. Using the defined SLR model (Fig. 14), the SI8 index explains 43% of the variance in ECw. The adjusted coefficient of determination (adjusted R2) was 0.42 and the standard error of regression was 1.40. The standard error of intercept was 0.30 and SI8 was 1.53, while Pvalues are less than 0.001. To build the MLR model, the best subsets method (Neter et al., 1996) was used with the adjusted R2 as the deciding parameter. The best subsets method works in a way that it creates all possible combinations of potential regressors for a given number of model dimensions. Adding one additional regressor to the SLR model (Fig. 14) increased the explained variance by 11%, while using two additional regressors explained 62% of the total variance of ECw. Further addition of regressors did not increase the explained variance significantly. The best MLR model with 3 regressors is defined by the equation:

EC w = 2.375 + 14.748 × S2 + 4.981 × SI7 − 16.876 × SI8 ,

(3)

where S2, SI7 and SI8 represent salinity indices as defined in Table 3. The coefficient of determination R2 was 0.62, while the adjusted R2 was 0.60. The standard error of the regression model was 1.16 dS m−1 and the standard errors of the intercept, S2, SI7 and SI8 were 1.29, 2.95, 0.87 and 1.49, respectively. The P-value for the intercept was 0.068, and for other regressors it was less than 0.001. For the proposed model we 11

Ecological Indicators 110 (2020) 105924

I. Racetin, et al.

Saltwater enters the NRD (Ljubenkov and Vranjes, 2012) by seawater intrusion through the riverbed (in the form of salt wedge) and from underground layers (where groundwater is already salinized because of direct contact with the sea via porous aquifer). Intensive agriculture requires large quantities of water, particularly during summers in arid and semi-arid climates Due to highly porous (carbonate) aquifer, sea water intrudes more than 10 km upstream in the NRD, which can lead to salinization of wells and drainage /irrigation canals. Additionally, water in drainage/irrigation canals in the NRD accumulates salts from the soil which are leached into the canals. Therefore, ECw of surface/ irrigation water should be regarded as a valuable alternative in situ measure of short-term soil salinity fluctuations in irrigated areas for the indirect remote sensing monitoring of salinity. This was confirmed for stations S-2, S-3 and S-5. Eventually, by substituting a monthly measurement of surface salinity with an average of daily measurements it would be possible to reduce the influence of precipitation on the data and achieve a more reliable estimate of overall salinity in the NRD. In August 2010, an exhaustive soil survey of the NRD was carried out by Zovko et al. (2018). The survey resulted in 246 soil samples of the shallow soil horizon (depth up to 25 cm), arranged in the regular grid of 500 by 500 m. In future studies, these in situ samples will be incorporated in investigations of the available satellite data for spatial salinity classification and mapping. Unfortunately, due to the Landsat-5 TM limitations and complex nature of the studied area, no very strong spatiotemporal correlation was achieved. The implementation of hyperspectral or multispectral satellites with better spatial and spectral resolution, and the analysis of longer time intervals, would possibly provide more precise and reliable data and models.

for August 2009, September 2010, and September and October 2011 (Fig. 10), while in the temporal domain only moderate correlations are reached for stations 2–4 and 6–9 (Fig. 12). Our results show that for indirect monitoring of salinity variations in the spatial domain, the most information is held in the indices which combine data from NIR and SWIR region and in the temporal domain for indices from Visible and NIR region. These findings could be useful and transferable to similar irrigated agricultural deltas in arid and semi-arid areas. The correlation between ECw and satellite data when using only the data from July to September (2009–2011) was more than double that of the full year (Fig. 13), confirming that the impact of saline surface waters (expressed by electrical conductivity) intensifies in this period. Using the data from ten surface water stations for July–September (2009–2011), which in total gives 90 observations, the strongest correlation was yielded with the SI8 index. The selected period represents the irrigation period in the Neretva River Delta (NRD) and approximately corresponds to the operation schedule of pumping stations, which regulate water tables in the NRD. The ground sample distance of the used Landsat-5 TM bands is 30 m, (one pixel covers an area of 900 m2) which can cause pixel values to be a mixtures of spectral signatures. In the NRD, due to small parcel sizes and agricultural practices, the problem of spectral mixtures for Landsat-5 TM bands was reflected in the increased coefficient of variation and standard deviation values of the reference spectra. In the limited NRD area (12,000 ha), soil substrates with various physicochemical properties are present. This influences the distribution of salinity in the soil profile (e.g. sandy soil texture impedes the capillary rise of salty water), as well as the plant and soil spectra (Zovko et al., 2018). In spite of these limitations, the Landsat-5 TM satellite data explained 43% of ECw variance using the Simple Linear Regression (SLR) model, and 62% using the Multiple Linear Regression (MLR) model with three regressors. The standard errors of regression were 1.40 dS m−1 for the SLR and 1.16 dS m−1 for the MLR model. Nevertheless, it was important for this study to address the problem of slight and moderate salinity, with the maximum recorded ECw value of 8.58 dS m−1 and the average value of only 3.12 dS m−1. Although the latter may seem negligible, the consequences in terms of recorded yield reductions are considerable, due to the sensitivity of specific crops grown in the NRD (Zovko et al., 2013). The achieved standard errors of regression (1.40 and 1.16 dS m−1) for prediction of ECw may be accurate enough if higher levels of salinity occur, but at present salinity levels and in relation to current agricultural crops they are of limited use. Agronomic production in the NRD mostly consists of citrus fruits (mandarins), peaches and apples (61%), followed by tomatoes, aubergines and peppers (17%). Other cultures include brassicas, strawberries, melons, watermelons, etc. (“Procjena rizika od katastrofa za Republiku Hrvatsku [Disaster Risk Assessment in Croatia]” 2015). Citrus fruits (mandarins), peaches, apples and strawberries are categorized as sensitive to salinity, while other agricultural crops fall into the moderately sensitive category (Ayers and Westcot, 1985). Therefore, with the salinity levels currently present in the NRD, except for surface water sampling station 1 (Fig. 4), these agricultural crops could be at risk of exhibiting yield reductions and, therefore, special care should be given to irrigation management to limit (or reduce) the agricultural as well as environmental damage. There is a general lack of research dealing with short-term salinity fluctuations (Scudiero et al., 2016), especially in irrigated areas. These areas are usually of high economic importance and great sensitivity to salinity hazard. In coastal irrigation schemes, such as the one in the lower Neretva River Valley, an extremely complex interaction between seawater, surface water, ground water, soil and human presence occurs. This presents a huge challenge to a wide variety of environmental scientists. In terms of the remote sensing approach, the addition of meteorological, geological and groundwater level data to the satellite data is needed to provide a more precise overview of salinity fluctuations. Climate, hydrological and hydrogeological processes, including engineering objects, control the overall salinity fluctuations in the NRD.

5. Conclusion Low to moderately salt-affected irrigated agricultural deltas, like the Neretva River delta, are exceptionally fragile landscapes which are often under risk of increased salinization and further degradation. The mosaic of extremely small parcels which are criss-crossed with different bodies of water: from seawater and natural surface water to artificial drainage canals (which can have different functions within the reclaimed area), pose a huge challenge for establishing high-quality environmental monitoring through in situ techniques which are usually organized at few sparse locations. Therefore, this research evaluates potential improvements in environmental monitoring through the use of Landsat-5 TM, as the addition of remotely sensed data offers complementary beneficial properties to the in situ data: it is collected over vast areas in a systematic way and regular intervals. The results indicate that there is weak correlation of multi-year Landsat-5 TM data and derived indices with ECw when whole year data are used. For the period of July to September (2009–2011), however, strong negative correlations with Salinity Index 8 and Global Vegetation Moisture Index were documented. This timing corresponds to the irrigation period in the NRD; pumping stations which regulate water tables in the NRD operate from mid-June until the end of September. Using the data from July to September (2009–2011) the Simple Linear Regression model with Salinity Index 8 managed to explain 43% of ECw variance, and the best Multiple Linear Regression model with three regressors: Salinity Index 2, Salinity Index 7 and Salinity Index 8 explained 62% of ECw variance. This indicates that Landsat-5 TM data, though perhaps not appropriate for year-round salinity monitoring, may prove useful for spatiotemporal monitoring of short-term salinity fluctuations. In future steps, freely available Sentinel-2, or perhaps Worldview series satellite data and additional in situ data (groundwater quality and levels, meteorological, crop type and soil properties) could be considered for salinity monitoring in the NRD. By increasing the spatial resolution of the data used, the problem of spectra mixing can be reduced. Implementation of linear mixture or linear subspace models in 12

Ecological Indicators 110 (2020) 105924

I. Racetin, et al.

the processing chain could potentially be used to tackle the same problem. A supervised classification approach should be evaluated to distinguish between various severity levels of salt effects. Ultimately, the high quality in situ and remote sensing monitoring data will enable decision-makers to implement timely measures to prevent degradation and to initiate the recovery of this environmentally and economically highly valued area.

Dehni, A., Lounis, M., 2012. Remote sensing techniques for salt affected soil mapping: application to the Oran region of Algeria. In: Mahmoudi, H., Merzouk, N.K., Spahis, N. (Eds.), Iswee'11. Elsevier Science Bv, Amsterdam, pp. 188–198. Douaoui, A.E.K., Nicolas, H., Walter, C., 2006. Detecting salinity hazards within a semiarid context by means of combining soil and remote-sensing data. Geoderma 134 (1–2), 217–230. https://doi.org/10.1016/j.geoderma.2005.10.009. El-Battay, A., Bannari, A., Hameid, N., Abahussain, A., 2017. Comparative study among different semi-empirical models for soil salinity prediction in an arid environment using OLI Landsat-8 data. Adv. Remote Sens. 6 (01), 23. El Harti, A., Lhissou, R., Chokmani, K., Ouzemou, J.-E., Hassouna, M., Bachaoui, E.M., El Ghmari, A., 2016. Spatiotemporal monitoring of soil salinization in irrigated Tadla Plain (Morocco) using satellite spectral indices. Int. J. Appl. Earth Obs. Geoinf. 50, 64–73. https://doi.org/10.1016/j.jag.2016.03.008. Eldeiry, A.A., Garcia, L.A., 2008. Detecting soil salinity in alfalfa fields using spatial modeling and remote sensing. Soil Sci. Soc. Am. J. 72, 201–211 10.2136/ sssaj2007.0013. Elhag, M., Bahrawi, J.A., 2016. Soil salinity mapping and hydrological drought indices assessment in arid environments based on remote sensing techniques. Geosci. Instrum. Methods Data Syst. Discuss. 1–31 doi: 10.5194/gi-2016-39. ENVI. 2018. Atmospheric Correction Module: QUAC and FLAASH User’s Guide. Accessed June 15. http://www.harrisgeospatial.com/portals/0/pdfs/envi/flaash_module.pdf. FAO. Salt-affected soil. Accessed 20-12-2017. http://www.fao.org/soils-portal/soilmanagement/management-of-some-problem-soils/salt-affected-soils/moreinformation-on-salt-affected-soils/en/. Farifteh, J., Farshad, A., George, R.J., 2006. Assessing salt-affected soils using remote sensing, solute modelling, and geophysics. Geoderma 130 (3–4), 191–206. https:// doi.org/10.1016/j.geoderma.2005.02.003. Fernández-Buces, N., Siebe, C., Cram, S., Palacio, J.L., 2006. Mapping soil salinity using a combined spectral response index for bare soil and vegetation: a case study in the former lake Texcoco, Mexico. J. Arid Environ. 65 (4), 644–667. https://doi.org/10. 1016/j.jaridenv.2005.08.005. Freedman, D.A., 2009. Statistical Models: Theory and Practice, 2 ed. Cambridge University Press. Ghassemi, F., Jakeman, A.J., Nix, H.A., 1995. Salinisation of Land and Water Resources: Human Causes, Extent, Management and Case Studies. CAB international. Gitelson, A.A., Kaufman, Y.J., Merzlyak, M.N., 1996. Use of a green channel in remote sensing of global vegetation from EOS-MODIS. Remote Sens. Environ. 58 (3), 289–298. https://doi.org/10.1016/s0034-4257(96)00072-7. Hamzeh, S., Naseri, A.A., AlaviPanah, S.K., Mojaradi, B., Bartholomeus, H.M., Clevers, J.G.P.W., Behzad, M., 2013. Estimating salinity stress in sugarcane fields with spaceborne hyperspectral vegetation indices. Int. J. Appl. Earth Obs. Geoinf. 21, 282–290. https://doi.org/10.1016/j.jag.2012.07.002. Huete, A.R., 1988. A soil-adjusted vegetation index (SAVI). Remote Sens. Environ. 25 (3), 295–309. https://doi.org/10.1016/0034-4257(88)90106-x. Huete, A.R., Liu, H.Q., Batchily, K., vanLeeuwen, W., 1997. A comparison of vegetation indices global set of TM images for EOS-MODIS. Remote Sens. Environ. 59 (3), 440–451. https://doi.org/10.1016/s0034-4257(96)00112-5. IUSS Working Group WRB, 2015. World reference base for soil resources 2014, update 2015: international soil classification system for naming soils and creating legends for soil maps. In: World Soil Resources Reports No. 106. FAO, Rome, pp. 192. Khan, N.M., Rastoskuev, V.V., Shalina, E.V., Sato, Y., 2001. Mapping salt-affected soils using remote sensing indicators—a simple approach with the use of GIS IDRISI. Paper Presented at the 22nd Asian Conference on Remote Sensing. Klačić, Ž., Romić, D., Romić, M., 1998. Salinity of agricultural soils in the river Neretva delta (Croatia). 1st Inter-regional Conference on Environment-water: Innovative Issues in Irrigation and Drainage. Lenney, M.P., Woodcock, C.E., Collins, J.B., Hamdi, H., 1996. The status of agricultural lands in Egypt: the use of multitemporal NDVI features derived from Landsat TM. Remote Sens. Environ. 56 (1), 8–20. https://doi.org/10.1016/0034-4257(95) 00152-2. Ljubenkov, I., Vranjes, M., 2012. Numerical model of stratified flow – case study of the Neretva riverbed salination (2004). Gradevinar 64 (2), 101–112. Lobell, D.B., Ortiz-Monasterio, J.I., Gurrola, F.C., Valenzuela, L., 2007. Identification of saline soils with multiyear remote sensing of crop yields. Soil Sci. Soc. Am. J. 71, 777 doi: 10.2136/sssaj2006.0306. Lugassi, R., Goldshleger, N., Chudnovsky, A., 2017. Studying vegetation salinity: from the field view to a satellite-based perspective. Remote Sens. 9 (2) doi: 10.3390/ rs90201022. Machado, R., Serralheiro, R., 2017. Soil salinity: effect on vegetable crop growth. management practices to prevent and mitigate soil salinization. Horticulturae 3, 30. https://doi.org/10.3390/horticulturae3020030. Markham, B.L., Barker, J.L., 1985. Spectral characterization of the LANDSAT Thematic Mapper sensors. Int. J. Remote Sens. 6 (5), 697–716. https://doi.org/10.1080/ 01431168508948492. Matthew, M.W., Adler-Golden, S.M., Berk, A., Richtsmeier, S.C., Levine, R.Y., Bernstein, L.S., Acharya, P.K., et al., 2000. Status of atmospheric correction using a MODTRAN4-based algorithm. In: Shen, S.S., Descour, M.R. (Eds.), Algorithms for Multispectral, Hyperspectral, and Ultraspectral Imagery Vi. Spie-Int Soc Optical Engineering, Bellingham, pp. 199–207. McGhie, S., 2005. Salinity indicator plants. In: Ryan, Melissa, McGhie, Sian (Eds.), Planning New South Wales. Department of Infrastructure and Resources Natural, Local Government Salinity Initiative ; no. 8. Dept. of Infrastructure, Planning and Natural Resources, Sydney. Metternicht, G., 1998. Analysing the relationship between ground-based reflectance and environmental indicators of salinity processes in the Cochabamba valleys (Bolivia). Int. J. Ecol. Environ. Sci. 24 (4), 359–370.

Declaration of Competing Interest The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. Acknowledgment The authors extend sincere thanks and appreciation to the Hrvatske vode (Croatian Waters) for kindly providing the in situ data and salinity monitoring reports for the river Neretva delta. Funding This research was partially supported through project KK.01.1.1.02.0027, a project co-financed by the Croatian Government and the European Union through the European Regional Development Fund–the Competitiveness and Cohesion Operational Programme. Appendix A. Supplementary data Supplementary data to this article can be found online at https:// doi.org/10.1016/j.ecolind.2019.105924. References Abbas, A., Khan, S., 2007. Using remote sensing techniques for appraisal of irrigated soil salinity. Paper presented at the Modsim 2007: International Congress on Modelling and Simulation: Land, Water and Environmental Management: Integrated Systems for Sustainability, Christchurch. Abbas, A., Khan, S., Hussain, N., Hanjra, M.A., Akbar, S., 2013. Characterizing soil salinity in irrigated agriculture using a remote sensing approach. Phys. Chem. Earth, Parts A/B/C 55–57, 43–52. https://doi.org/10.1016/j.pce.2010.12.004. Abrol, I.P., Yadav, J.S.P., Massoud, F.I., 1988. Salt-affected soils and their management. FAO Soils Bull. 39, 143. Allbed, A., Kumar, L., 2013. Soil salinity mapping and monitoring in arid and semi-arid regions using remote sensing technology: a review. Adv. Remote Sens. 2 (04), 373. Ayers, R.S., Westcot, D.W., 1985. Water Quality for Agriculture. Food and Agriculture Organization of the United Nations. Bannari, A., Guedon, A.M., El-Harti, A., Cherkaoui, F.Z., El-Ghmari, A., 2008. Characterization of slightly and moderately saline and sodic soils in irrigated agricultural land using simulated data of advanced land imaging (EO-1) sensor. Commun. Soil Sci. Plant Anal. 39 (19–20), 2795–2811. https://doi.org/10.1080/ 00103620802432717. Bastiaanssen, W.G.M., Molden, D.J., Makin, I.W., 2000. Remote sensing for irrigated agriculture: examples from research and possible applications. Agric. Water Manag. 46, 137–155. https://doi.org/10.1016/S0378-3774(00)00080-9. Ben-Dor, E., 2002. Quantitative remote sensing of soil properties. In: Sparks, D.L. (Ed.), Advances in Agronomy. Elsevier Academic Press Inc, San Diego, pp. 173–243. Birth, G.S., McVey, G.R., 1968. Measuring the color of growing turf with a reflectance spectrophotometer1. Agron. J. 60 (6), 640. https://doi.org/10.2134/agronj1968. 00021962006000060016x. Ceccato, P., Gobron, N., Flasse, S., Pinty, B., Tarantola, S., 2002. Designing a spectral index to estimate vegetation water content from remote sensing data: Part 1. Remote Sens. Environ. 82 (2–3), 188–197. https://doi.org/10.1016/s0034-4257(02)00037-8. Cenfetelli, Ronald T., Bassellier, Geneviève, 2009. Interpretation of formative measurement in information systems research. MIS Q. 33 (4), 689–707. https://doi.org/10. 2307/20650323. Chatfield, C., 2004. The Analysis of Time Series: an Introduction, 6th ed. CRC Press, Florida, US. Daliakopoulos, I.N., Tsanis, I.K., Koutroulis, A., Kourgialas, N.N., Varouchakis, A.E., Karatzas, G.P., Ritsema, C.J., 2016. The threat of soil salinity: a European scale review. Sci. Total Environ. 573, 727–739. https://doi.org/10.1016/j.scitotenv.2016. 08.177. Dehaan, R.L., Taylor, G.R., 2002. Field-derived spectra of salinized soils and vegetation as indicators of irrigation-induced soil salinization. Remote Sens. Environ. 80 (3), 406–417. https://doi.org/10.1016/s0034-4257(01)00321-2.

13

Ecological Indicators 110 (2020) 105924

I. Racetin, et al.

Romić, D., F. Tomić, J. Marušić, Ž. Klačić, and M. Vranješ. 1999. Kakvoća i mogućnost upotrebe vode za navodnjavanje u dolini Neretve. In, edited by Dragutin Gereš and Jure Margeta, 813-9. Dubrovnik, Croatia: Hrvatske vode. Rouse Jr, J.W., Haas, R.H., Schell, J.A., Deering, D.W., 1974. Monitoring Vegetation Systems in the Great Plains with Erts. NASA Special Publication 309-. Scudiero, E., Corwin, D.L., Anderson, R.G., Skaggs, T.H., 2016. Moving forward on remote sensing of soil salinity at regional scale. Front. Environ. Sci. 4, 1–5. https://doi.org/ 10.3389/fenvs.2016.00065. Scudiero, E., Corwin, D.L., Anderson, R.G., Yemoto, K., Clary, W., Wang, Z., Skaggs, T.H., 2017. Remote sensing is a viable tool for mapping soil salinity in agricultural lands. Calif. Agric. 71 (4), 231–238. https://doi.org/10.3733/ca.2017a0009. Scudiero, E., Skaggs, T.H., Corwin, D.L., 2014. Regional scale soil salinity evaluation using Landsat 7, western San Joaquin Valley, California, USA. Geoderma Regional 23 (Supplement C) 82–90. https://doi.org/10.1016/j.geodrs.2014.10.004. Stolte, J., M. Tesfai, J. Keizer, L. Øygarden, S. Kværnø, F. Verheijen, P. Panagos, C. Ballabio, and R. Hessel. 2016. Soil threats in Europe; EUR 27607 EN. Szabolcs, I., 1979. Review of Research on Salt-Affected Soils. Unesco. Tamas, J., Lenart, C., 2006. Analysis of a small agricultural watershed using remote sensing techniques. Int. J. Remote Sens. 27 (17), 3727–3738. https://doi.org/10. 1080/01431160500500474. Van Der Werff, H., Van Der Meijde, M., Jansma, F., Van Der Meer, F., Groothuis, G.J., 2008. A spatial-spectral approach for visualization of vegetation stress resulting from pipeline leakage. Sensors 8, 3733–3743. https://doi.org/10.3390/s8063733. Vermeulen, D., van Niekerk, A., 2016. Evaluation of a WorldView-2 image for soil salinity monitoring in a moderately affected irrigated area. J. Appl. Remote Sens. 10, 026025. https://doi.org/10.1117/1.JRS.10.026025. Zhang, H.J., Dong, H.Z., Li, W.J., Zhang, D.M., 2012. Effects of soil salinity and plant density on yield and leaf senescence of field-grown cotton. J. Agron. Crop Sci. 198, 27–37. https://doi.org/10.1111/j.1439-037X.2011.00481.x. Zhang, L., Zhou, Z., Zhang, G., Meng, Y., Chen, B., Wang, Y., 2017. Monitoring cotton (Gossypium hirsutum L.) leaf ion content and leaf water content in saline soil with hyperspectral reflectance. Eur. J. Remote Sens. 47 (1), 593–610. https://doi.org/10. 5721/EuJRS20144733. Zhang, T.T., Zeng, S.L., Gao, Y., Ouyang, Z.T., Li, B., Fang, C.M., Zhao, B., 2011. Using hyperspectral vegetation indices as a proxy to monitor soil salinity. Ecol. Ind. 11 (6), 1552–1562. https://doi.org/10.1016/j.ecolind.2011.03.025. Zovko, M., Romić, D., Colombo, C., Di Iorio, E., Romić, M., Buttafuoco, G., Castrignanò, A., 2018. A geostatistical Vis-NIR spectroscopy index to assess the incipient soil salinization in the Neretva River valley, Croatia. Geoderma 332, 60–72. https://doi. org/10.1016/j.geoderma.2018.07.005. Zovko, M., Romić, D., Romić, M., Ondrašek, G., 2013. Soil and water management for sustained agriculture in alluvial plains and flood plains exposed to salinity: a case of Neretva River Valley. In: Ahmad, Parvaiz, Azooz, M.M., Prasad, M.N.V. (Eds.), Ecophysiology and Responses of Plants under Salt Stress. Springer New York, New York, NY, pp. 473–494.

Metternicht, G., Zinck, A., 2003. Remote sensing of soil salinity: potentials and constraints. Remote Sens. Environ. 85 (1), 1–20. https://doi.org/10.1016/s00344257(02)00188-8. Metternicht, G., Zinck, A., 2008. Soil salinity and salinization hazard. In: Metternicht, Graciela, Zinck, Alfred (Eds.), Remote Sensing of Soil Salinization: Impact on Land Management. CRC Press, Boca Raton, pp. 1–20. Mougenot, B., Pouget, M., Epema, G.F., 1993. Remote sensing of salt affected soils. Remote Sens. Reviews 7 (3–4), 241–259. https://doi.org/10.1080/ 02757259309532180. Mulla, D.J., 2013. Twenty five years of remote sensing in precision agriculture: key advances and remaining knowledge gaps. Biosyst. Eng. 114 (4), 358–371. https://doi. org/10.1016/j.biosystemseng.2012.08.009. Muller, S.J., van Niekerk, A., 2016a. An evaluation of supervised classifiers for indirectly detecting salt-affected areas at irrigation scheme level. Int. J. Appl. Earth Obs. Geoinf. 49, 138–150. https://doi.org/10.1016/j.jag.2016.02.005. Muller, S.J., van Niekerk, A., 2016b. Identification of WorldView-2 spectral and spatial factors in detecting salt accumulation in cultivated fields. Geoderma 273, 1–11. https://doi.org/10.1016/j.geoderma.2016.02.028. Neter, J., Kutner, M.H., Nachtsheim, C.J., Wasserman, W., 1996. Applied Linear Statistical Models. Irwin, Chicago. Pitman, M.G., Läuchli, A., 2002. Global impact of salinity and agricultural ecosystems. In: Läuchli, André, Lüttge, Ulrich (Eds.), Salinity: Environment – Plants – Molecules. Springer Netherlands, Dordrecht, pp. 3–20. Procjena rizika od katastrofa za Republiku Hrvatsku [Disaster Risk Assessment in Croatia].“ In 463 2015 Republic of Croatia National Protection and Rescue Directorate. Rahmati, M., Hamzehpour, N., 2017. Quantitative remote sensing of soil electrical conductivity using ETM plus and ground measured data. Int. J. Remote Sens. 38 (1), 123–140. https://doi.org/10.1080/01431161.2016.1259681. Rasool, S., Hameed, A., Azooz, M.M., Muneebu, R., Siddiqi, T.O., Ahmad, P., 2013. Salt stress: causes, types and responses of plants. In: Ahmad, Parvaiz, Azooz, M.M., Prasad, M.N.V. (Eds.), Ecophysiology and Responses of Plants under Salt Stress. Springer New York, New York, NY, pp. 1–24. Rhoades, J.D., 1996. Salinity: electrical conductivity and total dissolved solids. In: Sparks, D.L., Page, A.L., Helmke, P.A., Loeppert, R.H. (Eds.), Methods of Soil Analysis Part 3—Chemical Methods. Soil Science Society of America, American Society of Agronomy, Madison, WI, pp. 417–435. Rhoades, J.D., Kandiah, A., Mashali, A.M., 1992. The Use of Saline Waters for Crop Production. FAO, Rome, Italy. Romic, D., Romic, M., Zovko, M., Bakic, H., Ondrasek, G., 2012. Trace metals in the coastal soils developed from estuarine floodplain sediments in the Croatian Mediterranean region. Environ. Geochem. Health 34 (4), 399–416. https://doi.org/ 10.1007/s10653-012-9449-z. Romić, D., M. Romić, M. Zovko, M. Bubalo, G. Ondrašek, S. Husnjak, N. Maurović, H. Bakić, L. Matijević, and M. Vranješ. 2014. Monitoring of soil and water salinity in the Neretva River valley-five-year (2009-2013) report.

14