Managing Flood Hazard in a Complex Cross-Border Region Using Sentinel-1 SAR and Sentinel-2 Optical Data: A Case Study from Prut River Basin (NE Romania)

Abstract

In this study, an alternative solution for flood risk management in complex cross-border regions is presented. In these cases, due to different flood risk management legislative approaches, there is a lack of joint cooperation between the involved countries. As a main consequence, LiDAR-derived digital elevation models and accurate flood hazard maps obtained by means of hydrological and hydraulic modeling are missing or are incomplete. This is also the case for the Prut River, which acts as a natural boundary between European Union (EU) member Romania and non-EU countries Ukraine and Republic of Moldova. Here, flood hazard maps were developed under the European Floods Directive (2007/60/EC) only for the Romanian territory and only for the 1% exceeding probability (respectively floods that can occur once every 100 years). For this reason, in order to improve the flood hazard management in the area and consider all cross-border territories, a fully remote sensing approach was considered. Using open-source SAR Sentinel-1 and Sentinel-2 data characterized by an improved temporal resolution, we managed to capture the maximum spatial extent of a flood event that took place in the aforementioned river sector (middle Prut River course) during the 24 and 27 June 2020. Moreover, by means of flood frequency analysis, the development of a transboundary flood hazard map with an assigned probability, specific to the maximum flow rate recorded during the event, was realized.

1. Introduction

When it comes to floods, Romania is among the most affected countries [1,2,3,4,5]. Starting at the beginning of 21st century, Romania witnessed flood events at increasingly high levels and frequencies, with floods occurring every two years (2005, 2006, 2008, 2010, 2013 and 2014) [4,6,7,8,9]. The negative consequences associated with the 40 major flood events recorded in Romania in the last 20 years (235,105 people affected, 10,731 people left homeless, 243 deaths, €2.6 billion worth of damage) ranks Romania as the second place in Europe after Russian Federation in terms of flood events incidence and people killed [10]. Many of these catastrophic flood events experienced by Romania in the last 20 years took place in the north-eastern part of the country, an area crossed by the catchments of Siret and Prut Rivers, some of the largest river basins in Romania [11]. In this area, historical flow rates for the entire Romanian territory were recorded: 4650 m³/s in the Siret River in 2005—the maximum value ever recorded, and 4240 m³/s on the Prut River in 2008—the second maximum value recorded for Romania [4]. The initial communicated flow rate for Prut River in 2008 was 7140 m³/s, which would have been the maximum ever recorded for Romania, but the Stânca-Costești Reservoir located on the Prut River decreased the discharge value at Stânca-Costești hydrometric station to 1260 m³/s [12,13].

With the accession of Romania to the European Union (EU) in 2007, the country was forced to align its flood risk management activities in compliance with the European Floods Directive (EFD) 2007/60/EC specifications [3,14]. Consequently, flood hazard maps (FHM), flood risk maps (FRM) and flood risk management plans (FRMP) were developed for all areas with potentially significant flood risk, identified within the country’s main hydrographic basins and by the district administrations [6]. Although the process of flood hazard assessment and production of flood hazard maps in Romania was carried out in concordance with current scientific standards and state of the art technology, comprising mainly of advanced hydrological and hydraulic modeling along with LiDAR digital elevation models (DEMs) and collaboration with other EU members (e.g., Danube Flood Risk Project, Monitoring of extreme flood events in Romania and Hungary using Earth Observation data Project) [15], a series of issues regarding the elaboration of these products along the transboundary Prut River (natural border between EU country Romania and non-EU countries Ukraine and Republic of Moldova) can be addressed [9,16]. Here, even though Prut River was responsible for catastrophic flood events that affected all three countries [4,12,17,18], the flood hazard was quantified by using simplified Geographic Information System (GIS) procedures that involved the reconstruction of the water levels produced at certain events, only for the Romanian territory and only for the 1% probability. Moreover, according to the EU Member State Report: RO—Romania on Assessment of Flood Hazard and Flood Risk Maps prepared in December 2014, no transnational cooperation project regarding the development of flood hazard and risk maps (FHRM) was initiated or implemented.

This aspect was favored by a cumulus of obstacles which consisted of the overlay of the Prut River with the European Union eastern border and the existence of no harmonized legislation concerning flood risk management [19]. Due to these issues, even if the EFD regulates the collaboration between EU member states and third-party countries in terms of effective flood prevention and mitigation, in practice, a cross-boundary cooperation it is difficult to implement [20]. In this particular case, one of the main disadvantages is represented by the fact that terrestrial laser scans and acquisition of LiDAR digital elevation models (essential data for hydraulic modeling techniques and elaboration of probabilistic flood hazard maps) cannot be performed at large scales. The existence of LiDAR data only for the Romanian territory is not able to offer an ample perspective regarding the flood hazard in the area along Prut River, the flood hazard maps developed here may underestimate the flood impact. However, despite all the mentioned impediments, a joint operational program between Romania, Ukraine and Republic of Moldova was implemented and completed in 2017. Following this project, called EASTAVERT PROJECT “The prevention and protection against floods in the upper Siret and Prut River Basins, through the implementation of a modern monitoring system with automatic stations”, significant progress has been made in terms of flood hazard and risk mapping within the transboundary area of Prut River.

Usually, the development of flood hazard maps involves hydrodynamic models, based on a series of data that describe statistical river discharge or rainfall values and high precision DEMs (e.g., LiDAR), along with remotely sensed maps of flood extents for calibration and validation purposes [8,21,22,23,24]. In the absence of long-term hydrometric monitoring data and high precision DEMs, flood hazard assessment and flood hazard maps can be realized on the basis of remote sensing data [25,26,27]. Flood probability maps can be obtained by correlating satellite images acquired regularly throughout the stages of a flood event, and flow rates values recorded at hydrometric stations located preferably upstream of the area affected by the flood [28]. In this context, one emergent solution to manage the flood hazard in complex cross-border regions could rely on the use of synthetic aperture radar (SAR) imagery [29,30]. Due to its active radar which operates day and night in all weather capabilities, this type of data is invaluable for uninterrupted monitoring and mapping of flood events [31,32,33,34]. Flood extent maps derived from SAR data can support the elaboration of flood hazard maps and therefore specific flood risk management measures and activities [35]. Until recently, the use of SAR data in similar applications was rather limited, due to the high costs and the complexity of their interpretation by non-expert users [36]. The launch of European Space Agency (ESA) Copernicus Program and its particular SAR mission Sentinel-1 in April 2014 brought up new opportunities in terms of near real-time flood disaster monitoring, flood extent mapping and emergency management [37,38,39]. Given the new improved capabilities in terms of global coverage, free availability, 3–6 days revisit time and 10 m spatial resolution, Sentinel-1 SAR imagery proved its efficiency in numerous disaster-related studies that focused on similar topics [40,41,42,43,44]. Moreover, the practice of using Sentinel-1 SAR data in combination with cloud-free ESA Sentinel-2 optical imagery (when available) can lead to more acurate flood hazard and flood risk analysis products, such as spatio-temporal inundation maps [45].

This paper aims to delineate the maximum inundation extent of a flood event that took place over a four day period (from 24 to 27 June 2020) in a complex cross-border region of the Prut River, which acts as a natural boundary between EU country Romania and non-EU countries Ukraine and Republic of Moldova. The main objectives consist of: (1) the exploration of Sentinel-1 SAR and Sentinel-2 optical data capabilities to capture all flood development stages, considering the image availability and repetitiveness (temporal resolution); (2) the development of quick and efficient methodologies to delineate flooded areas within Sentinel-1 SAR and Sentinel-2 optical images; (3) the development of a flood hazard map with an assigned probability for the entire study area which can offer an integrated perspective without any cross-border related limitations.

2. Study Area

The Prut River basin (total surface: 27,500 km²; total hydrological network length: 11,000 km) is located in the northeastern sector of Danube basin and overlaps with the territory of three countries: Ukraine (28%–7700 km²); Romania (39%–10,725 km²); and Republic of Moldova (33%–9075 km²) [4,18,46] (Figure 1a). It springs from the Ukrainian Wooded Carpathians at an altitude of 1580 m a.s.l. and flows into the Danube River near the Romanian village of Giurgiulești at an altitude of 2 m a.s.l. [12,18]. As the last major tributary of Danube River, it represents the natural border between Ukraine and Romania (31 km border length) and between Romania and Republic of Moldova (711 km border length) [3,4,9,12].

From a morphometric point of view, the Prut River has three sectors:

• Upper watercourse sector: from springs to Cernăuţi where the river flows out of the mountain region. This river sector presents itself as a typical mountain river, with a narrow and deep valley, medium slopes (4.5°) and steep banks. The riverbed has a width of 50–70 m, depth 0.5–1.5 m and rate flow of 1.0–1.5 m/s [13];
• Middle watercourse sector: from Cernăuţi to Ungheni (plain sector) with a floodplain width of 5–6 km with low banks and small slopes (0.1–0.2°). The riverbed has a width of 50–85 m, a depth of 2–3 m and a flow rate of 0.6–1.0 m/s [13];
• Lower watercourse sector: between Ungheni and the confluence with Danube River. The lower course is characterized by small slopes (0.08–0.1°), wide floodplain (10–12 km) and low flow rate (0.5–0.8 m/s). The river has a width of 60–100 m and a depth of 2–4 m [13].

Due to the lack of a legislation to ensure a transparent cooperation between these three countries, we obtained hydrological data just for the Romanian sector of Prut River (the common border with Ukraine and Republic of Moldova) [3] (Figure 1b). Therefore, within the Romanian territory, Prut River has nine gauging stations from the entry of Prut River to Prut-Danube confluence as follows: Oroftiana (included only water level measurements), Rădăuți-Prut, Stânca Aval, Ungheni, Prisăcani, Drânceni, Fălciu, Oancea and Șivița [3,4,12] (

Table 1. The main characteristics of the Prut River gauging stations located on the Romanian territory (see Figure 1a).

Gauge Station	Year of Inauguration	Elevation (m a.s.l.)	Latitude	Longitude	Max. Water Level Recorded (cm)	Max. Flow Rate Recorded (m3/s)	Date of Max. Flow Rate
Oroftiana 1	1976	123.47	48°11′12″	26°21′04″	876	-	-
Rădăuți-Prut	1976	101.87	48°14′55″	26°48′14″	1130	4240	28 July 2008
Stânca Aval	1978	62.00	47°47′00″	27°16′00″	512	1050	31 July 2008
Ungheni	1914	31.41	47°11′04″	27°48′28″	654	796	8–10 July 2010
Prisăcani	1976	28.08	47°05′19″	27°53′38″	622	900	9–10 July 2010
Drânceni	1915	18.65	46°48′45″	28°08′04″	718	736	17–18 July 2010
Fălciu	1927	10.04	46°18′52″	28°09′13″	650	722	19 July 2010
Oancea	1928	6.30	45°53′37″	28°03′04″	622	757	24 April 1979
Șivița 2	1978	1.66	45°37′10″	28°05′23″	-	-	-

¹ Oroftiana gauge station was designed for measuring only the water level. ² The Prut flow rate and water level at Șivița gauge station were directly influenced by Danube waters.

). The first installed gauging station was at Ungheni (1914) and the newest at Șivița (1978) [12]. On Ukrainian territory, one gauging station was installed in Cernăuţi city that included water level measurements (as in the case of Oroftiana gauging station) [12].

3. Database and Methodology

The Sentinel-1 SAR and Sentinel 2 data were obtained for the upper area within Cernăuţi–Ungheni sector, from Cernăuţi city (Ukraine) to Cotu Miculinți village (Romania)—Tețcani village (Republic of Moldova) (Figure 2). Along this river sector are located 26 Ukrainian settlements (of which two cities) with 323,788 inhabitants, 11 Romanian settlements (of which one city) with 20,059 inhabitants and seven Moldavian settlements (one city) with 13,406 inhabitants. The water level and flow rate for this sector was monitored by two of nine gauging stations: Oroftiana and Rădăuți-Prut gauging station (Figure 1b, Table 1).

3.1. Sentinel-1 SAR and Sentinel-2 Data

Synthetic Aperture Radar (SAR) is an imaging active radar system mounted on a moving platform that operates independently of Sun illumination and cloud coverage [47]. Such systems make use of an antenna to generate, transmit and collect electromagnetic pulses and backscattered echoes, respectively, in order to acquire information about Earth’s features and thus create SAR images [48]. ESA Sentinel-1 mission consists of a constellation of two polar-orbiting satellites (Sentinel-1A and Sentinel-1B), equipped with a C-band (wavelengths [λ] = 5.6 cm) active SAR sensor, which facilitates the acquisition of imagery regardless of day, night or weather conditions every 6 days or even 1–3 days for areas such as Europe [49,50,51]. The data were acquired under different configurations in four different modes: Stripmap (SM), Interferometric Wide swath (IW), Extra-Wide swath (EW), Wave (WV); and distributed freely at various product levels (Level-0, Level-1, and Level-2), each with specific resolutions, polarizations (single polarization HH or VV and dual polarization HH+HV or VV+VH) and extents [39,52].

Level-1 data are the generally available products intended for most data users and applications. Level-1 products are distributed as Single Look Complex (SLC) and Ground Range Detected (GRD) [53,54]. SLC data consist of focused complex imagery that uses the full C signal bandwidth, having amplitude and phase information preserved [55]. At the same time, an SLC dataset provides the highest possible spatial resolution [56,57]. Level-1 Ground Range Detected (GRD) comprises of focused SAR data that has been detected, multi-looked and projected to ground range using an Earth ellipsoid model, such as WGS84, storing only amplitude (pixel intensity) without any remaining phase information [44]. The information recorded in both Level-1 product types lies in the measured echoes of the backscattered signal at the C-band wavelength [58]. The backscattering coefficient, or sigma nought (σ⁰), is the measure of the incident microwave radiation scattered by the radiated terrain. The amplitude and phase of the backscattered signal depends on the land surface properties of terrain elements and their electromagnetic characteristics (e.g., materials) [59,60]. One of the main advantages of using SAR images in flood mapping applications is that they show a high contrast between ground and water surfaces [61]. Smooth open water or flooded areas act as ideal specular reflectors and are characterized by low SAR backscatter values [43,51]. Therefore, in SAR imagery, the water covered surfaces appear black [62].

Sentinel-2 consists of a constelation of two satellites (Sentinel-2A launched in June 2015 and Sentinel-2B launched in March 2017), equipped with a multi-spectral instrument (MSI). The multispectral sensor aquires optical imagery in 13 spectral bands from VIS/NIR to SWIR at a spatial resolution ranging from 10 to 60 m and at a high revisit time (5 days at the equator with two satellites and every 2–3 days at mid-latitudes.) The data are freely available under two distinct product levels (Level-1C and Level-2A) [63,64]. Level-1C data are distributed as Top of the Atmosphere (ToA) reflectance in cartographic geometry (100 × 100 km² orthoimages found in UTM/WGS84 projection). Level-2A data are based on Level-1C data and are atmospherically corrected to Bottom of Atmosphere (BoA) reflectance [65,66].

For the purpose of this study, freely available GRD dual-polarized (VV/VH) Sentinel-1 SAR images acquired by the sensor under interferometric wide-swath (IW) mode (the default acquisition mode over land) at a pixel spacing of 10 m and a swath of 250 km were downloaded using the Copernicus Open Access Hub [42,67,68,69,70]. Level-2A Sentinel-2 atmospherically corrected products were downloaded using the QGIS application and the semi-automatic classification plugin [71].

Due to their specific characteristics (especially the improved temporal resolution), the Sentinel-1 and Sentinel-2 imagery obtained managed to depict all development stages of a flood event that took place on a shared Prut River sector located between Ukraine, Romania and Republic of Moldova during 24 to 27 June 2020 (Figure 2). This aspect allowed near real-time flood monitoring as the images were available for download in less than 24 h since acquisition by sensor. A characterization of the images used is listed in

Table 2. Characteristics of the Sentinel-1 SAR and Sentinel-2 optical images used in the study (see Figure 2).

ID	Image Identifier	Satellite	Date Acquired/Time	Product Type	Mode	Polarization
A	S1A_IW_GRDH_1SDV_20200624T042904_20200624T042929_033154_03D73C_540B	S1A	2020-06-24/T04:29:04.880Z	GRD	IW	VV VH
B	S2B_MSIL2A_20200624T090559_N0214_R050_T35UMP_20200624T123246	S2B	2020-06-24/T09:05:59.024Z	LEVEL-2A	-	-
C	S1B_IW_GRDH_1SDV_20200625T042007_20200625T042036_022185_02A1B4_C73E	S1B	2020-06-25/T04:20:07.597Z	GRD	IW	VV VH
D	S1B_IW_GRDH_1SDV_20200626T160915_20200626T160940_022207_02A256_2724	S1B	2020-06-26/T16:09:15.797Z	GRD	IW	VV VH
E	S2B_MSIL2A_20200627T092029_N0214_R093_T35UMP_20200627T121756	S2B	2020-06-27/T09:20:29.024Z	LEVEL-2A	-	-
F	S1A_IW_GRDH_1SDV_20200627T160148_20200627T160213_033205_03D8C4_F195	S1A	2020-06-27/T16:01:48.625Z	GRD	IW	VV VH

.

3.2. SAR and Optical Data Pre-Processing

The GRD Sentinel-1 images downloaded were pre-processed using Sentinel Application Platform (SNAP) software, an ESA open-source architecture [72,73]. The pre-processing steps were based on a workflow proposed by Filipponi [74] and consisted mainly of (Figure 3a): data preparation, area and band subset, image calibration, speckle filtering, terrain correction and data conversion to dB.

In the first step, in order to reduce the processing time and memory use, the data were sub-sampled according to study area limit and polarization. Next, the precise orbit file, thermal and border noise were respectively applied and removed [37,75]. Even if both VH and VV polarized Sentinel-1 data can be employed in flood monitoring and mapping research [59,76], for the purpose of this study the VV polarization was preferred over VH polarization, as it is more suitable for delimiting flood waters within mixed, forest and agricultural land areas [35,50,77].

Image calibration consisted of the conversion of raw DN values to radiometrically calibrated SAR backscatter or sigma nought values (σ⁰) [78]. Speckle-filtering implies a procedure to increase image quality by reducing speckle noise. In this context, many speckle filters based on smoothing windows on a weighted summation of neighboring pixels were developed [77]. In this study, a Lee-Sigma filter with a window kernel size 5 × 5 pixels was applied [79]. Terrain corrections were performed in order to reduce the distortions determined by the varying viewing angles involved in the image acquisition process [80]. The SAR images underwent range Doppler terrain correction using the Shuttle Radar Topography Mission (SRTM) (1 arc second) DEM as a reference. During this operation, the images were corrected for topographical distortions determined by the side looking SAR mechanism, geolocated and re-projected to Universal Transverse Mercator (UTM) zone 35 N, with a 10 m ground sampling [81].

In order to ensure that all four datasets considered in this study were spatially aligned, an image co-registration was performed. This was carried out using the co-registration tool found in SNAP software and the Sentinel-1 image acquired on 24 June 2020 as a reference. In the last step, the calibrated backscatter coefficient values were converted to dB logarithmic units and exported for further water and non-water pixel delineation operations [76,82].

Level-2A Sentinel-2 imagery was downloaded using the the QGIS application and the semi-automatic classification plugin as Bottom of Atmosphere (BoA) reflectance and its pre-processing consisted of: data merging using only the 10 m spatial resolution spectral bands (band 2—blue, band 3—green, band 4—red and band 8—NIR), study area data clipping and data re-projection to Universal Transverse Mercator (UTM) zone 35 N [63]. To spatially match the Sentinel-2 optical images with the Sentinel-1 SAR images, a co-registration step was also performed using the SNAP application [83] (Figure 3b).

3.3. Water Pixel Extraction from SAR and Optical Data

Various methods have been developed for the detection of water surfaces and thus for flood mapping using radar SAR images [29,84]. Among them, considering water’s well-defined backscattering signature, backscatter histogram thresholding is the most widely used [85,86]. The threshold method can be used to differentiate between the water pixels, represented by low values from the no-water pixels characterized by high values within the image histogram [87].

Thresholding implies SAR backscatter intensity histogram analysis and the estimation of water and non-water pixel distribution. This method is based on histogram manipulation and on the distinction of two partially overlapping distributions that correspond to water and non-water pixels [88,89]. The intersection point between wet and dry distributions determines the threshold value needed for the two classes separation (i.e., open water surfaces and non-water surfaces) [62,90]. For Sentinel-1 SAR images considered in the present study, the extraction of the two distinct classes was conducted based on a threshold value which was initially obtained and then applied to all image’s histograms (Figure 3).

The threshold value was determined by analyzing the backscatter coefficient values, corresponding to a series of points manually digitized by means of visual interpretation, along the Prut River course and within the most visible flooded areas. Subsequently, for each Sentinel-1 image, a different set of points was acquired in areas that overlap with low intensity dark pixels that represent water. To make sure that the point locations coincide exclusively with water surfaces, a cross-examination with the cloud-free areas of two Sentinel-2 multispectral images dating from 24 to 27 June 2020 was performed. In this way, the extraction of the matching raster pixel values set an interval for the true water areas found within each image. The interval highest values that were equal to −15 dB (Image A), −16 dB (Image B), −17 dB (Image C) and −15 dB (Image D), established the threshold between water and land classes. Consequently, all values below these thresholds were considered representative for the water class and respectively for the flooded areas. Next, a binary geotiff layer showing the water and dry areas was obtained for all four images.

For data integration purposes, a shapefile vector format conversion was also realized. The detected water class consists of flooded areas and permanent water bodies. To highlight only the flooded areas, it was necessary to subtract the permanent water bodies from the initial water class dataset [91]. In this case, the permanent water bodies acquired from OpenStreetMap were used. The final step consisted of the export of four distinct datasets regarding the flooded areas, one dataset for each corresponding day of the flood (Figure 3).

The water pixel extraction from Sentinel-2 optical data was realized by means of Normalized Difference Water Index (NDWI) spectral index calculation and simple thresholding. The McFeeters NDWI [92] maximizes the reflectance of the water bodies in the green band and minimizes it in the NIR band [93,94]. NDWI was calculated according to the Equation (1):

$$\mathrm{NDWI} =\frac{\mathsf{\rho }\mathrm{Green}- \mathsf{\rho }\mathrm{NIR}}{\mathsf{\rho }\mathrm{Green}+ \mathsf{\rho }\mathrm{NIR}}$$

(1)

Considering this multispectral advantage, the NDWI is a simple and effective method to extract water bodies [95]. The NDWI is a dimensionless quantity which varies between −1.0 and 1.0. The positive NDWI values indicate the presence of water [96,97]. In this study, the water pixels were limited to the NDWI positive values. The flooded areas for the two Sentinel-2 datasets were delimited by permanent water bodies’ subtraction. The final result consisted of two layers, which illustrate the flooded areas within the cloud-free regions of the Sentinel-2 satellite images.

3.4. Water-Affected Pixel Extraction from Optical Data

The identification of the flooded areas can be realized even after the water has already withdrawn, by analyzing the flood’s effects on the soil [98]. Using the Sentinel-2 image acquired on 27 June 2020 (on the last day of the flood), the extraction of the flood affected pixels was possible. By means of Normalized Difference Vegetation Index (NDVI) and simple thresholding, we managed to explore the high contrast between the bare soil flood affected areas and the surrounding high vegetated agricultural land. NDVI was calculated according to the Equation (2):

$$\mathrm{NDVI} =\frac{\mathsf{\rho }\mathrm{NIR}- \mathsf{\rho }\mathrm{Red}}{\mathsf{\rho }\mathrm{NIR}+ \mathsf{\rho }\mathrm{Red}}$$

(2)

NDVI values range between −1.0 and +1.0. The NDVI values near zero and negative values indicate nonvegetated areas, such as rock, soil or water [99]. In this study, the water-affected areas were restricted through visual interpretation of the NDVI product and Sentinel-2 natural color composite image to values lower than 0.35. The threshold value was determined based on a set of points manually digitized in the visible flood-affected areas. As the result also included the areas covered with water, a Sentinel-1/Sentinel-2 combined flooded area from 27 June 2020 and permanent water bodies subtraction was realized. The final result consisted of a layer which depicted all the flood-affected areas from 24 to 27 June 2020.

3.5. Water Pixel Validation

The subtracted flooded areas and water-affected areas underwent a validation process in order to confirm the proposed methodologies and the determined threshold values. The validation for the Sentinel-1 extracted flood areas was realized by vector overlay analysis (Figure 3). Two vector datasets delineating the flooded areas on 26 June 2020 within the study area, one obtained by the proposed methodology and the other one offered by Copernicus Mapping Emergency Management Service (CMEMS-M), EMSR445 activation code product [100], were compared. The analysis revealed an overlay degree of over 95%, which was considered satisfactory for validation purposes of the Sentinel-1 water pixel extraction methodology adopted in the present study. Moreover, an aerial image acquired in the proximity of Rădăuți-Prut cross-border point between Romania and Republic of Moldova, which depicted the same flood event was also used for data visual interpretation and confirmation (Figure 4). The validation for Sentinel-2 flooded areas and water-affected areas was realized by means of visual interpretation.

3.6. Flood Frequency Analysis

In order to estimate the return periods of the flood events, their corresponding probabilities and probable discharges at Rădăuți-Prut gauging station, we considered the Gumbel’s distribution method, first introduced in 1941 [101]. Gumbel probability distribution [102] is expressed as Equation (3):

X_T = X + Kσ_x,

(3)

where, X_T is the probable discharge with a return period of T years (the magnitude of floods); X is the mean flood; σ_x is the standard deviation of the maximum annual flow rates; K is the frequency factor, expressed as Equation (4):

$$\mathrm{K}=\frac{{\mathrm{Y}}_{\mathrm{T}}-{\mathrm{Y}}_{\mathrm{n}}}{{\mathrm{S}}_{\mathrm{n}}},$$

(4)

where, Y_n is the expected mean and S_n is the expected standard deviations of reduced extremes found within Gumbel’s table, depending on the sample size, and Y_T is the reduced variate, expressed as Equation (5):

$${\mathrm{Y}}_{\mathrm{T}}=-\left[\mathrm{LnLn}\left(\frac{\mathrm{T}}{\mathrm{T}-1}\right)\right],$$

(5)

This is one of the most common distribution functions used to predict flood return periods and their probabilities based on time series discharge data [103,104].

4. Results and Discussion

4.1. Flood Hazard Map

All seven datasets corresponding to the flooded areas and flooded-affected areas were merged into one. In this way, a new flood extent map derived from satellite images was obtained (Figure 5).

As the map developed managed to capture the maximum flood extent, the conversion to a flood hazard map with an assigned probability was taken into consideration. To achieve this, a return period was established for the maximum flow rate of 2965 m³/s recorded at Rădăuți-Prut gauging station on 26 June 2020 at 10:00 AM (considering the whole flood event). This value was the only one considered representative for the flood wave, as at Oroftiana gauging station this type of data is not available. In order to carry out this statistical operation within the study area, a sample of the maximum annual flow rates recorded at Rădăuți-Prut gauging station in the last 42 years (from 1978 to 2019) was analyzed.

The flood frequency analysis consisted of: (1) the computation of the mean flood X and standard deviation σ_x considering the maximum annual flow rates (in our case these values were equal to: X = 996.21 and σ_x = 755.52); (2) selection of the corresponding values of Y_n and S_n from Gumbel’s table (in our case, given the 42 year sample size, the values were 0.5448 and 1.1458, respectively); (3) depending on the given return periods T, the reduced variate Y_T was calculated using Equation (5); (4) the calculation of the flood frequency factor K using Equation (4); (5) the probable flood discharges corresponding to different return periods and probabilities at Rădăuți-Prut gauging station were computed according to Equation (3). Detailed information about the computational process can be found in

Table 3. Flood frequencies estimated according to Gumbel’s distribution method at Rădăuți-Prut gauge station on the Prut River (see Figure 1b).

Return Period (T) in Years	Probability of Occurrence (%)	Reduced Variate (Yt)	Frequency Factor (K)	Computed Flood Discharges (XT) (m³/s)
1000	0.1	6.907255	5.5528495	5191.5417
100	1	4.600149	3.539316	3670.2634
33.3	3	3.491366	2.571624	2939.1456
20	5	2.970195	2.116770	2595.4911
10	10	2.250367	1.488538	2120.8451
5	20	1.499939	0.833600	1626.0222

.

Considering the computed flood discharges (Table 3) and the maximum flow rate of 2965 m³/s registered at Rădăuți-Prut gauging, a probability of occurrence and a return period was determined for the flood event that took place between 24 and 27 June 2020, on the middle course of Prut River. Consequently, this led to the conversion of the flood extent map derived initially from Sentinel-1 data into a flood hazard map with an assigned probability of 2.97% and a return period of 33.67 years (Figure 5).

4.2. Flood Wave Development Stages Based of Sentinel-1 SAR Satellite Images

4.2.1. Flood Status on 24 June 2020

On the first day of the flood, the most affected areas were located within the Ukrainian territory and corresponding Prut River sector (Figure 6). Here, during this day, 1609 ha and 1091 ha, respectively, were covered by water. At the same time, the flood wave started spreading across Romania, flooding 951 ha and the Republic of Moldova, where 303 ha were flooded. Meanwhile, the figures within the Romanian–Ukrainian and Romanian–Moldavian shared river sectors stood at 1003 ha and 772 ha. At 4:00 AM (the approximate time of image acquisition), the flow rate registered at Rădăuți-Prut gauging station was 990 m³/s (probability of occurrence 32.84%). The maximum value of the day, 1585 m³/s was recorded at 23:00 PM (probability of occurrence 20.51%) (

Table 4. The maximum flow rates recorded during the flood event at the Rădăuți-Prut gauge station and the corresponding probability of occurrence estimated using the Gumbel’s method.

Date	Flow Rate (m3/s)/Time of Image Acquisition	Probability of Occurrence (%)/Time of Image Acquisition	Maximum Flow Rate (m3/s) of the Day	Maximum Probability of Occurrence (%) of the Day
24 June 2020	990/4:29 AM	32.84/4:29 AM	1585	20.51
25 June 2020	1705/4:20 AM	19.07/4:20 AM	2610	4.97
26 June 2020	2860/4:09 PM	3.08/4:09 PM	2965	2.97
27 June 2020	2410/4:01 PM	5.38/4:01 PM	-	-

).

4.2.2. Flood Status on 25 June 2020

On 25 June, the situation escalated and the flood reached its peak. The most affected country was Ukraine with 2275 ha covered by water, followed by Romania with 1234 ha and Republic of Moldova with 532 ha. Regarding the three river sectors considered, the Romanian–Ukrainian was the most affected (1686 ha). On this day, the total area occupied by water reached 4042 ha (Figure 6). The flow rate recorded at 4:00 AM (the approximate time of image acquisition) at Rădăuți-Prut was 1705 m³/s (probability of occurrence 19.07%). The maximum value of the day, 2610 m³/s, was recorded at 23:00 PM (probability of occurrence 4.97%) (Table 4).

4.2.3. Flood Status on 26 June 2020

Starting on 26 June, the flood wave diminished over Ukrainian territory and shared Romanian–Ukrainian river sector. The flood that considerably affected Ukraine on 24 and 25 June decreased by half in terms of affected areas and remained constant over Romania. During this day, the flood waters began to cover territories located within the Republic of Moldova and Romania, the shared river sector being the most affected (Figure 6). The flow rate recorded at 4:00 AM (the approximate time of image acquisition) at Rădăuți-Prut was 2860 m³/s (probability of occurrence 3.08%). The maximum value of the day and also over the entire duration of the flood, 2965 m³/s, was recorded at 10:00 PM (probability of occurrence 2.97%) (Table 4).

4.2.4. Flood Status on 27 June 2020

On the last day of the flood event, few territories were still covered by water. As the flood wave followed its course, the downstream Romanian–Moldavian sector was the most affected (623 ha) (Figure 6). At 4:00 AM, the flow rate registered at Rădăuți-Prut had a value of 2410 m³/s (probability of occurrence 5.38%). Starting on 28 June 2020, the flow rates dropped below 1000 m³/s marking the end of the flood event (Table 4).

4.3. The Sentinel-1 SAR and Sentinel-2 Optical Data: Applicability and Limitations for Flood Assessment

In this study, an alternative solution for flood hazard map development in complex cross-border regions is presented. In this cases, due to different flood risk management legislative approaches, there is a lack of joint cooperation between the involved countries. As a main consequence, LiDAR high resolution digital elevation models and accurate flood hazard maps obtained by means of hydrological and hydraulic modeling are missing or are incomplete. This is also the case of Prut River middle course, which acts as a natural boundary between Romania and non-EU countries Ukraine and Republic of Moldova. Here, flood hazard maps were developed under EFD by the National Administration “Romanian Waters”, only for the 1% exceeding probability (floods that can occur once at every 100 years), using simplified Geographic Information System (GIS) procedures. For this reason, in order to improve the flood hazard management in the area, alternative solutions were considered. Consequently, such an alternative could rely on the use of Sentinel-1 SAR and Sentinel-2 optical imagery.

The launch of the Sentinel-1 mission in 2014 has significantly improved flood mapping capabilities mainly due to: the mission active radar sensor (SAR) that captures images both day and night (regardless of weather and its characteristics—clouds, fog, aerosols), high spatial resolution (10 m) and also due to the ability to acquire images at shorter intervals of time, even on consecutive days. Therefore, the mapping of floods that occurred in the past using Sentinel-1 SAR data can represent a viable solution in the attempt to identify areas with high flood potential and in the development of flood hazard maps. Sentinel-2 optical imagery has limited application in flood mapping, due to the presence of clouds and other atmospheric impediments during flood events. However, this type of data can complement the Sentinel-1 SAR data and improve the flood monitoring process. This can be realized through flood area delineation within the image’s cloud-free regions and flood-affected areas delineation within post-event images.

Even if the merging of the daily flood extents derived from Sentinel-1 SAR and Sentinel-2 optical imagery led to a map that illustrates a maximum flood extent, which subsequently was considered for conversion to a flood probability map, some limitations regarding the water pixel extraction from Sentinel-1 and Sentinel-2 data needs to be addressed. Among them which is the most challenging is the possibility that in the time interval between the acquisition of the satellite images, some areas may have been flooded and thus missed. This may lead to the underestimation of the flooded areas. To address this challenge and achieve a complete flood characterization, the flood-affected areas were extracted using the 27 June 2020 Sentinel-2 satellite image (acquired in the last day of the flood).

Realizing a flood hazard map with an assigned probability using exclusively satellite imagery is difficult to achieve. This is limited by correlations with statistical data from gauging stations (used to determine the probability of flooding) and by the availability and continuity of satellite images throughout the entire temporal development of the flood event. Regardless of these issues, in an attempt to map a flood event that occurred on the Prut River middle course during 24–27 June 2020, both conditions mentioned above were accomplished and Sentinel-1 SAR and Sentinel-2 optical satellite images managed to capture the daily development stages of the flood event and the flood-affected areas.

In this way, an overview regarding the start date of the flood, the withdrawal rate, the affected areas throughout the flood wave and the maximum extent was obtained. Depending on the surface of the flooded areas within the Sentinel-1 images, we managed to determine that the flood started on Romanian territory on 24 June 2020, reached its peak on 25 June 2020 and experienced a decrease in intensity on 26 June 2020.

Even if the merging of the daily flood extent and flood-affected areas derived from Sentinel-1 SAR and Sentinel-2 optical imagery led to a map that illustrates the maximum flood area, this cannot be considered a flood hazard map as it does not contain any information regarding any exceedance probability. Therefore, following the proposed methodology, a statistical flood frequency analysis based on the maximum flow rate values located at a nearby gauging station—Rădăuți-Prut—was conducted and a flood return period was assigned.

5. Conclusions

One of the most effective methods of reducing the catastrophic impact of floods is the prior scientific assessment of flood hazard, its prediction, modeling and mapping. Flood hazard maps can significantly contribute to reducing the negative effects associated with floods. In recent years, advances in remote sensing and Geographic Information Systems (GIS) have revolutionized flood hazard assessment methodologies, thus facilitating new interdisciplinary approaches. Under certain conditions, remote sensing and time series discharge data can offer an alternative solution for generating flood hazard maps, especially in complex cross-border regions such as in our study area.

Using Sentinel-1 SAR and Sentinel-2 optical data along with GIS, we managed to: (i) capture the maximum extent of a flood event that took place on the middle course of Prut River, a complex cross-border region shared by the EU member Romania and non-EU countries Ukraine and Republic of Moldova, where the flood hazard maps are either incomplete or missing; (ii) develop two simple and efficient water pixel extraction methodologies considering the high contrast between water, water affected areas and land within the Sentinel-1 SAR and Sentinel-2 data; (iii) explore the improved temporal resolution of combined Sentinel-1 SAR and Sentinel-2 optical data and obtain a transboundary flood hazard map with an assigned probability.

This study proved that alternative solutions for flood hazard maps obtained by means of hydrological and hydraulic modeling exist and can be successfully implemented. Considering the legislative and emergency situations issues, determined by the cross-border regions, this kind of methodology can significantly contribute to the improvement of flood risk management and can stimulate international cooperation. Moreover, given the open-source character of Sentinel-1 and Sentinel-2 data and the adaptability of the proposed methodologies, this workflow can be easily extrapolated to other areas where the flood hazard has not been yet regulated and no related information exists.