Multi-source remote sensing approach for monitoring the Caspian Sea shoreline dynamics

This study presents a scalable methodology for the semi-automated extraction of the Caspian Sea shoreline and water surface based on medium-resolution satellite imagery and cloud computing within the Google Earth Engine environment. The input datasets included very high-resolution PlanetScope imagery, Harmonized Landsat–Sentinel products, and DAHITI satellite altimetry data. The proposed approach consists of three stages: generation of reference and working shoreline datasets, validation-based selection of the optimal extraction algorithm, and its subsequent scaling to the entire Caspian Sea. To assess accuracy, twelve automated water–land boundary extraction scenarios based on different spectral indices and temporal compositing methods were evaluated. Validation on test sites in the northeastern part of the Caspian Sea was performed using transect analysis in Digital Shoreline Analysis System and reference shorelines derived from PlanetScope imagery. The results demonstrated that the highest positional accuracy was achieved by combining the MNDWI with the 25th-percentile temporal compositing approach. Meanwhile, the AWEI2 method with median aggregation produced a more spatially coherent water mask and was used as an additional quality-control layer to identify local artifacts and gaps in the MNDWI_p25 mask during subsequent manual correction. Based on the selected algorithm, a multi-year geospatial dataset of shorelines and water surface polygons of the Caspian Sea was generated for the period 2016–2025. The results showed that the sea level declined by 1.33 m, the water surface area decreased by 11,216.31 km 2, and the shoreline length was reduced from 19,907.94 km to 15,551.79 km. The most pronounced changes, reflecting the high sensitivity of shallow-water regions to fluctuations in water level, were observed in the northeastern Caspian Sea, where shoreline retreat exceeded 20 km in some areas. The obtained results confirm the effectiveness of the proposed approach as a basis for scalable monitoring of coastal-zone dynamics and demonstrate its potential for geoecological assessments, forecasting morphodynamic processes, and developing adaptation strategies under conditions of continued Caspian Sea level decline.

Coastal zones represent dynamic natural systems that are continuously transformed under the influence of hydrological, morphodynamic, and climatic factors and are subject to erosion and sediment accumulation processes, as well as the inundation and desiccation of shallow-water areas (Timmerman et al., 2021). In this context, the monitoring of water bodies cannot be limited solely to the analysis of surface area; it also requires detailed tracking of shoreline position as a fundamental indicator of ongoing changes (Saddam et al., 2025). The importance of such studies becomes particularly critical in shallow-water regions, where even minor sea-level fluctuations lead to substantial horizontal displacements of the water–land boundary and profound transformations of coastal ecosystems (Mentaschi et al., 2018; Thiéblemont et al., 2024). Under these geographical conditions, the shoreline serves as one of the most sensitive indicators of both short-term hydrometeorological processes and long-term climatic trends. Regular monitoring of shoreline position enables the quantitative assessment of coastal degradation or recovery rates and facilitates the timely identification of areas exposed to elevated environmental risks (Vitousek et al., 2023).

Modern remote sensing techniques provide unique opportunities for systematic and objective shoreline monitoring over large spatial and temporal scales (Chen et al., 2019; Vos et al., 2019). The use of optical satellite data with different spatial and temporal resolutions, together with radar and LiDAR observations, enables not only high-precision mapping of the waterline but also comprehensive analyses of its spatiotemporal dynamics and the identification of patterns of coastal transformation (Toure et al., 2019; Zhou et al., 2023). Consequently, the development and improvement of automated shoreline extraction methods have become one of the key tasks in geospatial analysis of water bodies and coastal zones (Xu, 2018), particularly in enclosed and semi-enclosed marine basins (Duan et al., 2025; Rahbani et al., 2026).

The primary source of data for long-term shoreline monitoring has traditionally been the Landsat program (missions 5, 7, 8, and 9), which provides the longest continuous observational record since the 1980s (Wulder et al., 2019). The satellites of this series offer a spatial resolution of 30 m in the multispectral bands, allowing the reconstruction of historical coastal dynamics over more than four decades. During the last decade, shoreline studies have been significantly advanced by the availability of Sentinel-2 (A/B) data within the Copernicus programme (Drusch et al., 2012). Its higher spatial resolution (up to 10 m) in the visible and near-infrared bands, combined with a high revisit frequency (up to 5 days), enables the detection of shoreline dynamics with improved spatial accuracy (Vos et al., 2023). In addition to optical data, Sentinel-1 synthetic aperture radar (SAR) imagery has been increasingly used. Owing to its capability to acquire observations independently of cloud cover and illumination conditions, these data are employed to fill gaps in observation series and to refine shoreline position under adverse weather conditions (Wei et al., 2021). To validate shoreline dynamics, satellite altimetry data derived from the successive TOPEX/Poseidon, Jason-1/2/3, and Sentinel-6 missions are widely used. These observations provide highly accurate estimates of absolute sea-level changes and are essential for establishing relationships between vertical water-level fluctuations and horizontal shoreline displacements (Donlon et al., 2021). Furthermore, shoreline extraction from very high and high spatial resolution imagery provides a more detailed and geometrically accurate representation of the water–land boundary, particularly in areas characterized by complex coastal morphology and shallow water environments. For this purpose, imagery from WorldView (Maxar) (Maglione et al., 2014), Pléiades (Airbus) (Diatta et al., 2025), and PlanetScope (Planet) (Tan et al., 2026) satellites are commonly employed.

The interpretation of optical imagery requires the application of specialized algorithms. Spectral-index-based methods remain the most widely used approach and, due to their computational simplicity and robustness under varying imaging conditions, constitute the standard framework for automated shoreline extraction (Zhou et al., 2023). The normalized difference water index (NDWI) and its derivatives are traditionally used and have demonstrated high effectiveness (Xu, 2006). In particular, the combined use of the water index NDWI and the vegetation index NDVI helps minimize classification errors in regions with dense vegetation cover and wetlands (Chowdhury and Yang, 2025). The reliability of water indices for shoreline extraction has been demonstrated for a wide range of aquatic environments, from river deltas to marine coastlines (Vos et al., 2023; Christofi et al., 2025). The automation of shoreline extraction is commonly achieved through the application of Otsu’s thresholding method, which determines the optimal threshold for class separation (Duan et al., 2025). To improve robustness against shadows and dark surfaces, the Automated Water Extraction Index (AWEI) is also widely applied (Feyisa et al., 2014).

Despite the effectiveness of index-based approaches, their accuracy is inherently limited by the physical pixel size, resulting in discretized shoreline contours. Overcoming this limitation has led to the development of sub-pixel analysis techniques (Palomar-Vázquez et al., 2023), aimed at reducing the constraints imposed by the spatial resolution of satellite imagery. Within this framework, the shoreline is treated as an isoline of a continuous spectral signal, allowing its position to be determined within an individual pixel through interpolation procedures. The effectiveness of this approach for monitoring shoreline dynamics has been confirmed by numerous studies demonstrating improved stability in water surface extraction (Pardo-Pascual et al., 2012).

To minimize the effects of cloud cover and ensure a spatially continuous representation of the shoreline, multi-temporal pixel-based aggregation methods based on the generation of composite images are increasingly employed. This approach is commonly implemented using the computational capabilities of the Google Earth Engine (GEE) platform (Gorelick et al., 2017) and provides datasets that are more robust to noise and missing observations. Its effectiveness for shoreline dynamics analysis has been demonstrated in a number of studies (Vos et al., 2019; Ding et al., 2021; Hamzaoglu and Dihkan, 2023; Karimzadeh, 2023).

For areas characterized by complex spectral properties, where index-based approaches may be insufficiently accurate, machine-learning algorithms are increasingly applied to achieve more flexible class separation. Among the most widely used methods are Random Forest (Liu et al., 2022) and Support Vector Machine (SVM) (Kalkan et al., 2013).

The current stage in the development of water–land boundary detection methods is associated with the active implementation of deep learning technologies based on convolutional neural networks. The use of deep segmentation architectures (such as U-Net and DeepLabV3+) enables high-precision delineation of water surfaces even under complex conditions and when using data from different sensors, as demonstrated in studies (Li et al., 2021; Blais and Akhloufi, 2025). Despite their high accuracy, these methods require substantial volumes of training data and significant computational resources.

An important methodological direction in recent research is the intercomparison of satellite data with different spatial resolutions, which allows for an objective assessment of shoreline extraction accuracy and improves the reliability of the final results. In this context, very high-resolution imagery is used as a reference basis (“ground truth”) for calibration and validation of results obtained from medium-resolution images (Landsat, Sentinel-2). The effectiveness of this approach is confirmed by studies (Liu et al., 2017; Dike et al., 2025). The application of high-detail imagery is particularly important for areas with complex coastal morphology and extensive shallow-water zones, where standard water–land boundary extraction methods are limited by pixel spectral mixing effects.

It is precisely under such extreme hydrodynamic conditions that the Caspian Sea is located, being the largest enclosed water body on Earth and serving as a representative case for studying such processes. Its level is subject to significant fluctuations driven by the hydrological regime of the Volga River basin, atmospheric processes over the sea surface, and climate change (Akbari et al., 2020; Duan et al., 2025). Since the late 20th century, a persistent declining trend in the sea level has been observed, which over the last decade has led to large-scale restructuring of the coastal zone (Chen et al., 2023). These processes are most pronounced in the northeastern part of the Caspian Sea, where the gentle bottom topography and shallow waters cause even minor decreases in water level to result in horizontal shoreline shifts of several kilometers (Akhmetov and Malakhov, 2026). This region is of high ecological value and is also an area of intensive oil and gas development, which makes shoreline monitoring critical for assessing environmental risks (Kamza et al., 2023).

Issues of Caspian Sea shoreline dynamics have been widely studied using satellite-based methods, primarily based on medium-resolution Landsat and Sentinel-2 data (Toorani et al., 2021; Safarov and Sertel, 2023; Safarov et al., 2024; Shamsaie and Ghaderi, 2025). In analyzing the conditions of complex coastal morphology and shallow waters in the Northern Caspian, it has been established that among spectral indices the modified normalized difference water index (MNDWI) demonstrates the highest effectiveness, providing more reliable water surface delineation through suppression of signals from land and anthropogenic features (Duan et al., 2025). This result is consistent with the findings of Meirambek et al. (Meirambek et al., 2024), where comparative testing of different indices (NDWI, NDMI, MNDWI, WRI, NDVI) demonstrated the superiority of MNDWI for water body extraction tasks. The application of such satellite image processing algorithms for the northeastern sector of the Caspian Sea enables not only the mapping of the current waterline position but also the prediction of potential flooding and desiccation zones. Such detailed assessment is crucial for environmental risk analysis and sustainable management of shallow-water areas (Kamza et al., 2023).

Alongside index-based methods, object-based image analysis has also been tested for Caspian shoreline mapping, allowing minimization of spectral mixing errors through the incorporation of object texture features (Rasuly et al., 2010). However, the development of cloud-based platforms (GEE) has enabled a return to index-based approaches at a new technological level, providing high accuracy through automated processing of long-term data series.

A critical limitation of most existing shoreline extraction approaches based on medium-resolution imagery remains the spectral mixing of pixels under conditions of extremely gentle bottom slopes and extensive silted shallow-water areas of the Caspian Sea. Under such extreme hydrodynamic conditions, standard threshold-based segmentation methods, when applied without external verification, inevitably generate spatial artifacts. To overcome this limitation, the present study implements an end-to-end verification framework, the key feature of which is the integration of medium-resolution satellite data with very high-resolution PlanetScope optical imagery (3 m spatial resolution) combined with comprehensive expert visual interpretation. The use of PlanetScope-derived reference shorelines as a calibration benchmark enabled rigorous validation of automated extraction methods and ensured the identification of the optimal shoreline extraction approach.

The methodological basis for the use of PlanetScope data was established in our previous study (Dedova et al., 2026), in which shoreline dynamics in the northeastern part of the Caspian Sea were investigated using high-spatial-resolution imagery for the period 2017–2024. In the present study, the obtained high-resolution shoreline contours are used as a reference for evaluating shoreline extraction methods based on medium-resolution data. At the same time, the focus is shifted from detailed mapping of a specific area to the development and evaluation of an approach suitable for long-term monitoring of the entire Caspian Sea.

The present study aims to provide a high-resolution quantitative assessment of the spatiotemporal dynamics of the Caspian Sea shoreline and water surface area over a 10-year period from 2016 to 2025. The research is based on the concept of integrating heterogeneous Earth observation datasets, including long-term archives of Sentinel-2 and Landsat 8/9 imagery within the framework of NASA’s Harmonized Landsat and Sentinel-2 (HLS) product, very high-resolution PlanetScope imagery, and continuous satellite altimetry records from the Database for Hydrological Time Series of Inland Waters (DAHITI) database.

To achieve this objective, a series of research tasks is addressed through a three-phase methodology. During the first phase, parallel data preparation is performed, where high-resolution PlanetScope imagery is used as the reference dataset and HLS data serve as the working dataset, with the formation of temporally consistent observation samples. During the second phase, the optimal shoreline extraction algorithm is selected through a validation procedure based on the comparison of twelve shoreline scenarios with the reference data using the Digital Shoreline Analysis System (DSAS) and an integrated accuracy assessment framework. During the third phase, the selected solution is methodologically scaled to the entire Caspian Sea basin. At this stage, satellite altimetry data are used to identify periods corresponding to annual sea-level maxima, which serve as the basis for defining temporally consistent windows for optical observations.

The scientific novelty of this study lies in the integration, within a unified methodology, of validation-based selection of a shoreline extraction algorithm using high-spatial-resolution data, the generation of hydrologically consistent temporal composites based on satellite altimetry, and the subsequent scaling of the developed approach to the entire Caspian Sea for the construction of a multi-year shoreline dataset based on medium-resolution satellite data. The practical significance of the study lies in the development of a reproducible methodological framework for regular satellite monitoring of the Caspian Sea coastal zone. The proposed methodology can be used to assess changes in the coastal zone and the condition of coastal ecosystems, support environmental monitoring and management activities, and assist in the development of adaptation strategies for the Caspian region countries in response to the ongoing decline in the Caspian Sea level.

Study area

The Caspian Sea is the largest enclosed water body on Earth, located at the boundary between Europe and Asia (Figure 1).

It borders Russia, Kazakhstan, Turkmenistan, Iran, and Azerbaijan. Its surface area is approximately 371,000 km 2, the water volume is about 78,000 km 3, and the maximum depth reaches 1,025 m. According to data from JSC “NC “Kazakhstan Gharysh Sapary” for 2023, the area of the Kazakh sector of the Caspian Sea was 105,745.23 km 2 (GHARYSH, 2023). Based on its physico-geographical characteristics, the sea is conventionally divided into the Northern, Middle, and Southern Caspian. It is formed within a deep tectonic depression and contains numerous bays and islands. Significant oil and gas fields are located in this region.

The Caspian Sea is situated within several climatic zones. The northern part is dominated by a temperate continental climate, the western coast experiences a moderately warm climate, the southwestern and southern parts are characterized by a subtropical climate, while the eastern coast is influenced by a desert climate. In winter, weather conditions in the northern and central regions are shaped by Arctic air masses, whereas the Southern Caspian is influenced by southern cyclones. In summer, the influence of the Azores High ridge over the western and northwestern regions leads to stable warm conditions, while the southeastern part is affected by the Iran–Afghanistan thermal low. Prevailing winds are mainly from the north and northwest, with an average wind speed of about 6 m/s, reaching 8–9 m/s on the Absheron Peninsula. In January, temperatures in the northeastern regions may drop to −30 °C, while in the south they range from 8 to 12 °C. In summer, air temperatures reach up to 44 °C along the eastern coast. Precipitation is highly unevenly distributed: in the east it is less than 100 mm per year, along the western coast up to 1,700 mm per year, and about 200 mm over the open sea. The components of the Caspian Sea water balance include river inflow, atmospheric precipitation, groundwater inflow, evaporation, and outflow to the Kara-Bogaz-Gol Bay.

Instrumental observations of the Caspian Sea level, available since the early 20th century, indicate significant variability in its long-term level regime, which can be divided into four key stages over the century-long period (Figure 2).

From 1900 to 1929, a persistently high sea level of about −26 m was observed. This was followed by a continuous decline from 1930 to 1977, reaching an absolute historical minimum of −29.01 m. From 1978 to 1995, an intense transgressive phase was recorded, with a maximum level of −26.49 m. The modern period (1996–2023) is characterized by a renewed deep regression, during which the sea level dropped to −28.74 m by the end of 2023. The observed acceleration in the rate of decline during 2005–2023 (16.2 cm/year) indicates a substantial disruption in the outflow component of the water balance under the influence of contemporary climate trends. Similar periods of rise and fall in the Caspian Sea level and their possible drivers are discussed in detail in (Ivkina et al., 2015; Vyruchalkina et al., 2020; Safarov et al., 2024; Kostianoy et al., 2025).

The contemporary decline in the Caspian Sea level is also confirmed by satellite altimetry data. According to the Database for Hydrological Time Series of Inland Waters (DAHITI), developed by the German Geodetic Research Institute (DGFI-TUM) (Schwatke et al., 2015), the water level of the Caspian Sea (Caspian Sea, 2026) is characterized by pronounced interannual variability and an overall downward trend since the mid-1990s. The DAHITI database provides a continuous and homogeneous time series of the Caspian Sea level from 1992 to 2026, with an accuracy on the order of a few centimeters and an average temporal resolution of approximately 10 days, determined by the repeat cycle of the Jason/Sentinel-6 satellite missions. The time series shown in Figure 3 exhibits a distinct sawtooth pattern, reflecting regular cyclical fluctuations in sea level throughout the hydrological year.

The spring–summer rise in the water surface is associated with a sharp increase in surface inflow to the Caspian basin, the largest contribution of which comes from the flood runoff of the Volga River. In contrast, the autumn–winter decline in sea level is caused by the seasonal reduction in river inflow during the low-flow period and by intensified evaporation from the sea surface.

At the beginning of the observation period (December 1992), the mean monthly level of the Caspian Sea was −26.54 m. Following a local maximum in 1995 (approximately −26.10 m), a period of relative stabilization occurred during 1996–2005, with seasonal fluctuations ranging between 30 and 40 cm. This phase was replaced in 2006 by a persistent downward trend, which became particularly pronounced after 2019, when the Caspian Sea level declined to −29.03 m by early 2025. The decreasing trend continued thereafter, and by early 2026 the sea level had reached its absolute minimum for the entire period of satellite observations, dropping to −29.11 m.

Under conditions of such instability and the morphological complexity of the coastline, it is essential to have a reliable tool for water surface identification based on spectral indices. Consequently, there is a need to develop an optimized remote sensing algorithm capable of rapidly and, most importantly, geometrically accurately delineating the water surface for continuous monitoring of shoreline dynamics.

Technological workflow structure

To address the research objectives, an end-to-end methodological framework was developed to describe the separation of spatial data streams and the structure of their integrated processing (Figure 4). In this scheme, blue rectangles represent data streams, whereas yellow rectangles indicate processing steps. In a generalized form, the methodological framework was based on the calibration–validation principle and included three key phases: parallel data preparation, validation-based selection of the optimal algorithm, and scaling of the selected solution to the entire Caspian Sea basin. The first phase of the study was implemented within a general data preparation block divided into two independent parallel information streams.

The first stream was aimed at generating a reference dataset of shorelines for the northeastern Caspian Sea using summer PlanetScope composites for 2022–2024 (Figure 4, Stream 1). A separate reference shoreline was produced for each year and subsequently used for comparison with the corresponding set of shorelines extracted from HLS data, thereby ensuring temporal consistency between the datasets being compared. The 3-year observation period was selected to provide representative hydrometeorological conditions while maintaining the comparability of the input data.

For these images, spectral indices (NDWI and NDVI) were calculated using the SNAP software package, after which the results were manually refined in ArcGIS. The outcome of this stage was the creation of highly accurate reference shorelines, which were subsequently used as the basis for assessing the quality of shorelines generated by semi-automated methods. In parallel, a working dataset was created based on the HLS products, which integrate observations from the Sentinel-2 and Landsat-8/9 satellites (Figure 4, Stream 2). To ensure data comparability, scenes with acquisition dates as close as possible to those of the corresponding PlanetScope images were selected. Using these data, spectral indices were calculated and temporal composites were generated within the Google Earth Engine environment, resulting in twelve shoreline extraction scenarios (a detailed description and structure of these scenarios are provided in Section 2.4).

The second phase (Optimal Algorithm Selection) consisted of comparing these twelve shoreline scenarios with the reference shoreline using the DSAS V6 module through transect-based analysis, which enabled the determination of spatial deviations along a large number of cross-shore profiles. The validation results were then integrated using an integral ranking assessment performed in Excel. Based on the combined accuracy metrics, the most effective shoreline extraction algorithm was identified.

In the final phase, the selected optimal solution was scaled to the entire Caspian Sea to generate a continuous multi-temporal record of shoreline positions. The final product of the proposed approach was the creation of a vector shoreline database covering a 10-year period (2016–2025), validated using high-resolution reference sites.

Creation of the reference dataset

In accordance with the proposed end-to-end processing architecture, the implementation of the first information stream (Stream 1) consisted of generating a verification dataset based on imagery from the PlanetScope satellite constellation with a spatial resolution of 3 m, required for the creation of the baseline shoreline. High-resolution baseline shorelines of the northeastern Caspian Sea were generated from PlanetScope data using a methodological approach described in greater detail in our previous study (Dedova et al., 2026), which analyzed shoreline dynamics in the northeastern Caspian Sea over the period 2017–2024. In the present study, PlanetScope data for the shorter period of 2022–2024 are used as a reference dataset for the quantitative validation of shoreline extraction methods applied to HLS imagery. Considering the geomorphological characteristics of the northeastern Caspian Sea coast, as well as the maximum scene area limitations imposed by Planet Explorer for image ordering and downloading, the study area was preliminarily divided into eight local sectors. To minimize the influence of cloud cover, atmospheric disturbances, and short-term wind-induced set-up and set-down events, summer composites were generated for each year during the 2022–2024 period. These composites were produced by aggregating scenes acquired within 2–3 consecutive days for each of the eight local sectors.

The initial stage of shoreline extraction involved the calculation of spectral indices. The Normalized Difference Water Index (NDWI) (NDWI, 2026) was used as the primary tool for water surface delineation, whereas the Normalized Difference Vegetation Index (NDVI) (EOS, 2026) served as an auxiliary criterion for verifying the water–land boundary. The resulting thematic index layers were subjected to threshold classification with an increment of 0.05 within the ranges of 0.2–0.3 for NDWI and −0.3 to −0.2 for NDVI. For each selected threshold value, binary rasters were generated and subsequently converted into polygons and then into vector lines. From the resulting set of alternative contour lines, a single representative shoreline was selected through visual comparison with the original satellite imagery. The selected line was the one that most accurately corresponded to the actual waterline position. Such an approach is widely applied in contemporary coastal monitoring practices to minimize spectral classification errors under complex coastal conditions (Alcaras et al., 2022; Usman et al., 2025).

The selected shorelines were then edited in the ArcGIS environment. Manual post-processing included vertex optimization using the Douglas–Peucker algorithm (Douglas and Peucker, 1973) with a tolerance of 10 m to eliminate the discretization effect arising from the vectorization of raster cells. At this stage, topological errors were corrected, vector continuity was verified, and spatial artifacts were removed. Based on the corrected linear features, final polygonal objects were generated, including the delineation of islands, with mandatory verification of topological consistency. The final output of this processing workflow was the creation of a polygon layer representing the water surface of the Northeastern Caspian Sea and a reference shoreline dataset. The individual processing steps are schematically illustrated in Figure 5.

Despite their high spatial resolution (3 m), the reference shorelines possess intrinsic positional uncertainty associated with image georeferencing accuracy, the spatial resolution of the source data, and the manual vectorization procedure. As a conservative estimate of the reference uncertainty associated with the spatial resolution of the source data, a value corresponding to one PlanetScope pixel (3 m) was adopted. This value should not be interpreted as an estimate of the total positional uncertainty of the reference shoreline. Due to the absence of independent high-precision geodetic, aerial survey, or LiDAR data, the total positional uncertainty of the reference shoreline was not quantitatively assessed. Since the spatial resolution of PlanetScope is approximately an order of magnitude higher than that of the HLS imagery (30 m), the influence of the uncertainty of the reference data on the comparative accuracy assessment is considered limited.

Creation of the working dataset

The working dataset was created using version 2.0 of NASA’s Harmonized Landsat and Sentinel-2 (HLS) products (HLS, 2026) within the cloud computing environment of Google Earth Engine (GEE) (Gorelick et al., 2017). The primary source of satellite data was the harmonized HLS-S30 product (Sentinel-2, 30 m). To fill spatial gaps caused by insufficient coverage of the study area, Sentinel-2 Harmonized data, resampled to a spatial resolution of 30 m, were additionally used, along with HLS-L30 data (Landsat 8/9). The main processing stages are presented in Figure 6. In the figure, blue rectangles represent data, while yellow rectangles indicate processing steps. The HLS data used in this study consist of Bottom-of-Atmosphere (BOA) reflectance values that have undergone atmospheric correction, geometric co-registration, and resampling to a common spatial resolution of 30 m and the UTM Zone 39N (WGS84) map projection. After importing the data into GEE, image scenes were selected according to their temporal proximity to the corresponding reference PlanetScope images (±7 days), the percentage of cloud cover (<40%), and their inclusion within the study sector. To remove clouds and cloud shadows, the Fmask quality mask, included as part of the HLS product, was applied. This mask enables the automatic identification and exclusion of pixels unsuitable for analysis based on the interpretation of the corresponding quality-assessment bit flags. For processing the original Sentinel-2 data, cloud masking was performed using the QA60 band.

The next processing stage involved the calculation of several spectral indices. In particular, the Normalized Difference Vegetation Index (NDVI), the Modified Normalized Difference Water Index (MNDWI) (Xu, 2006), and the Automated Water Extraction Indices (AWEI) were used, including the version designed for shadow-free conditions (AWEI1) and the version accounting for shadow effects (AWEI2) (Feyisa et al., 2014). Unlike the original NDWI, MNDWI incorporates the shortwave infrared (SWIR) band, which enables improved discrimination between water and land surfaces. The AWEI modifications, including AWEI1 and AWEI2, utilize multispectral information to enhance the reliability of water body extraction under different imaging conditions, with AWEI1 primarily intended for water detection and AWEI2 effectively suppressing the influence of cloud cover and terrain shadows.

The subsequent processing stage involved the temporal aggregation of satellite data. To reduce the effects of cloud contamination, atmospheric noise, and random outliers, temporal image composites were generated, including median values, mean values, and the 25th percentile.

To determine the optimal threshold separating the water and land classes, the automatic thresholding algorithm Otsu’s method was applied. This algorithm analyzes the histogram distribution of spectral index values within the study area and identifies the threshold value that maximizes the between-class variance of the two groups of pixels. The algorithm was applied sequentially to all calculated spectral indices, temporal composite variants, and individual sectors, enabling the automatic generation of image-specific classification thresholds for each image.

At the next stage, post-processing and morphological cleaning of the binary “water–land” mask were performed to remove noise and improve its topological consistency using connected-component filtering. Spatial thresholds of 1,500 and 500 m were specified for the removal of inland water bodies and small islands, respectively, corresponding to 50 and 17 pixels at a spatial resolution of 30 m.

The resulting binary raster of the water surface was converted into an ESRI Shapefile polygon vector format with a spatial resolution of 30 m. Shoreline extraction was performed using a polygon boundary extraction algorithm followed by geometry simplification (Simplify) to remove redundant vertices.

As a result, a working dataset consisting of 12 shoreline variants was generated for each analyzed site. These shoreline datasets were subsequently used in the transect analysis within the Digital Shoreline Analysis System (DSAS) software to quantitatively assess their spatial displacement relative to the reference shoreline.

Multi-criteria assessment and selection of the optimal shoreline extraction method

To determine the most robust and accurate method for shoreline delineation from HLS data, a validation procedure was performed using a reference shoreline digitized from high-resolution PlanetScope imagery. The quantitative assessment of the agreement between the working shorelines and the reference shoreline was carried out using the transect method in the stand-alone desktop application DSAS v6.1.177 (Digital Shoreline Analysis System), developed by the U.S. Geological Survey (USGS) (Henderson et al., 2026; Woods Hole Coastal and Marine Science Center, 2026).

Unlike the classical approach for analyzing long-term shoreline dynamics, in the present study the baseline was generated directly from the reference shoreline and represented the reference position of the water–land boundary, against which alternative shoreline variants from the working dataset, derived from HLS data, were evaluated.

For the northeastern coast of the Caspian Sea, the study area was preliminarily divided into eight sectors (Figure 7). For each year and each sector, a system of transects oriented perpendicular to the corresponding reference shoreline was generated using DSAS. Transects were generated with a nominal spacing of 5,000 m. Since the reference shoreline underwent virtually no generalization and preserved the natural sinuosity of the coastline, the transects were distributed along its actual geometry rather than along a simplified shoreline. As a result, in areas with complex coastal morphology, the effective spacing between adjacent transects became smaller than 5,000 m, providing a denser sampling of shoreline position. Depending on the morphological characteristics of each coastal sector, transect lengths varied from 1 to 10 km. Due to the complex geometry of the reference shoreline and the presence of small-scale features, the automatically generated transects did not always satisfy the perpendicularity condition. In such cases, manual corrections were performed, including adjustments to the position, orientation, and spatial spacing of the transects according to the local shoreline geometry and, when necessary, modifications of the transect spacing.

For each transect, the intersection point with the corresponding shoreline from the working dataset of test shorelines, generated for the same time window as the reference shoreline, was determined. Subsequently, the cross-shore distance along the transect from the reference shoreline to the shoreline derived from the HLS dataset was measured. In cases of multiple intersections, which are typical for lagoonal and marshy areas, the seaward-most intersection was retained. This approach ensured methodological consistency of the measurements and excluded the influence of inland water bodies and temporary water patches. Figure 8 illustrates an example of the DSAS analysis performed for a shoreline segment in the northern Caspian Sea.

The obtained distance values were interpreted as positional deviations of the shoreline from the working dataset relative to the reference baseline: where i is the transect number, y is the year, m is the processing scenario, x i m y is the shoreline position obtained using method m, and x i r e f y is the position of the reference shoreline.

For each processing scenario, aggregated deviation statistics were calculated for all transects, including the Median Absolute Error (MdAE), Mean Absolute Error (MAE), Root Mean Square Error (RMSE), and the 90th percentile of the absolute error (P90). The optimal shoreline extraction scenario was selected based on a comparative analysis of these metrics.

For each of the eight sectors, transect systems were generated and comparisons of shorelines from the working dataset were performed for three time periods (2022–2024). In total, 288 shoreline scenarios were analyzed and more than 21,000 distance measurements between shorelines were performed, providing a statistically robust assessment of the accuracy of different shoreline extraction methods.

To quantitatively evaluate shoreline extraction accuracy, an integral ranking score was employed. The score is based on a weighted sum of the ranks of four evaluation criteria and provides a comprehensive assessment of the performance of different satellite data processing scenarios. The integral ranking score R int was calculated as follows: where R a n k m e d i a n, R a n k R M S E, R a n k P 90 and R a n k C o v e r a g e denote the ranks assigned to each method based on MdAE, RMSE, P90 and the transect coverage metric, respectively. The transect coverage metric was calculated as follows: where N t o t a l is the total number of transects, and N u s e d is the number of transects included in the calculations for a given study sector.

For each study sector and each year, all methods were ranked independently. The best value for each evaluation criterion was assigned a rank of 1, while subsequent values received progressively higher ranks. In cases where two or more methods produced identical values, they were assigned the same rank, and the subsequent ranks were allocated according to the standard competition ranking procedure. For the MdAE, RMSE, and P90 metrics, lower values were assigned better (lower numerical) ranks, whereas for the Coverage metric, higher values (i.e., a greater proportion of valid transect intersections) were assigned better ranks.

The weighting coefficients were selected using an expert-based approach, taking into account the informativeness of the evaluation metrics and the objectives of the study. The highest weight (0.4) was assigned to the median absolute error, as it best characterizes the typical positional accuracy of the extracted shoreline while remaining robust to individual outliers. The root mean square error was assigned a weight of 0.3 because of its sensitivity to large local deviations. The 90th percentile of the absolute error was assigned a weight of 0.2 to account for the robustness of each method under the most challenging conditions. The lowest weight (0.1) was assigned to the Coverage metric, which was considered a supplementary criterion reflecting the completeness of shoreline extraction, expressed as the proportion of transects intersecting the reference shoreline. To evaluate the sensitivity of the integral ranking score to the selected weighting scheme, two alternative weighting configurations were also evaluated: (0.3, 0.3, 0.3, 0.1) and (0.25, 0.25, 0.25, 0.25).

Scaling the methodology and creating a multi-year geospatial dataset of the Caspian Sea water surface

The scaling of the methodology to the entire Caspian Sea was implemented through the semi-automated generation of a multi-year dataset of shorelines and water surface polygons within the GEE environment. The primary source of satellite data was the Harmonized Landsat and Sentinel-2 (HLS-S30) product. To eliminate spatial gaps in the coverage of the study area, additional Sentinel-2 Harmonized data, resampled to a spatial resolution of 30 m, were used. For the period 2016–2019, HLS-L30 data (Landsat 8/9) were additionally used.

To minimize the influence of seasonal variability and ensure the correct analysis of interannual observations, the generation of the shoreline and water surface polygon time series was performed for periods corresponding to, or close to, the annual maximum sea level. For the Caspian Sea, such conditions generally occur during the summer season. An additional advantage of this period is the relatively low cloud cover, which increases the availability of high-quality optical satellite imagery and facilitates more reliable delineation of the water–land boundary.

To determine the dates of annual sea-level maxima, satellite altimetry data from the DAHITI database, which provides long-term time series of Caspian Sea water levels [54], were used. Based on these data, the dates corresponding to the annual maximum sea levels were identified for each year, and the results are presented in Figure 9. Analysis of the time series showed that the overwhelming majority of annual maxima (31 out of 33 years) occurred in July (18 years) and June (13 years). The rare exceptions (May 2006 and May 2019) were most likely associated with an early spring flood period.

The selection of the 2016–2025 period for constructing the multi-year shoreline and water surface polygon dataset was determined by the operational deployment of the Sentinel-2 satellite constellation, which, together with Landsat 8/9 observations, substantially increased the density of satellite acquisitions. This provided a sufficient number of cloud-free images to achieve complete coverage of the Caspian Sea and to generate a consistent multi-year dataset. The use of earlier periods is limited by insufficient observation frequency and incomplete spatial coverage, which complicate the construction of a homogeneous time series for the entire water body.

The satellite data processing was performed in the Google Earth Engine (GEE) environment in accordance with the developed methodology, presented in the form of pseudocode (Algorithm 1). The complete Google Earth Engine code implementing the described algorithm is provided in the Supplementary Material (CaspS_shoreline.js). For each year, satellite images were selected within an individual temporal window linked to the date of the annual maximum sea level. Depending on the year, the temporal window included 10–25 days before and after the maximum level date. The extension of the temporal window and its asymmetric positioning relative to the maximum date made it possible to increase the number of available cloud-free observations and to form a more robust set of satellite data for subsequent shoreline delineation. For the selected images, the spectral indices MNDWI and AWEI2 were calculated, a spatially continuous MNDWI composite based on the 25th percentile was generated, water surface classification was performed using the Otsu method, and the resulting mask was post-processed by removing small inland water bodies and islands. The resulting output was a raster water–land classification mask with a spatial resolution of 30 m. The code also provides the option to convert the resulting raster mask into vector format and export the results as a Shapefile.

Due to limitations in computational resources and available memory when processing large areas in the GEE environment, the Caspian Sea water area was divided into three spatial sections. Each section was processed through a separate execution of the algorithm, followed by the merging of the obtained results in ArcGIS.

The calculation of the median AWEI2 composite (AWEI2_median) can be performed through a separate execution of the code by modifying the corresponding parameters for generating the output.

Input:

Area of Interest (AOI), spatial subregions of the Caspian Sea; date of the annual maximum sea level; annual observation time window; maximum allowable cloud cover; minimum size of inland water bodies to be removed; minimum size of islands to be removed.

Output:

A raster mask of the shoreline/water–land boundary with a spatial resolution of 30 m and, if required, vector shoreline polygons/lines.

1. Formation of input parameters and data

a. Divide the Caspian Sea water area into spatial subregions and select the subregion for processing.

b. Define the AOI for Otsu threshold calculation.

c. Specify the reference date corresponding to the annual maximum sea level and define the annual temporal window.

d. Set the allowable cloud cover.

e. Retrieve HLS-S30 and Sentinel-2 Harmonized data covering the AOI and the specified temporal window.

f. Apply spatiotemporal filtering and exclude specified orbits and tiles with coverage gaps or artifacts.

g. Generate a set of satellite scenes suitable for processing.

2. Preprocessing of satellite images

a. Apply an Fmask-based quality mask to the HLS-S30 images.

b. Apply the QA60 cloud mask to the Sentinel-2 images.

c. Apply scaling factors to the optical spectral bands to convert the values to reflectance.

d. Harmonize the spectral bands across the data sources used.

3. Calculation of water indices

a. Calculate the Modified Normalized Difference Water Index (MNDWI): MNDWI=(Green−SWIR1)/(Green + SWIR1)

b. Calculate the Automated Water Extraction Index with shadow suppression (AWEI2): AWEI2 = Blue+2.5×Green−1.5×(NIR + SWIR1)−0.25×SWIR2

4. Generation of the multi-temporal 25th-percentile MNDWI composite

a. Merge the HLS-S30 and Sentinel-2 Harmonized collections and identify unique observation dates.

b. Generate a daily mosaic of the HLS-S30 and Sentinel-2 Harmonized collections for each unique date.

c. Generate a spatially continuous 25th-percentile MNDWI composite for the selected study region.

Algorithm 1 Pseudocode of the Google Earth Engine workflow for annual shoreline extraction. 5. Determination of the water–land separation threshold using the Otsu method

a. Calculate the histogram of MNDWI values within the specified AOI.

b. Apply the Otsu automatic thresholding method to the resulting histogram.

c. Classify pixels with MNDWI values above the Otsu threshold as water surface.

6. Post-processing and morphological filtering

a. Remove inland water bodies smaller than the specified threshold.

b. Remove islands smaller than the specified threshold.

7. Export of results

a. Export the final raster mask of the shoreline/water–land boundary at a spatial resolution of 30 m.

b. If required, convert the raster boundary into vector polygons and export them in Shapefile format.

Results of reference shoreline generation from PlanetScope data

The implementation of the first information stream resulted in the creation of a reference dataset of shorelines and water surface polygons for the northeastern Caspian Sea based on PlanetScope satellite imagery with a spatial resolution of 3 m. For each year of the study period (2022–2024), a separate reference shoreline was generated and subsequently used for the validation of automated water–land boundary extraction methods based on HLS data. Figure 7 illustrates the spatial distribution of the reference shorelines for the three observation periods.

The figure also shows the division of the study area into eight local sectors adopted for the subsequent analysis. The northeastern Caspian Sea is characterized by a highly dynamic coastal zone due to the shallow-water nature of the basin and the high sensitivity of shoreline position to sea-level fluctuations.

A comparison of the shorelines for 2022–2024 indicates the continued retreat of the sea and a reduction in the water surface area. The most pronounced changes were observed in the northern and northeastern parts of the study area, where shoreline displacement reached up to 10 km in some locations. During the study period, the area of this part of the water body decreased by nearly 2,425 km 2. The largest reduction occurred between 2023 and 2024 (1,497.84 km 2), compared to the period between 2022 and 2023 (927.14 km 2). At the same time, the identified changes exhibited considerable spatial heterogeneity. In some sectors, the 2023 shoreline was displaced seaward relative to the 2022 shoreline, while in 2024, despite the overall trend of declining sea level, local areas exhibited the opposite pattern (Dedova et al., 2026). These features are most likely associated with the influence of short-term hydrometeorological factors, primarily wind-induced set-up and set-down events, which are characteristic of the shallow northeastern Caspian Sea. Following generalization and post-processing procedures, the total length of the reference shoreline was estimated at 3,345 km in 2022, 3,580 km in 2023, and 3,520 km in 2024. The variability of shoreline length compared with changes in water surface area is related to the sensitivity of linear morphometric parameters to local shoreline geometry. Under shallow-water coastal conditions, even minor sea-level fluctuations can lead to shoreline reconfiguration and the formation of temporary geomorphological features, thereby increasing shoreline variability.

Results of generating the working shoreline dataset from HLS data

The implementation of the second information stream resulted in the creation of a working shoreline dataset for each of the eight study sectors and each temporal snapshot (2022–2024) based on HLS data with a spatial resolution of 30 m. The dataset included 12 alternative scenarios for automated water–land boundary extraction.

The most challenging areas for automated shoreline delineation were Sectors 4 and 5, located along the eastern coast of the northern Caspian Sea (the sector locations are shown in Figure 7). These areas are strongly affected by wind-induced set-up and set-down events associated with the impact of wind forcing on the shallow-water basin. As a result, even when satellite images were acquired on closely spaced dates, the position of the waterline could vary substantially, leading to increased spatial discrepancies among shorelines derived using different extraction methods. As an example, Figure 10 presents the results of automated shoreline extraction for Sector 4 located on the eastern coast of the Northern Caspian.

Additional difficulties were identified in Sectors 1 and 3, situated near the Volga Delta and the Ural River Delta, respectively. Extensive reed beds and emergent coastal vegetation complicated the spectral separation of the water and land classes due to the presence of mixed pixels and the partial overlap of spectral signatures between water surfaces and vegetation cover. Under such conditions, some spectral indices exhibited instability in determining shoreline position.

It was found that the influence of the temporal compositing approach on shoreline extraction results depended on the number of available satellite observations within the selected temporal window. In cases where only a single cloud-free image was available for a particular sector and year, the application of different statistical aggregation methods (mean, median, and p25) had no effect on the resulting shoreline configuration, since the composite effectively corresponded to the original image. Despite the occurrence of such cases, multi-scene composites were available for the majority of the analyzed combinations, allowing the influence of different temporal aggregation approaches on water surface extraction accuracy to be evaluated.

The resulting working shoreline dataset encompassed a wide range of hydrometeorological and morphological conditions characteristic of the northeastern Caspian Sea and provided the basis for the subsequent quantitative assessment of shoreline extraction algorithms using PlanetScope reference data and transect-based analysis in the DSAS software package.

Results of selecting the optimal shoreline extraction method

To identify the most effective method for automated shoreline extraction, a comparative analysis was conducted between the 12 HLS-based shoreline extraction scenarios from the working dataset and the reference shoreline.

The preliminary analysis showed that although methods based on the NDVI index demonstrated superior performance for certain individual metrics in some coastal sectors, substantial errors were observed in other areas, including omissions of shoreline segments. This indicates the insufficient robustness of this spectral index for automated shoreline delineation. Therefore, the NDVI-based methods were excluded from the final comparative analysis, which was subsequently performed only for the specialized water indices MNDWI, AWEI1, and AWEI2.

The integral ranking scores calculated for each method, study sector, and year using the adopted weighting scheme (0.4, 0.3, 0.2, and 0.1) are presented in Table 1.

Based on the results presented in Tables 1, 2 summarizes the methods that achieved the lowest integral ranking score for each study sector and year. Where multiple methods obtained the same lowest score, all methods sharing the first rank were included.

Tables 1, 2 show that the lowest integral ranking scores were most frequently achieved by the MNDWI_p25 method. This method performed best for the majority of the investigated study sectors and years, particularly in sectors 1, 2, 5, and 6. Overall, MNDWI_p25 ranked first in 11 out of 24 cases (45.8%), including cases where multiple temporal compositing strategies (mean, median, and p25) produced identical results. MNDWI_median also demonstrated high performance, ranking first in several cases or yielding integral ranking scores very close to the minimum. The AWEI2 family of methods produced competitive results only in specific coastal sectors. In particular, AWEI2_mean and AWEI2_median achieved the lowest integral ranking scores in sector 5 (2023) and sector 7 (2022), while in sector 6 (2022) all AWEI2 variants shared the first rank. By contrast, the AWEI1 family generally exhibited higher integral ranking scores, indicating lower shoreline extraction accuracy than both MNDWI and AWEI2. When considering the spectral index families as a whole, the MNDWI variants (mean, median, and p25) ranked first in 17 out of 24 cases (70.8%), substantially outperforming AWEI2 (4 out of 24 cases, 16.7%) and AWEI1 (3 out of 24 cases, 12.5%). Analysis of the temporal compositing strategies further showed that the 25th-percentile composite (p25) most frequently provided the best performance, followed by the median composite, whereas the mean composite only rarely achieved the lowest integral ranking score. In several cases, different compositing strategies produced identical integral ranking scores, indicating that the choice of temporal compositing had only a minor influence on shoreline extraction accuracy under those specific conditions.

In addition to the adopted weighting scheme (0.4, 0.3, 0.2, 0.1), two alternative weighting configurations, (0.3, 0.3, 0.3, 0.1) and (0.25, 0.25, 0.25, 0.25), were also evaluated. In all cases, the overall conclusion regarding the superiority of the MNDWI_p25 method remained unchanged. Changes were observed only in a few study sectors where the values of the evaluated accuracy metrics were nearly identical. Under the weighting scheme (0.3, 0.3, 0.3, 0.1), the highest-ranked method changed only in sector 8 (2023), where AWEI2_p25 slightly outperformed MNDWI_p25, while in sector 3 (2022) MNDWI_mean and AWEI2_p25 received identical integral ranking scores. When equal weights (0.25, 0.25, 0.25, 0.25) were applied, an additional change occurred in sector 1 (2023), where AWEI2_p25 ranked first, whereas the changes observed in sectors 8 (2023) and 3 (2022) remained unchanged. Thus, changes in the weighting coefficients did not result in substantial changes in the overall ranking of the methods: the changes were limited to individual cases characterized by nearly identical values of the integrated score. Across all three weighting schemes, MNDWI_p25 retained its overall advantage, confirming the robustness of its selection as the primary method.

The obtained results reflect the evaluation of the methods solely in terms of the positional accuracy of the extracted shoreline. However, subsequent analysis demonstrated that the quality of the derived water masks may vary considerably among methods. In particular, the use of MNDWI_p25 composites often resulted in a more fragmented water surface with local discontinuities, likely caused by wave dynamics and variations in spectral reflectance along the water–land interface. In contrast, the AWEI2_median method generated a more spatially connected and homogeneous water mask without pronounced internal artifacts that could otherwise be interpreted as false island features.

The calculated integral ranks and related tabulated results are provided as Supplementary Material (CaspS_shoreline_rank.xlsx).

Therefore, considering both the positional accuracy assessment and the analysis of water surface spatial integrity, MNDWI_p25 was selected as the primary method for shoreline delineation due to its higher sensitivity to the water–land boundary. For the selected MNDWI_p25 method, the threshold values automatically determined using the Otsu method ranged from 0.045 to 0.056. AWEI2_median was used for the visual inspection of the water mask obtained using MNDWI_p25. Local gaps and “island” fragments within the water surface identified in the MNDWI_p25 mask were compared with the AWEI2_median mask to identify classification errors. The identified MNDWI_p25 classification errors were then corrected during manual vector editing in the ArcGIS environment. Thus, AWEI2_median was not used to directly determine the shoreline position, but served as an additional quality-control layer for detecting errors in the MNDWI_p25 water mask prior to manual polygon correction.

Results of generating a multi-year geospatial dataset of the Caspian Sea water surface

Based on the selected semi-automated water surface extraction method using Sentinel-2 and Landsat 8/9 data, a multi-year dataset of Caspian Sea water surface polygons was generated for the period 2016–2025. The resulting annual Caspian Sea water surface datasets for 2016–2025 are provided as Supplementary Material (Caspian_Water_Surface_2016_2025.gpkg).

An overview of the shoreline transformation is presented in Figure 11. A comparison of the water surface polygons derived for 2016 and 2025 revealed that the most significant reduction in water surface area occurred within shallow-water regions, which are highly sensitive to sea-level fluctuations.

To examine the spatial characteristics of coastal transformation in greater detail, several coastal sectors exhibiting the most pronounced changes during the study period were identified.

The most significant changes were observed in the northeastern Caspian Sea within the territory of Kazakhstan (Figure 12a). Here, the decline in sea level was accompanied by a substantial reduction in shallow-water areas, an increase in the extent of exposed land, and pronounced shoreline retreat. Along almost the entire eastern coast of the Northern Caspian, the shoreline retreated by more than 20 km during the observation period. Significant changes were also detected in the area of the Tyuleniy Islands (Figure 12b). As a result of the sea-level regression, several island sections in the eastern part of the archipelago became connected to the mainland, while the area of the remaining island portion increased by more than 2.5 times, from 200.47 km 2 to 504.98 km 2.

A marked reorganization of the coastal zone was also observed near the Khazar Peninsula, Turkmenistan (Figure 12c), where the decline in sea level resulted in shoreline displacements of several kilometers in certain areas and was accompanied by changes in the configuration of coastal promontories, spits, and shallow bays. The most intensive changes occurred between 2022 and 2023. Unlike the northeastern Caspian, where extensive low-lying territories became exposed, changes in this region were primarily expressed through the redistribution of coastal landforms and the restructuring of local shoreline morphology.

Substantial shoreline changes were also identified in the vicinity of the Volga Delta and along the northwestern Caspian coast (Figure 12d). Analysis of the extracted shorelines revealed that Chechen Island, located off the Russian coast, became connected to the mainland in 2023, effectively becoming part of the Agrakhan Peninsula.

Table 3 summarizes the changes in the Caspian Sea water surface area, shoreline length, and annual maximum sea level for the period 2016–2025, derived from satellite monitoring data.

The absolute values of the water surface area and shoreline length are determined by the spatial resolution of the input satellite data (30 m), the selected methodology for water polygon extraction, subsequent post-processing procedures, and the vector generalization rules applied. Therefore, the reported absolute values should be regarded as approximate estimates. However, the consistent application of the same processing methodology and unchanged spatial resolution throughout the 2016–2025 time series ensures the methodological consistency and reliability of the identified long-term trends.

In 2016, the Caspian Sea water surface area was 390,644.77 km 2, whereas by 2025 it decreased to 379,428.46 km 2. Thus, over the study period the water surface area declined by almost 11,216.31 km 2, while the maximum sea level decreased from −27.13 m to −28.46 m, i.e., by 1.33 m relative to the 2016 level. A 1 m decrease in the Caspian Sea level was accompanied by an approximate reduction of 8,433 km 2 in water surface area.

Figure 13 illustrates the rate of sea-level decline and the reduction in water surface area over the period 2016–2025. The results indicate that comparable changes in sea level may be accompanied by substantially different reductions in the water surface area. For example, a 0.18 m drop in sea level in 2023 was accompanied by a reduction of 2,135.19 km 2, whereas an equivalent 0.18 m decrease in 2025 resulted in a much smaller area change of 1,284.48 km 2. At the same time, the largest observed sea-level decrease of 0.32 m in 2022 corresponded to a reduction of 932.64 km 2 in water surface area. The results demonstrate that changes in the Caspian Sea water surface area are controlled not only by the magnitude of sea-level decline but also by the geomorphological characteristics of the coastal zone. This effect is particularly pronounced in the shallow northern part of the sea, where even relatively small water-level fluctuations result in substantial shoreline retreat.

The peak rate of area loss occurred during 2022–2024, when 4,359.83 km 2 of water surface was lost in total (maximum annual loss of 2,224.65 km 2 in 2024). Overall, after 2022 an acceleration in the rate of water surface reduction is observed, which led to shoreline retreat exceeding 20 km in some sections of the northeastern coast, transformation of island systems, and changes in the configuration of shallow bays. This indicates that several coastal zones are transitioning into a phase of active morphodynamic reorganization and suggests that parts of the shallow northern and northeastern Caspian Sea are entering an active desiccation stage.

An additional indicator of coastal morphological restructuring is the dynamics of shoreline length (Table 3). Over the study period, shoreline length decreased from 19,907.94 km to 15,551.79 km. At the same time, interannual variability in this parameter is higher compared to water surface area. This is related to the high morphological complexity of shallow coastal zones, where sea-level decline leads to intensive shoreline reconfiguration, including changes in the geometry of spits, lagoons, islands, and other depositional landforms.

The obtained results are consistent with previous studies addressing recent changes in the Caspian Sea. In particular, the observed sea-level decline agrees with established trends associated with changes in the water balance, increased evaporation, and climate-driven alterations in atmospheric circulation over the region (Vyruchalkina et al., 2020; Kostianoy et al., 2025). The spatial concentration of the most intensive changes in the shallow northern and northeastern Caspian Sea confirms satellite-based shoreline monitoring results reported in (Ivkina et al., 2015; Akhmetov and Malakhov, 2026), where these areas are characterized by maximum sensitivity to sea-level fluctuations.

Thus, the generated multi-year geospatial dataset enables quantitative assessment of spatial and temporal patterns of coastal transformation in the Caspian Sea and can serve as a basis for further monitoring, modeling, and forecasting of morphodynamic processes under continuing sea-level decline.

A distinctive feature of the present study is the implementation of an end-to-end methodological framework based on the validation-driven selection of an optimal algorithm followed by its methodological scaling to the entire Caspian Sea. The construction of optical image acquisition time windows was performed using satellite altimetry data, which ensured their hydrological consistency. This approach provides spatiotemporal comparability of conditions between different years and reduces the influence of seasonal variability of Caspian Sea levels on the analysis results.

The present study is based on the data and methodological approaches developed in our previous work (Dedova et al., 2026); however, its objectives and results are substantially broader. The previous study focused on a detailed analysis of the shoreline of the northeastern part of the Caspian Sea using PlanetScope data for 2017–2024. In the present study, PlanetScope is used as a reference dataset for validating shorelines derived from HLS data, while the main focus is placed on comparing 12 processing scenarios, assessing shoreline extraction accuracy using DSAS, and selecting the most optimal method. The use of the selected MNDWI method with statistical aggregation based on the lower 25th percentile enabled the generation of a consistent shoreline time series for the entire Caspian Sea for the period 2016–2025.

The selection of the test sites was based on two main criteria: the diversity of coastal conditions and the complexity of delineating the water–land boundary. The sample included deltaic and wetland areas of the Volga (sector 1) and Ural (sector 3) rivers, shallow-water coasts (sectors 4–5), sandy and sandy-muddy sections (sectors 5–6), morphologically complex and highly indented shorelines (sector 7), as well as areas with a clearly defined water–land boundary, including rocky and erosional coastal features in the Tyub-Karagan area (sector 8). In addition, areas subject to pronounced anthropogenic influence were included, including the Bautino port infrastructure (sector 8) and the dam of the Kalamkas field (sector 6). At the same time, priority in site selection was given to the most dynamic and challenging areas for automated shoreline extraction, particularly in the northeastern part of the Caspian Sea. This is because the low seabed slopes and substantial variability of the waterline in this region create the most challenging conditions for shoreline delineation. In contrast, most of the Caspian Sea coastline is characterized by a more clearly defined water–land boundary, making shoreline extraction less challenging in these areas.

The obtained results demonstrated the superiority of the MNDWI family of methods, particularly the MNDWI_p25 variant, which most frequently produced the lowest integral ranking scores. This is likely related to the combination of the characteristics of the spectral index itself and the temporal compositing approach. Unlike the linear combinations of spectral bands used in the AWEI family of indices, MNDWI is based on the normalized contrast between the green and shortwave infrared bands, providing a more robust separation of water and land. An additional advantage was provided by the use of the 25th-percentile composite, which reduces the influence of mixed pixels along the water–land boundary and minimizes the impact of short-term anomalies associated with surface wave conditions and hydrometeorological factors when determining the Caspian Sea shoreline. This is particularly important because extensive shallow waters, waterlogged soils, and frequent wind-driven set-up and set-downs events substantially complicate the automatic extraction of the water–land boundary. This is particularly important for the northeastern Caspian Sea, where wind-driven set-up and set-down events are a characteristic feature of the hydrological regime (Yeltay et al., 2026). At the same time, analysis of the spatial structure of the extracted water masks showed that the AWEI2_median method produces a more homogeneous water surface without local discontinuities. Therefore, it was used as an additional quality-control tool for evaluating the results of automated water surface extraction.

The obtained results are consistent with recent studies reporting accelerated shoreline regression and a reduction in the water surface area of the northern Caspian Sea in recent years (Court et al., 2025; Duan et al., 2025), and they also have important practical implications. In particular, the continuing reduction of the water surface area in the northern Caspian Sea leads to a decrease in the breeding habitats of the Caspian seal and the habitats of sturgeon species. Another illustrative example is the gradual expansion of the Tyuleniy Islands archipelago, which may potentially alter the habitat conditions of the endemic Caspian seal, colonies of waterbirds, and other island ecosystems (Court et al., 2025). Although the assessment of ecological consequences was beyond the scope of the present study, the proposed methodology can be used for their future monitoring.

A similar relationship between changes in water level and coastal zone transformation has also been reported for other shallow-water bodies, such as the Aral Sea, although the quantitative relationships are determined by the specific morphometric characteristics of each basin (Shibuo et al., 2006; Ivkina and Shivareva, 2010). However, the quantitative relationships between water level and water surface area are unique to each water body and cannot be directly transferred without considering its morphometric characteristics (Duan et al., 2025).

Despite the obtained stable and consistent results, the proposed methodology has several limitations.

First, the semi-automated delineation of shorelines and water polygons was performed using HLS data (and additionally Landsat 8/9 and Sentinel-2 in cases of insufficient coverage) with a spatial resolution of 30 m. This limitation leads to a partial loss of fine-scale morphological information: narrow channels, small islands, and complex shoreline features may be smoothed or not detected. Consequently, this affects the absolute estimates of water surface area and shoreline length, particularly in regions with highly fragmented coastal geometry. An additional limiting factor is the dependence of optical data on cloud cover and atmospheric conditions, which may introduce temporal sampling heterogeneity.

Second, the DSAS transect-based method is designed exclusively for assessing the positional displacement of the shoreline and does not account for the internal structure of the water mask. As a result, even when shoreline positional accuracy is high, the continuity of the water surface may not be fully represented. In addition, within the framework of the present study, transect analysis was limited to the 2022–2024 period due to the high computational and manual effort required for transect correction in areas of complex shoreline morphology. Although this period is sufficient for comparative evaluation of methods, in some cases a reduced number of shoreline–transect intersections was observed due to local displacements or missing shoreline segments, which may affect the statistical robustness of the estimates.

Third, the reference shorelines derived from PlanetScope data are not error-free and possess their own positional uncertainty, primarily associated with the 3 m spatial resolution, image geolocation accuracy, spectral classification, and subsequent manual editing. The 3 m value should therefore be regarded as a conservative estimate of the uncertainty specifically associated with the spatial resolution of PlanetScope, rather than as an estimate of the total positional uncertainty of the reference shoreline. Additional sources of uncertainty may include image geolocation errors, manual vectorization and interpretation of the water–land boundary, as well as the effects of spectral mixing in shallow-water areas, water turbidity, wind-induced surface waves, and short-term hydrometeorological changes. Owing to the absence of independent high-accuracy geodetic, aerial photography, or LiDAR data, the absolute positional uncertainty of the reference shorelines was not quantitatively assessed. However, this uncertainty is substantially smaller than the spatial resolution of the HLS imagery (30 m) used for shoreline extraction. Therefore, the reported accuracy metrics should be interpreted as a relative assessment of the performance of the evaluated methods with respect to a common reference dataset rather than as an absolute measure of shoreline positional error.

Fourth, the quantitative validation of the methods under consideration was conducted only for test sites in the northeastern part of the Caspian Sea. Although the selected sites encompass a range of substantially different coastal conditions, the validation sample does not represent all types of coastal conditions of the Caspian Sea equally. Therefore, the application of the methodology to the entire sea area should be regarded as justified methodological scaling rather than as a fully comprehensive quantitative validation for every type of coastal environment. To assess the consistency of the results beyond the validation sites, visual inspection of the extracted shorelines was additionally performed.

Fifth, the use of time windows tied to annual maximum sea levels helps to partially reduce the influence of seasonal water-level variability, but does not eliminate the effects of short-term hydrometeorological factors, including wind-driven set-up and set-down, which are particularly characteristic of the shallow areas of the northeastern Caspian Sea. As a result, the observed shoreline position reflects not only long-term changes in sea level but may also be affected by the short-term hydrometeorological conditions at the time of image acquisition. Therefore, the derived shoreline positions should be interpreted with consideration of the possible contribution of short-term wind-driven set-up and set-down processes, which cannot be completely eliminated.

Sixth, an additional limitation of the study is its temporal coverage (2016–2025), which was determined by the need to use a homogeneous time series of high-frequency Landsat 8/9 and Sentinel-2 observations. The temporal compositing window (±15 days relative to the date of the annual maximum sea level) was selected empirically. Initially, the analysis was performed using narrower time intervals; however, these did not provide complete cloud-free coverage of the entire Caspian Sea. Analysis of earlier periods showed a substantial reduction in the availability of cloud-free scenes. Achieving comparable spatial coverage would have required extending the compositing window to nearly 2 months, making semi-automated shoreline extraction near annual hydrological maxima methodologically unstable. Expanding the temporal window to compensate for the lack of available imagery reduces the comparability of hydrological conditions and, consequently, decreases the reliability of interannual comparisons.

Despite these limitations, the proposed methodology provides an efficient framework for the semi-automated extraction of the Caspian Sea shoreline and water polygons, as well as for the generation of consistent multi-year time series through the combined use of optical satellite imagery and satellite altimetry. The proposed approach is not limited to retrospective analysis. The use of freely available satellite data together with the cloud-based Google Earth Engine platform enables regular shoreline updates without the need for additional field surveys. This makes the proposed methodology a promising tool for the operational monitoring of contemporary coastal changes, the assessment of coastal ecosystem degradation, and support for environmental management and decision-making.

As a result of this study, an end-to-end methodological framework for the semi-automated monitoring of coastal systems was developed, validated, and scaled using the Caspian Sea as a case study. The integration of cloud computing within the Google Earth Engine (GEE) environment, the NASA Harmonized Landsat–Sentinel (HLS) archive, high-resolution PlanetScope imagery, and satellite altimetry data enabled the automated generation of multi-year shoreline and water surface datasets, followed by manual quality control, ensuring their spatial and temporal consistency while reducing the influence of seasonal Caspian Sea level variability on the analysis.

Based on the conducted analysis, the following key conclusions can be formulated: Multicriteria validation of 12 processing scenarios using transect analysis in the DSAS v6.1.177 module, the integral ranking score, and expert visual assessment demonstrated the superiority of the modified water index MNDWI combined with 25th-percentile statistical aggregation (MNDWI_p25). The methods were validated on reference sites in the northeastern part of the Caspian Sea using PlanetScope data for 2022–2024. This algorithm provided the most accurate delineation of the water–land boundary under the conditions of the shallow shelf and developed vegetation cover of the most dynamic shallow northeastern part of the Caspian Sea. At the same time, AWEI2_median produced more spatially homogeneous water masks without erroneous “island” artifacts and can be used as an additional quality-control tool for automated classification.

The use of the DAHITI satellite altimetry archive to determine the dates of the annual maximum sea level (predominantly June–July) made it possible to generate hydrologically consistent temporal composites of satellite imagery. This approach ensured reliable interannual comparability of the results and minimized the influence of seasonal water-level fluctuations when analyzing long-term shoreline dynamics.

Based on the generated homogeneous geospatial dataset for the 2016–2025 period, a decline in the maximum Caspian Sea level of 1.33 m (from −27.13 m to −28.46 m) was recorded, resulting in a substantial transformation of its coastal zone. Over the 10-year period, the total water surface area decreased by 11,216.31 km 2 (from 390,644.77 km 2 to 379,428.46 km 2), while the total shoreline length decreased from 19,907.94 km to 15,551.79 km.

It was established that the most pronounced transformation occurred in the Kazakh sector of the northeastern Caspian Sea, where, due to the extremely gentle bottom slopes, even relatively small decreases in sea level result in substantial horizontal shoreline displacements. Shoreline retreat exceeded 20 km in some areas. Regression of the sea caused the drying of lowland territories (Tyuleniy Islands Archipelago), the connection of large island systems (Chechen Island) to the mainland, and the disappearance of extensive shallow bays (along the northeastern coast). The peak of water surface degradation occurred during the 2022–2024 period, when 4,359.83 km 2 of the water surface was lost. The obtained results demonstrate the high sensitivity of the northeastern Caspian Sea to the ongoing sea-level decline and emphasize the need for its regular satellite monitoring.

The developed methodology for shoreline extraction, combining automated processing with manual quality control, and the resulting multi-year geospatial dataset provide a foundation for regular satellite monitoring of the Caspian Sea coastal zone. The obtained results can be used for studying contemporary coastal transformation processes, assessing the condition of coastal ecosystems, and supporting the management of natural resources in the Caspian region.