Surface Water Change Detection via Water Indices and Predictive Modeling Using Remote Sensing Imagery: A Case Study of Nuntasi-Tuzla Lake, Romania

Water body feature extraction using a remote sensing technique represents an important tool in the investigation of water resources and hydrological drought assessment. Nuntasi-Tuzla Lake, a component of the Danube Delta Natural Reserve, is located on the Romanian Black Sea littoral. On account of an event in summer 2020, when the lake surface water decreased significantly, this study aims to identify the variation of the Nuntasi-Tuzla Lake surface water over a long-term period in correlation with human intervention and climate change. To this end, it provides an analysis in the period 1965–2021 via hydrological drought indices and data mining classification. The latter approach is based on several water indices derived from Landsat TM/ETM+/OLI and MODIS full-time series datasets: Normalized Difference Vegetation Index (NDVI), Normalized Difference Vegetation Index (NDVI), Modified NDWI (MNDWI), Weighted Normalized Difference Water Index (WNDWI), and Water Ratio Index (WRI). The experimental results indicate that the proposed classification methods can extract relevant features from waterbodies using remote sensing imagery with a high accuracy. Moreover, the study shows a similarity in the evolution of surface water cover identified with the data mining classification and the drought periods detected in the flow data series for the Nuntasi and Sacele Rivers that supply the Nuntasi-Tuzla Lake. Overall, the results of our investigation show that human intervention and hydrological drought had an extensive impact on the long-term changes in surface water of the Nuntasi-Tuzla Lake.


Introduction
Climatic changes can induce modification of lake surfaces and riverbeds with severe implications for agricultural and economic activities [1][2][3]. The scientific literature demonstrates the "sensitivity" of physical, chemical, and biological lakes' properties to climate change and human intervention [4][5][6]. Since 1960 permanent surface water has been shrinking, drying or has completely disappeared: Aral Lake, Urmia Lake (70% of its surface between 2012-2017 [7]), Chad Lake, Lop Nur, etc. [4,5]. Lake surface area was used by Benson and Paillet [8] in order to establish an indicator of lakes' hydrological response to climatic change. In comparison with classic methods, Remote Sensing (RS) technology presents a viable alternative to improve the observation of lake water surface in hydrological studies.
Initially, the NDVI (Normalized Difference Vegetation Index) was an indicator used to efficiently monitor vegetation condition [9] and drought activity [10,11]. This index was used also to extract water features [12]. McFetters [13] developed a new method to delineate water bodies based on the spectral characteristics of water which refers to the capability of a water body to absorb near infrared radiation (NIR) and allows visible green light to penetrate the water body. This new indicator, NDWI (Normalized Difference Water Figure 1. Map of the study area, the Nuntasi-Tuzla Lake. In this context, the purpose of this paper is to: (i) analyze the hydrological data in order to understand the response of surface water to climate change or/and human intervention and (ii) evaluate the performance of the most commonly used water indices In this context, the purpose of this paper is to: (i) analyze the hydrological data in order to understand the response of surface water to climate change or/and human intervention and (ii) evaluate the performance of the most commonly used water indices using time series Landsat data in order to capture the variation in surface water during 1984-2021 period. These results will contribute to a better understanding of the lake's evolution in the context of global changes.

Study Area
The Nuntasi-Tuzla Lake is situated in the well-known Histria region (Histria or Istros acropolis is the first urban settlement, founded by ancient Milesian in the 7th century BC). From a geomorphological point of view, the Histria region is a coastal lowland developed at the contact between the Central Dobrogea plateau and the Black Sea [19]. The average altitude is 25 m with a slope that decreases towards the sea. Geologically, green schist covered with a loess layer of varying thickness (2-15 m) represents the basement. The climate of the region is temperate continental with marine influences (precipitation of 400 mm and 11 • C average temperature). From a hydrological point of view, the lake is part of the Razim-Sinoe lagunar complex, a component of the Danube Delta Natural Reserve. Today, the Razim-Sinoe complex (Figure 2), initially a marine bay (the former Halmyris gulf), is composed of the following lakes: Razim, Sinoe, Golovita, Zmeica, Babadag, Nuntasi-Tuzla, Istria, and Ceamurlia si Agighiol. The average altitude is 25 m with a slope that decreases towards the sea. Geologically, green schist covered with a loess layer of varying thickness (2-15 m) represents the basement. The climate of the region is temperate continental with marine influences (precipitation of 400 mm and 11 °C average temperature). From a hydrological point of view, the lake is part of the Razim-Sinoe lagunar complex, a component of the Danube Delta Natural Reserve. Today, the Razim-Sinoe complex (Figure 2), initially a marine bay (the former Halmyris gulf), is composed of the following lakes: Razim, Sinoe, Golovita, Zmeica, Babadag, Nuntasi-Tuzla, Istria, and Ceamurlia si Agighiol. Concerning the lake's genesis, the latest scientific arguments clearly demonstrate that neo-tectonic activity was the main factor that led to the Nuntasi-Tuzla and Istria lakes' development 1200-1600 years ago [19]. The other lakes were separated from the Black Sea by sand-belts [17], which formerly had outlets allowing the penetration of sea water. The lakes were interconnected by small creeks.  Concerning the lake's genesis, the latest scientific arguments clearly demonstrate that neo-tectonic activity was the main factor that led to the Nuntasi-Tuzla and Istria lakes' development 1200-1600 years ago [19]. The other lakes were separated from the Black Sea by sand-belts [17], which formerly had outlets allowing the penetration of sea water. The lakes were interconnected by small creeks.
According to Breier [20], Nuntasi-Tuzla Lake has the following morphometric features: the length is 2 km, the breadth varies between 1.7 and a maximum of 3 km, the depth of the lake varies between 2.15 m and 6.15 m, and the surface area is 1050 hectares. The Nuntasi-Tuzla Lake is supplied by two rivers: Nuntasi and Sacele. From the economic point of view, according to town hall documents, agriculture and fisheries are the predominant economic sectors in the region. The mud baths were a great tourist attraction at the beginning of the 20th century due to its therapeutic qualities. Nowadays, the most important touristic attraction in the region is Histria ancient city.
In a recent work [21], we described the irrigation system built in 1971-1975 and operable till 1990. After this year, the demand for irrigation decreased to 20% of its capacity.
In order to ensure water demand, the Razim-Sinoe lagoons were transformed into a large freshwater reservoir. The three greatest hydraulic works were constructed in that period: (i) closure of the natural gateway between the lagoons and Black Sea, (ii) closurevia sluices of the connections between the lakes and (iii) resizing of three canals which supply the reservoir with fresh water from the Danube (St George Arm) so as to ensure a total discharge of 260 mc/s (for a probability of occurrence of 1%) with an average of 80 mc/s [21]. These works led to important changes in the hydrological regime of the region and to a great disturbance of the ecological balance of all the lakes. Nuntasi-Tuzla Lake was connected to Istria Lake by a canal and sluice. This canal was silted up with alluvium.
A detailed description of all human intervention in this area is provided in Section 3.1.

Hydrological Datasets
Two hydrological parameters are monitored at the two hydrometric stations (Nuntasi and Sacele). This study is based on two time data series: annual average discharge for Nuntasi and Sacele spanning the 1965-2020 period, and daily discharge for the period 2008-2020. The annual average discharge was used in order to investigate changes in the time series. The time series data were collected by the National Administration "Romanian Water" Dobrogea-Littoral Branch at the hydrometric station situated in the Nuntasi-Tuzla lake basin.

Remote Sensing Datasets
The time-series datasets of Landsat surface reflectance are acquired directly from the GEE (Google Earth Engine) platform for the period 1984-2021 in the following way: (i) from 1984 to 1999, and from 2003 to 2011, the data were collected from Landsat 4-5 Thematic Mapper (TM); (ii) for 2000, 2001 and 2002, Landsat 7 Enhanced Thematic Mapper (ETM+) images were collected; and (iii) for 2013-2021, the data were collected from Landsat 8 Operational Land Imager (OLI). The datasets belong to Landsat Collection 1-Level 1 corrected data, which have the highest radiometric and positional quality and are suitable for time-series analysis. The spatial resolution of all Landsat scenes is 30 m.
Due to Landsat 7 ETM+ data gaps from 2012, the time-series dataset of the moderateresolution imaging spectroradiometer (MODIS) TERRA surface reflectance (MOD09GA) was used to determine the surface water for the seasons of 2012. MODIS images with 500 m spatial resolution were acquired from the GEE platform. The study area is entirely contained within path 181 and row 29 for all Landsat TM/ETM+/OLI and Modis images ( Figure 1).

Methods
The methodology used in this paper consists of three parts: A. a literature investigation to understand human intervention in the Nuntasi-Tuzla Lake basin; B. a hydrological anal-ysis of two tributaries (Nuntasi and Sacele) concerning hydrological drought occurrence, which refers to (i) detecting changes in time series data through KhronoStat software and (ii) calculating hydrological drought; and C. RS analysis, which refers to: (i) processing of the Landsat and MODIS datasets so that only the images with Cloud Cover of less than 30% are investigated (ii) calculation of RS indices (RSI) NDVI, NDWI, MNDWI, WNDWI, WRI and their average for each season of each year of the period considered; (iii) in order to extract the water surface features, a decision tree was built. We chose the Classification and Regression Tree (CART) model developed by Breiman et al. [22], as a very popular statistical learning tool for analyzing large complex data.

Hydrological Analysis
KhronoStat (HydroSciences, Montpellier, France) is statistical software able to identify the sudden changes in a hydro or climatic data series using a set of statistical tests: the Pettitt rank-based test (non-parametric test), Buishand "U" test, Lee and Heghinian Bayesian method (parametric test) and Hubert procedure segmentation. The Buishand "U" and Lee and Heghinian tests are applicable if the time series investigated are normal. The Hubert procedure is the only one able to detect multiple breaks. The series is divided into "m" segments and, in order to limit the segmentation, the averages of two contiguous segments must be significantly different. This constraint is satisfied by Sheffe's test [23,24].
Considering hydrological drought analysis, an overview of the principal scientific literature [25][26][27][28] demonstrates that methodologies to characterize hydrological droughts can be divided into two categories (   [22], as a very popular statistical learning tool for analyzing large complex data.

Hydrological Analysis
KhronoStat (HydroSciences, Montpellier, France) is statistical software able to identify the sudden changes in a hydro or climatic data series using a set of statistical tests: the Pettitt rank-based test (non-parametric test), Buishand "U" test, Lee and Heghinian Bayesian method (parametric test) and Hubert procedure segmentation. The Buishand "U" and Lee and Heghinian tests are applicable if the time series investigated are normal. The Hubert procedure is the only one able to detect multiple breaks. The series is divided into "m" segments and, in order to limit the segmentation, the averages of two contiguous segments must be significantly different. This constraint is satisfied by Sheffe's test [23,24].
Considering hydrological drought analysis, an overview of the principal scientific literature [25][26][27][28] demonstrates that methodologies to characterize hydrological droughts can be divided into two categories (  Based on the SPA method, two droughts could be pooled if the reservoir has not totally "recovered" after the first drought event. Our interest is to find minor droughts also and not only major ones. For this study, we decided to use the TLM method. The TLM method, well known as the "method-of-crossing theory", was introduced by Yevjevich in 1967 [29]. First, a threshold level, Qo, is chosen. The deficit begins where the flow values are below this Qo level and ends where the flow is above Qo value. We propose to use as threshold a value in the range of Q70%-Q95% resulting from the FDC analysis. To conclude, in order to evaluate the hydrological drought in the Nuntasi and Sacele Rivers, the methodology proposed is: (i) in order to select the driest period, an Based on the SPA method, two droughts could be pooled if the reservoir has not totally "recovered" after the first drought event. Our interest is to find minor droughts also and not only major ones. For this study, we decided to use the TLM method. The TLM method, well known as the "method-of-crossing theory", was introduced by Yevjevich in 1967 [29]. First, a threshold level, Qo, is chosen. The deficit begins where the flow values are below this Qo level and ends where the flow is above Qo value. We propose to use as threshold a value in the range of Q70%-Q95% resulting from the FDC analysis. To conclude, in order to evaluate the hydrological drought in the Nuntasi and Sacele Rivers, the methodology proposed is: (i) in order to select the driest period, an analysis of annual discharge time series data Water 2022, 14, 556 6 of 20 (1965-2020) will be performed with KhronoStat software which detects a sudden change in time series data; (ii) an assessment threshold value with the FDC curve; (iii) hydrological drought assessment with TLM methods using daily discharge data.

RS Analysis
Using the special spectral reflectance (RS) properties of water, it is possible to differentiate between water and other surface materials. A multiple RS indices algorithm based on supervised classification was used to identify the surface water cover of Nuntasi Tuzla Lake. The spectral indices measurements acquired from Landsat TM/ETM+/OLI and MODIS datasets from 1984-2021 were analyzed by the algorithm that trains the randomly chosen samplings to construct a binary decision tree, and finally decodes this in order to classify the surface water ( Figure 4). The algorithm was fully implemented and run on the GEE platform. analysis of annual discharge time series data (1965-2020) will be performed with KhronoStat software which detects a sudden change in time series data; (ii) an assessment threshold value with the FDC curve; (iii) hydrological drought assessment with TLM methods using daily discharge data.

RS Analysis
Using the special spectral reflectance (RS) properties of water, it is possible to differentiate between water and other surface materials. A multiple RS indices algorithm based on supervised classification was used to identify the surface water cover of Nuntasi Tuzla Lake. The spectral indices measurements acquired from Landsat TM/ETM+/OLI and MODIS datasets from 1984-2021 were analyzed by the algorithm that trains the randomly chosen samplings to construct a binary decision tree, and finally decodes this in order to classify the surface water ( Figure 4). The algorithm was fully implemented and run on the GEE platform. Cloud filtering and masking of remote sensing data is a necessity to monitor the land surface cover more accurately. We used the cloud flag from the Level 1 Quality Assessment (QA) band of the Landsat datasets to limit the cloud cover at less than 30% in each selected image. Then, cloudy pixels were masked and not taken into consideration during the modeling of the surface water. This method is considered robust against errors [30,31]. After applying the cloud filter, we selected 233 Landsat RS images covering three seasons (spring, summer and autumn) of each year for the period investigated (Table 1). Similarly, the MOD09GA dataset for each season of 2012 was selected. In winter, most of the RS images did not pass the cloud filter, thus we decided to omit this season from our research. Cloud filtering and masking of remote sensing data is a necessity to monitor the land surface cover more accurately. We used the cloud flag from the Level 1 Quality Assessment (QA) band of the Landsat datasets to limit the cloud cover at less than 30% in each selected image. Then, cloudy pixels were masked and not taken into consideration during the modeling of the surface water. This method is considered robust against errors [30,31]. After applying the cloud filter, we selected 233 Landsat RS images covering three seasons (spring, summer and autumn) of each year for the period investigated (Table 1). Similarly, the MOD09GA dataset for each season of 2012 was selected. In winter, most of the RS images did not pass the cloud filter, thus we decided to omit this season from our research. seasonal RSI in summer 2020 is taken from the RSI indices summed from 1 June to 31 August 2020 divided by the valid data (number of RS images) from the same period.
Water is >1 1 weighted coefficient; if a = 0 than WNDWI = MNDWI, if a = 1 than WNDWI = NDWI. For our study, the weighted coefficient a was set to 0.50, as the tests show that this value allows a high overall accuracy of the index performance [17].
In the formula, Green denotes the reflectance of the green band (Band 3 of the Landsat OLI data, Band 2 of the Landsat TM/ETM+ imagery, Band 4 of the MODIS data), Red represents the reflectance of the red band (Band 4 of the Landsat OLI data, Band 3 of the Landsat TM/ETM+ imagery, Band 1 of the MODIS data), NIR designates the reflectance of the near-infrared band (Band 5 for Landsat OLI, Band 4 for Landsat TM/ETM+, Band 2 of the MODIS data), MIR indicates the reflectance of the middle-infrared band (Band 6 for Landsat OLI and MODIS, Band 5 for Landsat TM/ETM) and SWIR stands for the reflectance of the SWIR1 band, which corresponds to Band 6 for Landsat OLI and MODIS data, and Band 5 for Landsat TM/ETM+ datasets.
The CART models lie on the supervised branch of the machine learning (ML) algorithms and may be used for both classification and regression problems. They find homogeneous subsets by recursively partitioning the input data, based on independent variable splitting criteria employing variance-minimizing algorithms. Being a non-parametric algorithm, intuitive and with a significant ability to characterize complex interactions among variables, the CART models have been widely used in different practical remote sensing applications [32,33]. For this study, we employed the CART algorithm provided by the GEE, with the seasonal averaged values of the RS indices (NDVI, NDWI, MNDWI, WNDWI and WRI) for each year as the input data. The surface water results for each season of the period investigated were then predicted based on the trained classification model.
To evaluate the performance of the RSI indices and the trained classification algorithm, three Landsat images for each period studied were randomly selected (Table 3). In order to represent the water and the land features in the study area, we used the stratified random sampling method to select 240 test samples (120 water and 120 non-water samples). The water samples consist of river, lake and sea water pixels, and the non-water samples are composed of soil, vegetation and built-up area pixels. The validation of the accuracy assessment has been achieved by referring to the original satellite images and Google Earth images.

Human Intervention
A review of the principal scientific literature [27,[34][35][36][37][38][39][40] concerning human intervention in the Nuntasi-Tuzla Lake region demonstrates huge landscape modification that began more than a century ago. Two series of works were carried out in the region: the first began in the period 1903-1905 and the second started in 1970. The first series was carried out for fisheries and finalized in four stages [37]: (i) 1903-1916, when a link between the Sf Gheorghe arm and the Razim-Sinoe Lake was made in order to improve the fresh water circulation [41]; work began in 1905 with the Dunavat canal and was continued with the Dranov (1912) and Crasnicol canals; (ii) in 1930-1948 and (iii) 1950-1963, a great number of canals were excavated in order to connect the Danube with the lakes (Lipoveni, Fundea, Mustaca canals) and these lakes with each other; (iv) in 1963-1985, the existing canals were broadened and deepened [37] and a water gate at Portita mouth was built [39]. In 1963, a study was conducted concerning the territorial planning of the Dranov-Razelm-Sinoie Lacuster Complex [40]. This document recommended the closure of the Razim Sinoe Lagoon complex and the transformation of this system into an irrigation reservoir. In this context, in the period 1971-1975 the Portita mouth was closed by a dam [28] and a series of hydraulic works were made to control the communication between Razim Lake and the Black Sea through Sinoe Lake. From 1971 to 1976, the irrigation system was built and consisted of six subsystems. The Sinoe irrigation subsystem, of which the study area is a part, came into operation in 1976. The Sinoe SPA (Supply Pumping Station) covered a surface of 57,162 ha and had a flow rate of 46.1 mc/s. The distribution network had a length of 157.3 km [36]. The irrigation system was in operation from 1976 to 1990. After this year, water demand for irrigation decreased. Under the new property law, a conservation process for the irrigation system began [41,42], resulting in their destruction (over 75% of these arrangements were inoperable in 2016 [42,43]. According with the National Agency for Land Improvement, in 2021 only 10,930 ha were contracted (19%) in the study area and only 1250 ha have been irrigated. All this human intervention modified the hydrological condition of the Nuntasi-Tuzla ecosystem. The former Research and Design Institute for Water Management (the current Aquaproiect Company) shows a high demineralization of lake water from 57.8 g/L (1934) to 1.6 g/L (1982), a decrease in production of therapeutic mud, which was 0 in 1980, and an increase in the fresh water supply through the two tributaries, from 5 mil/mc in 1977 to 15 mil/mc in 1980 [44]. Braier [20] specified in her study that the salinity of this lake decreased due to the freshwater input and to the fact that the connection between this lake and the Istria lake was interrupted.

Hydrological Drought Analysis
Unfortunately, there are no meteorological or hydrological drought studies for this region. The idea to investigate this area appeared when, in 2020, the Istria town hall informed the authorities that Lake Nuntasi-Tuzla had dried up.
The Nuntasi Tuzla Lake is supplied by the Nuntasi and Sacele Rivers and we assume that there is a breakpoint in the time series data based on the variation of annual discharge presented in Figure 5, which shows a descending trend. As we showed in Section 2.2.1, in order to detect changes in the time series data, we used the annual discharge series for the 1965 to 2020 period. The multiannual average discharge for Nuntasi River is 0.348 m 3 /s, and is 0.081 m 3 /s for the Sacele River. It should also be specified that, for the Nuntasi River, the values for 2015 and 2017 are missing. The Nuntasi Tuzla Lake is supplied by the Nuntasi and Sacele Rivers and we assume that there is a breakpoint in the time series data based on the variation of annual discharge presented in Figure 5, which shows a descending trend. As we showed in Section 2.2.1, in order to detect changes in the time series data, we used the annual discharge series for the 1965 to 2020 period. The multiannual average discharge for Nuntasi River is 0.348 m 3 /s, and is 0.081 m 3 /s for the Sacele River. It should also be specified that, for the Nuntasi River, the values for 2015 and 2017 are missing. The hydrological data are not normally distributed, even after several transformations. For the Nuntasi and Sacele annual discharge series, Buishand and Lee and Heghinian tests were not performed. The Pettitt test rejected the null hypothesis H0 = no break for all confidence levels (90%, 95% and 99%) and admitted a break in 1997. The Hubert segmentation procedure detected four break points (1980, 1989, 1997 and 2006).  (Table 4), the average calculated for each subseries period detected by the Hubert procedure corresponded to the different stages of hydraulic work implemented in the study region and described in the previous Section 3.1. For the Nuntasi River, during the maximum operating period (1981-1989) of the irrigation system, the average annual discharge is 1.4 times greater than that from the previous period (1965)(1966)(1967)(1968)(1969)(1970)(1971)(1972)(1973)(1974)(1975)(1976)(1977)(1978)(1979)(1980). This increase is due to the infiltration resulting from the irrigation water. The drill installed in the region by the former Land Improvement Office provided an increasing in the aquifers level [36]. Moreover, after 1989 the multiannual discharge for the period 1990-1997 returned to the baseline (1965-1980) and continued to decrease, reaching a dangerous level (6.75 times less than the average annual discharge for the 1965-1980 period). The same behavior was seen in the Sacele River, but the decrease was only by 3.1 times. Based on this situation, in order to detect hydrological drought using the two procedure described The hydrological data are not normally distributed, even after several transformations. For the Nuntasi and Sacele annual discharge series, Buishand and Lee and Heghinian tests were not performed. The Pettitt test rejected the null hypothesis H 0 = no break for all confidence levels (90%, 95% and 99%) and admitted a break in 1997. The Hubert segmentation procedure detected four break points (1980, 1989, 1997 and 2006). For the Sacele discharge series, the Pettitt test rejected the null hypothesis for all confidence levels and admitted a break point in 2003, whilst the Hubert segmentation detected two break points, one in 1996 and one in 2003. For both Nuntasi and Sacele annual discharge, the breakpoints detected by the Pettitt and Hubert tests corresponded (1997 for the Nuntasi time series and 2003 for Sacele, respectively). As shown in the following table (Table 4), the average calculated for each subseries period detected by the Hubert procedure corresponded to the different stages of hydraulic work implemented in the study region and described in the previous Section 3.1. For the Nuntasi River, during the maximum operating period (1981-1989) of the irrigation system, the average annual discharge is 1.4 times greater than that from the previous period (1965)(1966)(1967)(1968)(1969)(1970)(1971)(1972)(1973)(1974)(1975)(1976)(1977)(1978)(1979)(1980). This increase is due to the infiltration resulting from the irrigation water. The drill installed in the region by the former Land Improvement Office provided an increasing in the aquifers level [36]. Moreover, after 1989 the multiannual discharge for the period 1990-1997 returned to the baseline (1965)(1966)(1967)(1968)(1969)(1970)(1971)(1972)(1973)(1974)(1975)(1976)(1977)(1978)(1979)(1980) and continued to decrease, reaching a dangerous level (6.75 times less than the average annual discharge for the 1965-1980 period). The same behavior was seen in the Sacele River, but the decrease was only by 3.1 times. Based on this situation, in order to detect hydrological drought using the two procedure described in the Section 2.3 (TLM and SPA, respectively) we used the daily discharge time series data for the 2007-2020 period. Firstly, we construct the FDC by plotting the empirical cumulative frequency (EFQ) of the river daily discharge against exceedance frequency. Figure 6 shows the FDC for the Nuntasi River. Firstly, we construct the FDC by plotting the empirical cumulative frequency (EFQ) of the river daily discharge against exceedance frequency. Figure 6 shows the FDC for the Nuntasi River. The difference between Q95%, Q90%, Q80% and Q75% is small (Table 5). In this context we chose as threshold (Qo) the value corresponding to Q90% for both rivers.  The difference between Q95%, Q90%, Q80% and Q75% is small (Table 5). In this context we chose as threshold (Qo) the value corresponding to Q90% for both rivers. For a threshold value of 0.03 (Q90%), the TLM identified not only major droughts but also minor droughts (Figure 7). There are 11 periods of hydrological drought in total (Table 6) and the driest years are 2013 and 2020 with 16% and 20% deficit, respectively, of the period investigated. For the Sacele river, only two hydrological drought periods have been detected, in 2019 (4 days) and 2020 (4 days).  (Table 6) and the driest years are 2013 and 2020 with 16% and 20% deficit, respectively, of the period investigated. For the Sacele river, only two hydrological drought periods have been detected, in 2019 (4 days) and 2020 (4 days).   Table 7 presents the confusion matrix of seasonal surface water detection results for the selected imagery introduced in Table 3. The overall accuracy (OA), Kappa coefficient, producer accuracy (PA) and user accuracy (UA) [45] of the surface water mapping are above 90%. On that account, the approach presented in this paper can be considered one of the optimal methods for surface water extraction for the study area (Figure 8).    Similarly, an accuracy verification procedure was applied on MODIS extracted surface water results ( Table 8). The confusion matrix provided accuracies of higher than 90%, proving that the trained classification algorithm could effectively detect surface water from the MODIS dataset. Further, in order to evaluate the effectiveness of the surface water detection results from the MODIS dataset for 2012, a correlation analysis was performed for the period 2000-2011 between the MODIS and Landsat results. The multiple correlation coefficient of surface water between Landsat and MODIS was 0.97 and the coefficient of determination (R-squared) was 0.95 (Figure 9), showing that the surface water results from the MODIS dataset closely correspond to those from Landsat. On that account, the surface water detected from MODIS can effectively be used for seasonal change analysis of the surface water cover in the study area.
2000-2011 between the MODIS and Landsat results. The multiple correlation coefficien of surface water between Landsat and MODIS was 0.97 and the coefficient of determina tion (R-squared) was 0.95 (Figure 9), showing that the surface water results from th MODIS dataset closely correspond to those from Landsat. On that account, the surfac water detected from MODIS can effectively be used for seasonal change analysis of th surface water cover in the study area.  Figure 10 shows the variation of the seasonal average surface water from 1984 t 2021 of Nuntasi-Tuzla Lake.  (Figure 11d). For the last period considered, 2013-2021, the largest seasonal average surface water was determined in spring 2019 (Figure 11e), and the smallest in autumn 2020 (Figure 11f). Figure 12a,b show the lake surface maps overlaid to produce lake surface area changes maps for the periods investigated.

Seasonal Water Surface Variation
The Hubert test detected five subseries in the seasonal water surface data series (in summer 1991, spring 2012, spring 2013, summer 2019 and summer 2020) which coincide with the detection in the flow data series of the Nuntasi and Sacele Rivers and full irrigation system operation.  (Figure 11d). For the last period considered, 2013-2021, the largest seasonal average surface water was determined in spring 2019 (Figure 11e), and the smallest in autumn 2020 (Figure 11f). Figure 12a,b show the lake surface maps overlaid to produce lake surface area changes maps for the periods investigated.
The Hubert test detected five subseries in the seasonal water surface data series (in summer 1991, spring 2012, spring 2013, summer 2019 and summer 2020) which coincide with the detection in the flow data series of the Nuntasi and Sacele Rivers and full irrigation system operation.
the lake surface maps overlaid to produce lake surface area changes maps for the periods investigated.
The Hubert test detected five subseries in the seasonal water surface data series (in summer 1991, spring 2012, spring 2013, summer 2019 and summer 2020) which coincide with the detection in the flow data series of the Nuntasi and Sacele Rivers and full irrigation system operation. The seasonal average surface water was approximately 992 ha from 1984 to summer 1991, and increased by 4%, on average, after 1991. It did not then suffer significant changes until spring 2012 ( Figure 10). Figure 12c shows the maximum growth of surface water area, 11.34%, from the smallest area in autumn 1989 to the largest in summer 2001. In summer 2012, the extracted surface water area decreased by 16.57% compared to the previous season, and it further reduced considerably till summer 2013 (Figure 11g), when the lake surface reached 670 ha, approximately. Figure 12d shows the decrease in surface water area from spring 2013 to summer 2013. Following 2014, the seasonal average surface water area increased steadily until spring 2019 (Figure 11e). Figure 12e shows the decrease in water area from spring 2019 to autumn 2019.
On average, the surface water from 2013 to spring 2020 decreased by about 9.26% compared to values from 2000 to 2012 and by 8.86% compared to the 1984-1999 period. In summer and autumn of 2020, it dramatically decreased by 18.76% and 57.83%, respectively, compared to the average surface water of previous years.
It is noted that the periods with the lowest values for lake surface coincide with the prolonged drought periods determined according to the TLM method (Table 6). The seasonal average surface water was approximately 992 ha from 1984 to summer 1991, and increased by 4%, on average, after 1991. It did not then suffer significant changes until spring 2012 ( Figure 10). Figure 12c shows the maximum growth of surface water area, 11.34%, from the smallest area in autumn 1989 to the largest in summer 2001. In summer 2012, the extracted surface water area decreased by 16.57% compared to the previous season, and it further reduced considerably till summer 2013 (Figure 11g), when the lake surface reached 670 ha, approximately. Figure 12d shows the decrease in surface water area from spring 2013 to summer 2013. Following 2014, the seasonal average surface water area increased steadily until spring 2019 (Figure 11e). Figure 12e shows the decrease in water area from spring 2019 to autumn 2019.
On average, the surface water from 2013 to spring 2020 decreased by about 9.26% compared to values from 2000 to 2012 and by 8.86% compared to the 1984-1999 period. In summer and autumn of 2020, it dramatically decreased by 18.76% and 57.83%, respectively, compared to the average surface water of previous years.
It is noted that the periods with the lowest values for lake surface coincide with the prolonged drought periods determined according to the TLM method (Table 6).
On average, the surface water from 2013 to spring 2020 decreased by about 9.26% compared to values from 2000 to 2012 and by 8.86% compared to the 1984-1999 period. In summer and autumn of 2020, it dramatically decreased by 18.76% and 57.83%, respectively, compared to the average surface water of previous years.
It is noted that the periods with the lowest values for lake surface coincide with the prolonged drought periods determined according to the TLM method (Table 6).

Conclusions
The main purpose of this study was to detect via remote sensing techniques the long-term variation of the Nuntasi-Tuzla Lake surface water during the 1984-2021 period in correlation with human intervention and climate change. We have reached several conclusions: In natural conditions, the lake system is connected to the Black Sea via the Portita, Periboina and Edighiol outlets ( Figure 2) and with the St George Arm of the Danube River via a system of canals and marshes. The nine lakes which compose the Razim Sinoe System are interconnected by the above mentioned system of canals. Human intervention has led to a deterioration of the ecosystem of the entire Razim Sinoe system, which has been isolated from the Black Sea, and the connections between the lakes have been cut by sluices. In this context, the water became freshwater. The cannal between Nuntasi-Tuzla Lake and Istria Lake has been silted since 1976. In

Conclusions
The main purpose of this study was to detect via remote sensing techniques the long-term variation of the Nuntasi-Tuzla Lake surface water during the 1984-2021 period in correlation with human intervention and climate change. We have reached several conclusions: (i) In natural conditions, the lake system is connected to the Black Sea via the Portita, Periboina and Edighiol outlets ( Figure 2) and with the St George Arm of the Danube River via a system of canals and marshes. The nine lakes which compose the Razim Sinoe System are interconnected by the above mentioned system of canals. Human intervention has led to a deterioration of the ecosystem of the entire Razim Sinoe system, which has been isolated from the Black Sea, and the connections between the lakes have been cut by sluices. In this context, the water became freshwater. The cannal between Nuntasi-Tuzla Lake and Istria Lake has been silted since 1976. In view of the fact that the Dobrogea region is the driest region of Romania [35] and in order to ensure agricultural development, the state authorities from the period 1968-1975 built an important irrigation system in this region. In the irrigation period, the two tributaries (Nuntasi and Sacele rivers) have supplied the lake constantly, but after 1990, when the irrigation was stopped, the river flow decreased over time reaching its lowest level in the 2004 (2007)-2020 period. In a recent publication, the authors [21] showed that "the precipitation increased starting with 2012 but the evapotranspiration losses are much larger than the precipitation increase". We can conclude that the budget is negative. It is apparent that the surface lake may be a subject of irreversible changes. (ii) The analysis of the daily flows of the two rivers during the 2007-2020 period detected several important drought events. Among these, two drought periods with long duration were determined, in 2013 and 2020. In this context, we further investigated if there was any influence of flow decreasing on the lake water surface and if there have been any other similar situations in the past. (iii) Using the CART model with the seasonal averaged values of RS indices (NDVI, NDWI, MNDWI, WNDWI and WRI) as the input data, we assessed the seasonal water lake surface variation during the period 1984-2021 with over 90% mapping accuracy, user accuracy and overall accuracy. The results of the proposed classification method revealed that the evolution of surface lake water is correlated with human intervention and the hydrological drought identified with the TLM method. Significant decrease during the 2003-2020 period was identified in the surface water lake's evolution, thus the hydrological drought identified in 2011, 2012, 2013 and 2020 corresponds with the lowest values for water lake surface. In our opinion, the method based on remote sensing data and the CART model is calibrated, due to the results obtained. Unfortunately, we have only a short series of daily records, which limits this study.
Overall, the findings of this study provide some insights into the evolution of the Nuntasi Tuzla Lake and the factors driving it through a long period of time, which may be further used to formulate scientific measures to prevent other drought disasters similar to that of August 2020. Future work may imply improvement of the surface water changes' detection method, by investigating other spectral indices and comparing with state-of-theart classification methods. Although the CART model is a very intuitive and powerful classifier, it often involves longer time to train and fails to meet the best performance if the problem has many uncorrelated variables. Other classification methods, such as support vector machines (SVM) or multilayer perceptron (MLP) neural network, which work well with large input data, can be adopted to identify surface water in Landsat satellite images.
According to the Romanian Water National Administration-Littoral Branch (ABADL), an extensive restoration campaign started in September 2020 and the canal between Istria and Nuntasi-Tuzla lakes was restored ( Figure 13).
A fragile equilibrium was reestablished and in spring of 2021 a group of flamingo birds was spotted.
An Integrated Drought Management Programme (IDMP) was launched in 2013 by GWP (Global Water Partnership) for Central and Eastern Europe. The main objective [21] is to support the responsible authorities to prepare and integrate the Drought Management plan (DMP) within the future Basin Management Plan (BMP). In this context, more efforts are needed to validate the method proposed in this study.
formance if the problem has many uncorrelated variables. Other classification methods, such as support vector machines (SVM) or multilayer perceptron (MLP) neural network, which work well with large input data, can be adopted to identify surface water in Landsat satellite images.
According to the Romanian Water National Administration-Littoral Branch (AB-ADL), an extensive restoration campaign started in September 2020 and the canal between Istria and Nuntasi-Tuzla lakes was restored ( Figure 13).