Assessment of Spatio-Temporal Changes in Water Surface Extents and Lake Surface Temperatures Using Google Earth Engine for Lakes Region, Türkiye

Abstract

This study aims to extract water surface area and lake surface water temperature (LSWT), and to present long-term spatio-temporal analysis of these variables together with meteorological parameters. Three lakes in Türkiye’s Lakes Region, namely, Lake Burdur, Egirdir, and Beysehir, were considered as test sites. The normalized difference water index (NDWI) was applied to Landsat 5 and 8 data from 2000 to 2021 to extract the water extent in the Google Earth Engine (GEE) cloud-based platform. In addition to the lake surface area, Landsat thermal images were used to examine the LSWT. The findings indicated that water pixels could be extracted rather accurately using NDWI, with an overall accuracy of 98%. Between 2000 and 2021, the water surface area value of Lake Burdur decreased by more than 22%, while Lake Egirdir has dropped by less than 4%, and Lake Beysehir has not changed noticeably. LSWT of Burdur and Egirdir Lakes increased by more than 2.13 °C and 0.32 °C, respectively, while it decreased about 1.5 °C for Beysehir Lake. The obtained results were evaluated with meteorological parameters and our findings indicated that human-induced activities were more dominant than climate effects over Lake Burdur, unlike the others.

1. Introduction

Lakes are among the most significant water resources, crucial to the ecosystem and regional climate because of their several roles, such as providing drinking water and regulating climate change, in addition to functioning as natural rehabilitation sites, fishing waters, reservoirs for groundwater, and their promotion of sustainable development. Every country and each nation are attempting to preserve water and to use the least amount of water as the lives of all living organisms depend upon it. Surface waters such as rivers, lakes, oceans, and wetlands are vital to the environment’s ecological balance and hydrological cycle [1,2]. Climate change and human activities may have a dramatic impact on inter-annual and intra-annual variations of surface waters, which may also have profound influence on human society and ecosystems [3,4]. The primary causes of inland water body deterioration and loss include land use changes, expansion in agricultural areas, water diversion by dams, sewerage infrastructure (particularly in river valleys and coastal regions), air and water pollution, and nutrient dumping [5]. Changes in water surface area caused by human activities are primarily related to water withdrawals for agricultural reasons and household purposes. Spatio-temporal mapping of lake surface water extent is critical for ecological and environmental impact assessment. Thus, obtaining reliable data on the dynamics of such water surfaces assists decision-makers in protecting them due to their numerous benefits, stated above, and enables them to develop long-term management strategies.

Remote sensing (RS) has been used as a cost-effective tool to monitor both spatial and temporal changes in surface water resources without direct contact with the object or ground [6,7]. Wide coverage, high return frequency, rich information content, and relatively low cost are the primary features of RS image data, which are utilized in land use and land cover mapping with diverse time-space spectra. Satellite images enable the detection of surface water bodies across a broad region at a certain moment. Since satellite data are readily available and time efficient, RS technology has become one of the most frequently utilized tools for identifying surface water bodies in recent years [1,8]. Satellite sensors, commonly used for water surface area mapping and monitoring, are The Advanced Very High-Resolution Radiometer (AVHRR) [9], Moderate Resolution Imaging Spectroradiometer (MODIS) [10,11], Systeme Probatoire d’Observation de la Terre (SPOT) [12], Visible Infrared Imaging Radiometer Suite Visible Infrared Imaging Radiometer (VIIRS) [13], Landsat [14,15,16], and Sentinel-2 [17,18] satellite images. One of the most used optical sensors in surface water resources and other environmental research is the Landsat satellite series, with its long-term data record and extensive return coverage. The sensors of all Landsat missions have provided one of the most useful, free, and continuous datasets with 16 days temporal resolution to the scientific community for long-term detecting and monitoring of the Earth’s surface since 1972 [19].

Since 2008, researchers have been using all accessible Landsat images to better obtain surface water change information thanks to the open sharing of Petabyte-Scale Landsat data [20]. Simultaneously, in recent years, several cloud computing platforms have been created to analyze large-scale geospatial data without requiring much technical knowledge or effort, and they have indicated a large amount of potential in large-scale land cover mapping. Cloud computing platforms such as Google Earth Engine (GEE) are typical cloud-based platforms that can handle data in an innately parallel, high-performance manner, and GEE comprises of multi-petabyte distant sensing data that has been preprocessed to be ready-to-use freely and effectively. It is also performed to investigate the variations in water extents of the lakes over time [21].

Water surface area detection is closely related to RS-based studies. Various techniques have been developed to extract the surface water body information because it is required to convert remotely sensed images into information for practical applications like water surface area mapping. Many studies used a combination of spectral water index methods in which the water surface areas can be extracted accurately, quickly, and easily. According to Buma et al. [22], the Modified Normalized Difference Water Index (MNDWI) was an effective indicator for mapping recent area changes in Lake Chad, with an overall accuracy of about 96%, and it was used to estimate the changes in the lake’s area from 2003 to 2016. Xia et al. [4] analyzed the data on the GEE platform using the available 16,760 scenes of Landsat images in the Huai River Basin from 1989 to 2017. Normalized difference vegetation index (NDVI), normalized difference water index (NDWI), and enhanced vegetation index (EVI) were used to quantify the spatio-temporal variability of surface water area changes over time. Wang et al. [23] used all available Landsat images (7534 scenes) in order to understand the long-term surface water area changes in the Hetao Plain by adopting the MNDWI, NDVI, and EVI to map the open water bodies from 1989 to 2019 in the GEE cloud platform. They studied the variability of water bodies over 32 years and identified the changing trend by evaluating the area and number of various water bodies. Bai et al. [24] examined the dynamic changes in three distinct types of inland lakes in Central Asia. In that study, Landsat images were employed for change detection due to their extensive archive of over 40 years, and the findings were also analyzed with respect to both climate change and human activities. In another study, the surface area of Lake Gregory in Australia was predicted using Landsat satellite data and a stochastic technique by Soltani et al. [25].

Apart from the surface water area information, surface water temperature is another important variable that affects the Earth’s energy budget, and many environmental processes. Water temperatures of lakes have been rising dramatically around the world, and this trend is likely to have a negative impact on aquatic ecosystems [26]. Therefore, several studies on lake’s water temperature had been conducted. Yamak and Comert [27] used Landsat images for the extraction of the Land Surface Temperature (LST) map, and low LST values were observed in vegetated areas, forest areas, and wetlands. In the study of Xie et al. [28], MODIS LST data were used to examine the spatial and temporal variations in lake surface water temperature (LSWT) of China’s lakes between 2001 and 2016. Their findings indicated that the 169 large lakes included in the study had an overall increasing trend in LSWT, with an average rate of 0.26 °C/decade.

Because of drought and misuse of water resources, half of Türkiye’s inland water bodies have lost their ecological and economic qualities for the past 40 years [29]. The main threats to inland water bodies in Türkiye, which are the primary causes of the loss, include unplanned water infrastructure projects such as dams built on rivers feeding lakes, in addition to inter-basin water transfers, excessive irrigation use in agriculture, pollution, illegal hunting and fishing, and other large-scale investments such as highways [30]. The RS data have been used to examine Türkiye’s inland water bodies, as well. Using Landsat imagery, Fethi et al. [31] calculated the surface area of Lake Beysehir as 652 km² in 1984, 617 km² in 2005, and 639 km² in 2014. Besides, for Lake Egirdir, the extent of the water surface was 463 km² in 1984, 455 km² in 1990, 458 km² in 2004, and 454 km² in 2013. Gozukara et al. [32] utilized Landsat imageries and determined that Lake Burdur’s water level decreased by 13.74 m (−2%) and the lake surface area decreased from 211 km² to 130 km² (−38.39%) over 42 years (1975–2017). Additionally, Dervisoglu [30] assessed the temporal changes in Türkiye’s inland Ramsar sites using Landsat satellite images and the NDWI. It was discovered that inland wetlands declined at varying rates between 1985 and 2020. Moreover, Abujayyab et al. [33] collected 570 Landsat scenes from GEE for the surface area extraction of Lake Burdur during a 34-year period between 1984 and 2019, where a seasonal fluctuation in the area of water bodies was determined by combining their suggested extraction technique with NDWI. Considering the study of Aslan and Koc-San [34], the water surface area and LSWT of Burdur and Egirdir lakes were analyzed using Landsat satellite images acquired in 1998, 2008, and 2018. Thus, after a 20-year period, it was discovered that the increase in LSWT values was about 2.2 °C for Lake Burdur and approximately 1.3 °C for Lake Egirdir, and both lakes’ water surface area decreased.

This study initially aimed to detect and analyze the spatio-temporal variations in water surface areas and LSWTs of Lake Burdur, Egirdir, and Beysehir in the Lakes Region using Landsat images through the GEE platform. These changes in the water surface area of the selected lakes were evaluated based on a period from 2000 to 2021. Many of the aforementioned studies utilized one satellite image for each analyzed year instead of a time-series dataset. Some of them analyzed only surface water extent, while the others only examined the LSWT. Furthermore, few of them used meteorological parameters to understand the effect of climate on the variations. In this study, we present detailed long-term analyses of the variations in water surface area and LSWT, and their relationship with meteorological parameters. Firstly, in order to understand the long-term changes, all available cloud-free Landsat images of Landsat TM and OLI sensors, which were acquired between the years 2000 and 2021, were determined to map the open water bodies. Secondly, the study employed the NDWI to detect the water pixels to calculate the lake water extent of each lake. Thirdly, we extracted LSWT images and examined the alleged connections between the lake water extents and the LSWT variable. Finally, a correlation analysis of the obtained results was conducted with meteorological parameters (i.e., precipitation, temperature, and evaporation) gathered from the meteorological satellites to reveal the climatic effects on the water surface area and LSWT. Our findings will serve as valuable outputs for other regional countries seeking to preserve water supplies.

2. Study Area

Türkiye’s 200 natural lakes have a combined surface area of about 906,000 hectares, representing over 1.6% of the country’s total land area administered by the General Directorate of Nature Conservation and National Parks (GDNCNP), which is under the Ministry of Agriculture and Forestry (MoAF). The national inventory includes 302 distinct lakes and wetlands with surface areas ranging from 8.0 to 357,269 hectares [35]. The selected study area is in the Lakes Region which is situated in Türkiye’s Mediterranean region between latitudes of 37°01′ N and 38°30′ N, and longitudes of 29°33′ E and 32°21′ E. It is composed of basins and 16 lakes, the majority of which are tectonic in origin. The western half of the area includes several closed basins, the most notable of which is Lake Burdur. The eastern region (Lake Beysehir and Lake Sugla basins) contains sub-basins of the greater Konya basin. The middle part of the region is linked to the Mediterranean Sea by a natural strait across Lake Kovada that allows surplus water from Lake Egirdir to flow into the Aksu River. In this region, the elevation varies from 1000 to 1500 m. The climate of this region is dry and hot in the summer, cold in the winter, and wet in the spring with an average yearly temperature between 11 °C and 14 °C at 1000 m [36]. In the meteorological records of 2021, the lowest temperature was observed to be 3.5 °C in Goksun, whereas the highest one was 46.1 °C in Koycegiz. Conversely, the mean annual precipitation was 507.5 mm in September 2021 while the long-term mean was 666.5 mm [37].

In this study, the three primary lakes in Lakes Region were selected as test sites in order to examine the changes in water surface areas during the period from 2000 to 2021 (Figure 1). Lake Burdur (brackish) is in the provinces of Burdur and Isparta. Lake Egirdir (freshwater) is in the province of Isparta, and Lake Beysehir (freshwater) is located in the provinces of Konya and Isparta. The maximum extent, elevation, maximum depth, and other related properties of the examined lakes are portrayed in

Table 1. Various properties of study lakes [38].

Lake	Maximum Extent (km2)	Type	Maximum Depth (m)	Elevation (m)	Purpose of Use	Protection Status
Burdur	195.37	Brackish	100	857	-	Ramsar Site
Egirdir	476.92	Freshwater	15	917	Irrigation	Natural Protected Area Drinking Water Reservoir
Beysehir	706.46	Freshwater	15	1124	Drinking water, Irrigation	Natural Park Wetland Protected Area Drinking Water Reservoir

.

3. Materials and Methods

3.1. Materials

In the study, two satellite sensors, including Landsat 5 TM and Landsat 8 OLI/TIRS, were used to cover as much of the study period as feasible from 2000 to 2021. While Landsat 5 TM provides coverage from 2000 to 2011, Landsat 8 OLI/TIRS provides coverage from 2013 to 2021. The characteristics of each sensor are provided in Table 2. The GEE platform was used to query and analyze all these images.

GEE is a cloud-based computing platform that allows the execution of geospatial analysis on Google’s cloud infrastructure. For creating and executing scripts, a web-based integrated development environment (IDE) is utilized. Additionally, this IDE is used to visualize geographical studies through the JavaScript Application Programming Interface (API). GEE libraries could be used to build JavaScript and Python programs [39]. The GEE allows for the filtering and removal of low-quality images using shadows, snow, and clouds. All images used in this research were masked based on the test sites’ borders. The time series continuity was improved by analyzing as many images as possible. Image collection from Landsat 5 TM and Landsat 8 OLI/TIRS satellites was acquired by applying a cloudiness ratio of less than 10% (no cloud over the lakes) over the scenes.

Table 2. Information about used satellite imageries [40].

Satellite	Spectral Resolution (µm)	Spatial Resolution (m)	Radiometric Resolution (bit)
Landsat 5 TM (2000–2011)	6 Optical Bands (0.45–2.35) 1 Thermal Band (10.40-12.50)	Blue, Green, Red, NIR, SWIR1-2: 30 m Thermal Infrared: 120 m	8
Landsat 8 OLI/TIRS (2013–2021)	9 Optical Bands (0.43–2.30) 2 Thermal Bands (10.6-12.51)	Cirrus, Coastal/Aerosol, Blue, Green, Red, NIR, SWIR1-2: 30 m Thermal Infrared 1–2: 100 m Panchromatic: 15m	16

Table 2. Information about used satellite imageries [40].

Satellite	Spectral Resolution (µm)	Spatial Resolution (m)	Radiometric Resolution (bit)
Landsat 5 TM (2000–2011)	6 Optical Bands (0.45–2.35) 1 Thermal Band (10.40-12.50)	Blue, Green, Red, NIR, SWIR1-2: 30 m Thermal Infrared: 120 m	8
Landsat 8 OLI/TIRS (2013–2021)	9 Optical Bands (0.43–2.30) 2 Thermal Bands (10.6-12.51)	Cirrus, Coastal/Aerosol, Blue, Green, Red, NIR, SWIR1-2: 30 m Thermal Infrared 1–2: 100 m Panchromatic: 15m	16

In this research, all 606 accessible imageries between 2000 and 2021 were used to examine the variations in lake extents at finer temporal scales (Figure 2), except for 2012 due to no images being available in that year. Previous studies [30,33,35,41] have used short- and long-term data in monitoring the changes that occurred in the water surface area of the lakes. Due to describing long-term patterns of surface water fluctuation, images of different years must be captured at the same time zone (the same period of a year), which might lead to inconsistency in surface water extraction and uncertainty in the results. In addition, identifying the best time of year for a single image and finding the images that are acquired in the same period of different years are also problematic. Thus, using all accessible images within a year helps avoid these mapping difficulties, discrepancies, and ambiguity. Consequently, comprehensive research is required to ascertain the ongoing change of surface water areas utilizing adequate images acquired over several seasons and years.

The LSWT refers to the temperature of the lake’s surface water, which is a critical indicator for the lake’s hydrology and biogeochemistry. Surface temperature patterns throughout time could provide insight into how climate change is affecting the lake area. A relationship existed between the extent of lake surface area lost and the rise in LSWT measurements, according to Aslan and Koc-San [34]. LST studies commonly use Landsat images, which are among the most widely utilized types of RS data. In this study, LSWT data for the three lakes were calculated using Landsat 5 TM and Landsat 8 OLI/TIRS thermal bands and inversion of Planck function through the GEE. These bands have native spatial resolution of 120 m and 100 m for Landsat 5 and Landsat 8 satellites, respectively, but they are provided after resampling to 30 m by the United States Geological Survey (USGS).

The average surface temperature of the world has been rising over the previous century, and it is projected that future climate change will have a detrimental effect on aquatic systems [42]. Climate change has a negative impact on water bodies in Türkiye and the rest of the world [43]. According to Abujayyab et al. [33]’s study, there is a substantial positive or negative relationship between the reduction of surface water area and meteorological factors. Therefore, to indicate whether there are correlations between water surface area and LSWT with meteorological parameters (temperature, evaporation, and precipitation), ERA5-Land and TerraClimate monthly satellite dataset were used in this study. ERA5-Land dataset is produced from ERA5 land variables with enhanced resolution. The processing steps consist of modeling data using globally consistent observations considering the laws of the physics, and they have been providing data (50 variables) in a uniform grid since 1981 with approximately 11 km spatial resolution. Monthly data have also been produced by taking an average of hourly collected data [44]. The TerraClimate monthly dataset is produced by combining interpolated time-varying data CRU Ts4.0 and the Japanese 55-year Reanalysis (JRA55) with a high spatial resolution WorldClim dataset. Unlike ERA5-Land, TerraClimate grid has nearly 5 km spatial resolution and 14 variables from the years 1958 to 2021 [45]. Both datasets exist in the database of the GEE cloud platform, thus the corresponding data were retrieved from the GEE database for this research.

3.2. Methods

Various techniques have been developed to extract surface water area information since it is required to convert remotely sensed images into information for practical applications, such as land cover and surface water mapping. Depending on the number of the utilized bands, thematic water surface area extraction algorithms can be split into two main categories: single-band and multi-band methods. The multi-band methods, which are based on the spectral water index, are better than single-band water indices in water extraction because of the advantage of using various reflectivity variations in the spectrum that enables separating water and other land cover classes. Meanwhile, water indices have been successfully utilized to extract surface water resources using remotely sensed data, and they have the benefits of accurate, simple, quick, and repeatable extraction of surface water information to capture the significant intra-annual and inter-annual water fluctuation [4,46]. The methodology’s process flow of this study is portrayed in Figure 3.

3.2.1. Water Surface Area Extraction

The water index method is used in conjunction with other RS indices to determine the substantial difference between a water body and a non-water body where water extraction is quite precise. When extracting large surface water extent, this approach has the benefits of speed, precision, simplicity, resilience, and reproducibility [23]. The NDWI was introduced by McFeeters [47] as in Equation (1). It is a spectral index that maximizes water reflection in the green band while minimizing water reflectance in the near-infrared (NIR) band [48]. The NDWI scale ranges from −1 to +1, where NDWI > 0 denotes water, while NDWI ≤ 0 represents non-water.

This study used NDWI to extract water surface areas from Landsat satellite images, since many studies [46,49,50,51] indicated that NDWI provided satisfactory results. The threshold value of zero was suggested for extracting the surface water pixels with all positive NDWI values classified as water and all negative and zero NDWI values classified as non-water.

NDWI = (Green − NIR)/(Green + NIR)

(1)

where Green is the reflected green visible light and NIR is the reflected near-infrared energy.

3.2.2. Accuracy Assessment of Extracted Water Surface Area Map

The purpose of this section is to demonstrate the efficiency of the proposed surface water body extraction technique and to verify the quality of the lake boundary extraction data. To determine the accuracy assessment of the NDWI zero threshold in GEE for each lake, the NDWI map was clipped according to the boundaries. Then the images, which have the smallest and largest lake area, were processed. Random points were selected for each NDWI-based lake map by taking considering that the dots covered all lake surface area and its boundaries (Figure 4). Thus, a total of 318 points were chosen for Lake Burdur, 403 points for Lake Egirdir, and 495 points for Lake Beysehir from high-resolution Google Earth images. The accuracy points were then compared with the lake over multi-spectral Landsat image to verify whether they were in the water or not.

A confusion matrix was generated as portrayed in Table 3, to compute the producer accuracy, user accuracy, and overall accuracy given by Equation (2) [52]. Additionally, the Kappa coefficient was computed to evaluate the degree of agreement between the reference data and the random classifier. It could vary between −1 and +1, where 1 represents complete agreement [53]. In order to perform accuracy assessment using confusion matrix, four variables below should be identified:

TP: True positive value: the number of correct water pixels.

FP: False positive value: the incorrect numbers of the water pixels.

FN: False negative value: undetected number of the water pixels.

TN: True negative value: the number of non-water pixels that were correctly rejected.

Overall Accuracy (%) = (TP + TN) × 100/Total Number of Pixels (T)

(2)

Table 3. The confusion matrix.

	Category/Class	Reference Data
		Water	Non-Water
Classified Data	Water	TP	FP
	Non-water	FN	TN

Table 3. The confusion matrix.

	Category/Class	Reference Data
		Water	Non-Water
Classified Data	Water	TP	FP
	Non-water	FN	TN

Producer accuracy, user accuracy, and overall accuracy reflect accurate predictions and have a range of 0 to 100%, with the value close to 100% indicating perfect accuracy.

3.2.3. The Determination of LSWT Variable

In the study, the inversion of the Planck function, which is one of the most accurate LST retrieval methods with an average RMSE equal to 2.29 K [54], was used to extract surface temperature from the Landsat imageries [55]. For the satellite imageries, a cloud filter was applied on the GEE, and atmospheric correction (which corrects for brightness/contrast, the effects of solar elevation angle, and Earth-sun distance effect) was implemented via radiance and reflectance conversions. Furthermore, the utilized images were checked via quality assessment (QA) bands in the GEE. The initial step in this approach was to compute the brightness temperature values after calculating spectral radiance using formulas provided by a previous study [56] based on the thermal bands of each sensor. After that, the normalized vegetation index threshold method and the emissivity adjusted LST values were used to calculate the LSWT values [57]. Land surface emissivity (ε) is the effectiveness of thermal energy that is transferred from the surface to the surrounding atmosphere [58,59]. All applications were applied using the GEE code for each Landsat sensor based on its thermal band. Landsat 5 has one thermal band (band 6) whereas Landsat 8 has two thermal bands (band 10 and 11). Thus, band 6 of Landsat 5 and band 10 (TIR1) of Landsat 8 were used in the application. The following fundamental equations are utilized in LSWT determination [Equations (3)–(7)]:

L_(TOA) = (M_L × Qcal) + A_L

(3)

where L_(TOA) is Top of Atmospheric spectral radiance value in watts/m² × sr × µm; M_L is band specific multiplication rescaling factor equal to 0.055375 for Landsat 5, and 0.0003342 for Landsat 8; Qcal is quantized and calibrated standard product pixel values (DN); A_L is resizing factor equal to 1.18243 for Landsat 5, and 0.1 for Landsat 8.

T_B = K₂/ln(K₁/(L_(TOA) + 1))

(4)

where T_B is the Brightness Temperature; K₁ (watts/m² × sr × µm) and K₂ (K) are band-specific thermal conversion constant obtained from the metadata file: K₁ = 607.76 and K₂ = 1260.56 for band 6 of Landsat 5 TM, while K₁ = 774.89 and K₂ = 1321.08 for band 10 of Landsat 8 OLI/TIRS.

LSWT = T_B/(1+ ((λ × T_B/α) × ln(ε)))

(5)

where LSWT is lake surface water temperature; ε is surface emissivity; λ is average wavelength of band; α = 1.438 × 10⁻² mK. Surface emissivity ε can be estimated using NDVI-based threshold method, as portrayed in

Table 4. Surface emissivity values assigned based on NDVI thresholds.

NDVI Threshold	Land Cover Type	Surface Emissivity
NDVI < 0	Water	0.985
0 ≤ NDVI ≤ 0.1	Bare Soil	${\rho }_R$ (red reflectance band)
0.1 ≤ NDVI ≤ 07	Vegetation Mixed with Soil	0.990 × Pv + 0.984 × (1 − Pv) + 0.04 × Pv × (1 − Pv)
NDVI > 0.7	Vegetation	0.990

[57,60,61].

${\rho }_R$ is the reflectance band in the Red region of the electromagnetic spectrum. NDVI is obtained using the Red and NIR reflectance bands as below:

NDVI = (NIR − Red)/(NIR + Red)

(6)

Pv [62] is vegetation proportion, and calculated by:

Pv = ((NDVI − NDVImin)/(NDVImax − NDVImin))²

(7)

In order to evaluate the LSWT results, MODIS LST data were utilized for cross-validation since there is no ground station in the test sites. For the evaluation, the average LSWT time series were compared with MODIS LST and the correlation results were presented for each lake.

4. Results and Discussion

4.1. Accuracy Assessment of Extracted Water Surface Area Map

The extracted results from the Landsat images, clipped based on each lake border, were compared point-by-point with the NDWI map for the minimum and maximum water surface area maps to demonstrate the accuracies of our results.

Table 5. The accuracy assessment results of extracted water pixels with NDWI.

	Burdur Lake 10 October 2021 Landsat 8	Burdur Lake 20 May 2004 Landsat 5	Egirdir Lake 10 October 2021 Landsat 8	Egirdir Lake 17 March 2004 Landsat 5	Beysehir Lake 17 September 2001 Landsat 5	Beysehir Lake 27 April 2013 Landsat 8
	Smallest Area	Largest Area	Smallest Area	Largest Area	Smallest Area	Largest Area
Area (km2)	120.53	159.29	439.85	467.00	605.89	685.01
User Acc. (Water)	99%	99.1%	98.3%	98%	100%	97.7%
Producer Acc. (Water)	100%	100%	100%	100%	99%	100%
Overall Acc.	99.3%	99.45%	98.7%	98.5%	99.3%	98.7%
Kappa	0.983	0.988	0.959	0.963	0.984	0.973

displays the smallest and the largest lake area values in addition to their accuracy for users, producers, and overall.

After generating the confusion matrix for maximum and minimum water surface area in each lake during a 21-year period, the overall accuracy of the results was more than 98%, as the producer accuracy and user accuracy for water areas was more than 97%. Moreover, Kappa coefficient values were close to 1, representing perfect agreement. Considering the suggested extraction technique outcomes in this study, the evaluation indicates a high degree of accuracy in water extraction, which is indicative of its effectiveness and simplicity. It also demonstrates that this approach can be used in other relevant geographic regions with identical backgrounds for future research.

4.2. Spatial and Temporal Changes of Lake Water Extent from 2000 to 2021

4.2.1. Lake Burdur

Lake Burdur has a brackish character and is high in arsenic, making it unsuitable for residential and agricultural irrigation [63]. According to the examination of long-term remotely sensed images, water surface area has decreased considerably between 2000 and 2021. Temporal changes in lake water surface area and its trend between 2000 and 2021 are portrayed in Figure 5. There is no doubt that throughout the course of the last two decades, the area of the lake has steadily dropped, except for the period from 2003 to 2004. After 2004, the water surface area continued to decrease at a steady rate, and the annual area trend has been decreasing at the same rate until 2021. From the first to the end period of the study, the water surface area decreased by more than 22% (35 km²) from average area of 156.1 km² in 2000 to 121.4 km² in 2021. To achieve the obtained results, time series trend analysis with ANOVA test was applied to indicate if this trendline result was statistically significant, and the test results proved this value at a significance level of p < 0.01.

Figure 6a illustrates the seasonal water surface area changes between 2000 and 2021. The lake reached its maximum water surface area value in the spring of 2004 by 159.3 km² and it has its lowest value in the summer of 2021 by 120.5 km². From Landsat images, it can be observed that the northeastern portion of the lake experiences the greatest drought (Figure 6a). The missing bars in Figure 6a indicate that there is no image or cloud-free image during the given season in that year.

Generally, except for the year 2004 when the lake’s surface area achieved its highest values and recovered by 0.68%, the lake surface area decreased most of the time as indicated in Figure 6b where the variations in annual water surface areas were compared with respect to the water surface area in 2000. The highest reduction occurred in 2021 when the lake’s surface area decreased by 22.21%.

The primary causes for such a massive change would be the impacts of climate change on the basin, and the exaggerated use of water sources that were meant to feed the lake as the lake is close to the city center. Lake Burdur, located in a closed basin, is fed by precipitation, water flowing into the lake from rivers/streams and groundwater. Over the years, dams and ponds were erected on streams/rivers that deliver water to Lake Burdur, thus the only water source of the lake is precipitation. As a result, after the construction of the dams and ponds over the rivers/streams that are the main sources feeding the lake, the water level of the lake decreased fast due to low precipitation and high evaporation rates.

4.2.2. Lake Egirdir

Lake Egirdir is a freshwater lake in the Isparta region that was created by tectonic and karstic forces. Along with serving as a natural drinking water basin, it is a significant wetland considering its biological benefits and irrigation uses. In the study period for Egirdir Lake, it was discovered that the surface area decreased slightly. During the first few months of 2004, the lake’s water surface area reached its peak level of 464 km², and it reached its lowest surface area value of 440 km² in the autumn of 2021 as observed in Figure 7, which presents a time series for the water surface area changes during the 21-year period. The water surface area value of the lake, which was at its maximum in 2004, decreased by more than 3% (15 km²), from average area of 458 km² in 2000 to 443 km² in 2021. For this result, time series trend analysis with an ANOVA test was applied to demonstrate if this trendline result was statistically significant, and the test results proved this value at a significance level of p < 0.01.

The seasonal changes in water surface area are portrayed in Figure 8a. The lake’s surface area increased in the winter and spring seasons and decreased in the summer and autumn, indicating a seasonal fluctuation in the lake’s extent. This was attributed to the fact that during the summer months, irrigation purposes led to the reduction in surface water area. Figure 8a also portrays the spatial changes in water surface area for the first and last assessment intervals. The missing bars in Figure 8a indicate that there is no data or cloud-free image during the given season in that year.

The variations in Egirdir Lake extent are just 5–15 km² during the study period, which only indicates a 3–4% deviation from the extent in 2000. As observed in Figure 8b, the annual variations in lake water surface area for each year were compared with respect to the water surface area in 2000. The highest reduction within the study period was observed in 2021 when the lake’s surface area decreased by 3.39%. When comparing Burdur and Egirdir Lakes, Lake Egirdir was almost three times larger than Lake Burdur in 2000, and it is not as close to the city center as Lake Burdur. Additionally, the number of constructed ponds and dams over the rivers/streams feeding the lake are not as many as Lake Burdur. These circumstances may have contributed to the surface water area of Lake Egirdir by being less impacted by human activities than Lake Burdur.

4.2.3. Lake Beysehir

The third-biggest lake in Türkiye, Lake Beysehir, is a freshwater lake located in the province of Konya. In 1992, the lake was designated as a natural protected area, and it is used for irrigation and domestic purposes. In Lake Beysehir, after analyzing the RS images, almost no change was determined in the water surface area. During the study period between the years 2000 and 2021, the average area trend increased slightly by less than 1% at a statistically significant level (p < 0.01) as portrayed in Figure 9. An area of 685 km² has been determined as the maximum value occurring in the spring of 2013, which is the year that experienced a significant increase throughout study period. This increase was due to the rise in the amount of precipitation in that year. Despite the slight increase in the direction of water surface extent during the 21-year period, the water surface area, which was 627 km² in 2000, decreased to 623 km² in 2021.

In contrast to the summer and autumn seasons, the lake surface area values increased to its highest levels ever recorded in the winter and spring seasons, as demonstrated in Figure 10a, where between spring and autumn, there have been a few modest changes. Spring results have generally been somewhat greater than autumn ones; this was because water was used for irrigation purposes. The missing bars in Figure 10a indicate that there is no data or cloud-free image during the given season in that year.

Generally, reductions in lake surface area occurred most frequently before 2013 as observed in Figure 10b, which portrays the percentage of annual variation for water surface area in comparison to lake surface area in 2000. The lake’s surface area decreased by 2.25% with water surface area value of 606 km² in 2021, which was recorded in the autumn season. While the highest recovery rate that occurred in the surface water area of the lake was in 2013, it increased by 3.18% annually.

4.3. Relationship between LSWT and Lake Water Extent Changes

The LSWT values for the three lakes were examined, and an upward trend was observed for Lake Burdur and Lake Egirdir, whereas a downward trend was observed for Lake Beysehir. The annual average value of Lake Burdur was 19.00 °C in 2000, whereas it increased to 21.10 °C in 2021 while increasing by 2.13 °C. While it was 15.8 °C in 2000, it reached 16.1 °C in 2021 for Lake Egirdir. For Lake Beysehir, the value was 17.00 °C in 2000, while it decreased to 15.50 °C in 2021 as observed in Figure 11, which portrays the annual mean of LSWT during the study period. No images were available for the year 2012. As a result, the year’s water temperature value could not be computed. Some images were identified as outliers due to their high deviation from the mean values, which were likely caused by cloud presence or atmospheric conditions.

In 1998, based on Aslan and Koc-San [34], Burdur Lake’s area was roughly 161 km², but by 2018, it had shrunk by 127 km². Therefore, it has lost more than 20% of its surface area in the last two decades, while the LSWT value increased by 2.2 °C throughout the same period. According to our study after evaluating the data, the surface water temperature in Burdur Lake has increased and the trend tends to upward. Comparison of the 2000 and 2021 datasets indicates that the lake surface area decreased by 35 km², with an increasing LSWT of 2.13 °C as portrayed in Figure 12a, which indicates the changes in the lake surface extent with LSWT.

The average values of LSWT for the summer and autumn seasons in the years 2000 and 2021 have been evaluated and extracted based on the NDWI maps of the corresponding year to indicate the spatial changes in LSWT. The LSWT values varied between 20.4 °C and 35.7 °C in the summer of 2000, while they increased between 22.4 °C and 35.9 °C in the summer of 2021 as portrayed in Figure 12b,d, which indicate the LSWT maps for summer season of years 2000 and 2021, respectively. In the autumn season, the LSWT values varied between 11.3 °C and 20.5 °C in 2000, while they increased between 15.5 °C and 24.1 °C in the same season of 2021, as portrayed in Figure 12c,e. The lake lost more than 22% of its surface area during the 21-year period, indicating that the LSWT and the surface extent of the lake were found to be in a directly proportional relationship.

Aslan and Koc-San [34] analyzed the relationship between the lake surface area and the change in LSWT, stating that the LSWT values increased as the lake extent decreased between 1998 and 2018. The LSWT values of Egirdir Lake rose by more than 1.3 °C and the lake’s surface area reduced by less than 4 km². Based on the comparisons of the 2000 and 2021 datasets in our study, the LSWT of Egirdir Lake has increased by more than 0.32 °C, while the water surface area decreased by 15 km² throughout the study period from 2000 to 2021, as portrayed in Figure 13a, which demonstrates the changes in lakes surface area with LSWT.

To portray the spatial changes that occurred in water surface area due to a decrease or increase in the water temperature, the average values of LSWT for the summer and autumn seasons in the years 2000 and 2021 have been calculated. The LSWT values varied between 17.4 °C and 34.8 °C in the summer of 2000, while they increased between 20.9 °C and 39.00 °C in the summer of 2021 as portrayed in Figure 13b,d. In the autumn season, the LSWT values varied between 7.9 °C and 24.1 °C in 2000, whereas they increased between 9.8 °C and 22.7 °C in the same season of 2021, as can be observed in Figure 13c,e. The lake has lost less than 3% of its surface area during the 21-year period, indicating that LSWT could have a slight impact on the lake’s extent.

For Lake Beysehir, according to the 2000 and 2021 datasets, the LSWT has a downward tendency. Figure 14a depicts a decreasing trend in the LSWT by about 1.5 °C at a statistical significance level (p < 0.01), and the water surface area also decreased by less than 3 km² throughout the study period.

The spatial changes in water surface area with the average LSWT for the summer and autumn seasons of 2000 and 2021 have been extracted. The LSWT values varied between 20.6 °C and 39.9 °C in the summer of 2000, while they ranged between 21.3 °C and 31.9 °C in the summer of 2021, as portrayed in Figure 14b,d. In the autumn season, the LSWT values varied between 11.9 °C and 26.7 °C in 2000, whereas they increased between 15.4°C and 35.3 °C in the same season of 2021 as portrayed in Figure 14c,e. The surface area variation was around 1% during the 21-year period, indicating that LSWT could not have more impact on the lake’s extent.

Since there was no ground-based LST station in the region, we compared the average LSWT time series with MODIS LST for each lake. Figure 15 represents the scatter plot and correlation values between LSWT and MODIS LST for all three lakes. The correlation results were above 0.93, revealing that the obtained LSWT values are highly compatible with MODIS LST data.

4.4. Meteorological Parameters and Correlation Analysis with LSWT and Lake Water Extent

To present meteorological effects on the obtained time series results, correlation analyses of LSWT and surface water area with meteorological parameters (temperature, precipitation, and evaporation) were performed. For this purpose, the correlation coefficient (R) was determined between LSWT and water surface area extent obtained from Landsat satellite images and meteorological parameters gathered from ERA5-Land and TerraClimate datasets for each lake. ERA5-Land has 50 available variables, but only 4 parameters (temperature above 2 m, temperature of water surfaces (lake mix layer temperature), total evaporation, and precipitation) were used in the study. The TerraClimate dataset has 14 bands, whereas four bands, namely, total precipitation, evapotranspiration, and temperature bands (minimum and maximum temperatures) were utilized in the study. The average temperature was also calculated from minimum and maximum temperature bands.

With the aim of evaluating time series statistics of temperature datasets between the years 2000 and 2020, boxplots were created for each lake, which are portrayed in Figure 16. The boxplots are a summary of data distribution including maximum value, first quartile of data, median and mean values, third quartile of data, and maximum value. Mean value is indicated with ‘x’ in the boxplots, and the points above the minimum/maximum values (Figure 17a) are the outliers in the dataset. In the boxplots, LSWT and lake mix layer temperature are the surface temperature of the lake, whereas the temperature (2 m) of ERA5-Land and the average temperature of TerraClimate represent the air temperature.

While maximum LSWT values of Lake Burdur demonstrate nearly close values to the lake mix layer temperature, the maximum values of LSWT in Egirdir and Beysehir Lakes are lower than lake mix layer temperature. Additionally, maximum LSWT from largest to smallest can be sorted as Lake Burdur, Lake Beysehir, and Lake Egirdir, respectively. The primary reason could be the geographic locations (elevation) of the lakes, provided in Table 1. Concerning the elevations, Beysehir Lake is in the highest location among the lakes with 1124 m. Egirdir and Burdur Lakes are situated at 917 m and 857 m, respectively. Conversely, the secondary reasons may be related to the shallowness of the water and mixed pixels on the edge of the shoreline.

The other important parameters are precipitation and evaporation. According to statistics of these parameters for each lake (Figure 17), Burdur Lake has lower precipitation and higher evaporation than Egirdir Lake. Likewise, Egirdir Lake has lower precipitation and slightly higher evaporation than Beysehir Lake. However, the maximum precipitations of the lakes are quite different and Beysehir Lake has the highest precipitation value. These findings also reveal that precipitation and evapotranspiration statistics are meaningful when considering the water surface area results.

The monthly correlation analysis was performed with produced data (water surface area and LSWT) and climate parameters (temperature, precipitation, and evaporation/evapotranspiration). Obtained R values from correlation analysis are provided in

Table 6. The monthly correlation analysis of meteorological parameters with LSWT and water surface area.

	Landsat		ERA5-Land				Terraclimate
Lake	Monthly Correlation	LSWT	Temperature (2 m)	Lake mix layer temperature (C°)	Total evaporation (m)	Precipitation (mm)	Total Precipitation (mm)	Evapotranspiration (mm)	Average Temp
Burdur	Water surface area	−0.207 **	0.018	0.012	−0.086	−0.107	−0.061	0.095	0.028
	LSWT		0.507 **	0.478 **	−0.567 **	−0.224 **	−0.194 **	0.576 **	0.524 **
Beysehir	Water surface area	−0.418 **	−0.602 **	−0.584 **	0.529 **	0.487 **	0.430 **	−0.332 **	−0.480 **
	LSWT		0.806 **	0.858 **	−0.929 **	−0.349 **	−0.651 **	0.857 **	0.944 **
Egirdir	Water surface area	−0.768 **	−0.729 **	−0.734 **	0.554 **	0.386 **	0.346 **	−0.564 **	−0.716 **
	LSWT		0.964 **	0.972 **	−0.872 **	−0.599 **	−0.530 **	0.872 **	0.960 **
			very high	high	moderate	low	very low
			** Correlation is significant at the 0.01 level (2-tailed).

. R values were classified as: very high (0.8 ≤ R ≤ 1.0), high (0.6 ≤ R < 0.8), moderate (0.4 ≤ R < 0.6), low (0.2 ≤ R < 0.4), and very low (0 < R < 0.2). Additionally, a 2-tailed statistical significance test was applied to reveal the degree of the relationship among the parameters. Statistical significance is often referred to as the p-value (probability value). The p value allows us to determine if there is a statistically significant correlation between two variables. As a general assumption, a p < 0.05 indicates statistical significance while p < 0.01 is explained as high statistical significance. It was discovered that there is negative correlation between the LSWT and water surface area. The monthly correlation between water surface area and LSWT was discovered highly statistically significant [(p < 0.01) for all three lakes]. However, the annual correlation was discovered to be statistically significant (p < 0.05) for Egirdir Lake, and the annual correlation for Burdur and Beysehir Lakes was statistically insignificant (p > 0.05). Because increasing air temperature and evaporation have affected Burdur Lake’s water surface area [64], and shrunk lake area has slightly increased the LSWT in time. One should also remember the inertia of hydrological systems that have an impact on the obtained correlation results.

Water surface area has also a statistically highly significant (p < 0.01) monthly correlation with temperature and evaporation in Beysehir and Egirdir Lake, unlike Burdur Lake. Although it has the lowest precipitation and the highest temperature, there is no correlation between water surface area and all meteorological parameters of both datasets in Lake Burdur, although there was a 22% loss in water surface area. However, low and moderate correlations were determined on the annual scale.

TerraClimate average temperature indicates a highly statistically significant (p < 0.01) correlation with LSWT for all lakes. The precipitation and temperature dataset in TerraClimate indicates higher correlation than ERA-5 Land in Beysehir Lake. LSWT indicated a highly statistically significant (p < 0.01) correlation with evaporation and evapotranspiration values, except for Beysehir lake, and the same pattern was observed for its correlation with precipitation in both datasets for all lakes. According to the results, it can be concluded that the sensitivity degree of Egirdir Lake to temperature and evaporation is identical with Beysehir Lake. It could be related to the similar loss amount of water in these lakes. For all lakes, the relationships between water surface area and other meteorological parameters are high for Beysehir Lake and Egirdir Lake, unlike Burdur Lake.

5. Conclusions

The primary purpose of this research is to use RS images to examine the dynamic changes in surface water areas of natural lakes. Additionally, the relationship between changes in water extent and the LSWT was explored to facilitate a better understanding of the causes of surface water extent changes. RS analysis was performed to facilitate the rapid and effective processing of many satellite images within the GEE cloud environment that enables creating a new database for the spatio-temporal changes of natural lakes. The advantages of generating such a database via the use of new technological tools such as RS include the long-term changes in water surface areas for the attention of decision-makers and other authorities responsible for the sustainable management of water resources. Analyses were conducted in this manner on the GEE platform without the need to download the satellite imageries. RS images from the Landsat 5 TM and Landsat 8 OLI/TIRS collections were used to illustrate the research region and 10% cloud threshold was applied to these images based on each lake’s borders. Additionally, the quality of all images was checked using the QA bands.

This study includes a total of 606 Landsat images. Landsat 5 TM imagery was accessible from 2000 to 2011, whereas Landsat 8 OLI imagery was available from 2013 until 2021, thus, there were no available images in 2012. To identify which pixels should be assessed as water, the NDWI was used, where NDWI greater than zero was considered to be water. In the GEE code editor, NDWI values for all images in the collection were calculated and utilized to extract the water surface area. The results indicate that the water pixels can be extracted quite well using NDWI with overall accuracy above 98% for all the lakes’ areas. It is also indicated that the decrease in water surface area has a direct correlation with the rise in LSWT values, especially for Lake Burdur and Lake Egirdir. Based on the comparison between the 2000 and 2021 datasets, Lake Burdur has lost almost 22% of its water surface area while the LSWT has increased by more than 2.13 °C. In Lake Egirdir, the water extent has decreased slightly by less than 3%, and the LSWT increased by more than 0.32 °C from 2000 to 2021. The water surface extent did not change significantly for Lake Beysehir at the end of the study period, where the annual variations took a fluctuating character; however, the LSWT decreased by about 1.5 °C.

Correlation analyses were conducted for the water extent of each lake and LSWT with meteorological parameters of two different datasets, namely, ERA-5 Land and TerraClimate. Correlation results indicate that the LSWT has a very high correlation with air and surface temperature for Lake Beysehir and Lake Egirdir, but a moderate correlation for Lake Burdur. Considering the water surface area, Lake Beysehir and Lake Egirdir have a high correlation with ERA5-Land temperature datasets, unlike Lake Burdur. Burdur Lake’s water surface area has no correlation with meteorological parameters, meaning the drought in Lake Burdur is not directly related to the meteorological parameters. The common issues for water surface area changes were excessive water withdrawals from the lakes themselves, primarily those of freshwater nature, lake usage for industrial reasons, and climate change effects. Overall, the variations were more likely due to human-induced activities, especially for Lake Burdur. As a future perspective, all lake regions in Türkiye and in the world can be analyzed with the same time series analysis to reveal their spatio-temporal changes. New datasets and methodologies may be compared to discover the best solutions for this research design.