Unveiling Torrential Flood Dynamics: A Comprehensive Study of Spatio-Temporal Patterns in the Šumadija Region, Serbia

Abstract

This paper presents a comprehensive analysis of flood frequency and a spatio-temporal characterization of historical torrential floods in the Šumadija region using water discharge datasets and documented events. A chronology of 344 recorded torrential flood events, spanning from 1929 to 2020, illustrates the region’s vulnerability, with a death toll exceeding 43. The study defines the intra-annual primary and secondary peaks of torrential flood occurrences and explains their spatial distribution. Furthermore, the identification of suitable probability distribution functions underscores the necessity of tailored approaches for effective flood risk management in this diverse geographical environment. The study employed Flood Frequency Analysis (FFA) and goodness-of-fit tests, including the Kolmogorov–Smirnov (K-S) and Cramér–von Mises (CvM) tests, to assess the frequency and magnitude of flood events and evaluate diverse distribution functions. The main results include the identification of suitable probability distribution functions for each river within the region, emphasizing the need for tailored approaches in flood risk management. Additionally, discharge values for various return periods offer crucial insights for informed decision-making in flood risk management and infrastructure planning.

1. Introduction

Currently, two primary concerns regarding water resources stand out prominently, both intricately linked to hydrological extremes: water scarcity and flooding [1]. Among the different types of floods, torrential floods are the most frequent and costly natural hazards in Serbia [2]. Torrential floods arise from the extreme rainfall showers in slope-dominant watersheds, inducing rapid maximal discharges and a high amount of sediment in the river beds due to intense soil erosion on the slopes, often coinciding with landslide movements or rock-falls [3,4]. This natural phenomenon with sudden onset is a consequence of natural conditions and factors, and its severe societal impacts can be significantly multiplied by anthropogenic influence [5]. In this regard, the historic torrential floods with disastrous material damage and the high death toll of April, May, and September 2014 in small and medium watersheds in central Serbia are the latest testimonies that the character of torrential floods appears to be a real challenge for flood risk management in Serbia [6].

The inventorization of floods is highly important for knowledge gain and is a first step towards tackling a problem. In Serbia and a large number of countries, torrential floods are not documented systematically and centralized—only the information on recent floods is scattered across different sectors and can become lost over time. However, researchers worldwide have endeavored to build good examples of flood documentation encompassing the different periods and spatial extent. In this way, torrential flood research can benefit from well-documented datasets. Archer and Fowler [7] (2021) provided a flash flood chronology for Britain with nearly 8000 events for the period between 1700 and 2020. Haltas et al. [8] (2021) compiled and analyzed an inventory of 2101 flood events in Turkey recorded between 1930 and 2020. Vennari et al. [9] (2016) presented a database for the Campania region in southern Italy, which includes around 500 flash flood events from 1540 to 2015. Diakakis [10] (2013) presented a database of floods caused by torrential rainfall in the Athens region, which happened during 1988–2010 and took 182 human lives, at the same time indicating a positive trend in terms of their frequency. Finally, Adhikari et al. [11] (2010) compiled a digitized Global Flood Inventory—GFI for the period 1998–2008 with geo-referenced flood spots.

In Serbia, Petrović et al. [12] (2014) and Petrović [4] (2021) created and analyzed the Inventory of torrential floods in Serbia over a period of 99 years, 1915–2013, with nearly 850 recorded events and more than 133 deaths caused by floods, while the revised and updated version, covering 105 years, reveals valuable information on more than 2100 flood events with more than 193 flood deaths. In this paper, a dataset of torrential floods for the Šumadija region is derived, and an extended regional chronology for the period up to 2023 is presented. In Šumadija, the suburbs of Belgrade, the capital of the Republic of Serbia, as well as the cities Kragujevac and Čačak, were endangered by torrential floods many times in history. Torrential flood waves cause more damage in lower deforested, agricultural, or urbanized parts of a watershed than in higher lands due to a stronger terrain dissection. Šumadija is known as one of the central and economically developed areas of Serbia and has been the subject of research from the aspect of natural hazards [13]. However, there is a gap in detailed research that focuses on torrential floods and maximal discharges as hydrological extremes occurring in small- and medium-sized watersheds.

Hydrological monitoring from the second half of the 20th century allows the statistical analysis of time series of discharge and water level data. Knowledge of water quantity, water discharge, and water levels appears to be a starting point for the management of water resources. Trend changes related to discharge, especially concerning hydrological extremes—floods and droughts—can be detected by analyzing data series. Milanović Pešić [14] (2023) analyzed the annual discharge variability and revealed that there is a decrease in discharge in all studied watersheds when comparing discharge data for periods 1961–1990 and 1991–2020. Despite this, large-scale torrential flooding is expected due to more frequent extreme rainfall events. The rivers in the Šumadija region have larger fluctuations in discharge—from the spring months with snow melting and extreme rainfall episodes resulting in high discharges to the summer months characterized by low water levels and even dry river beds in drier years. This leads to a conclusion that also takes into account the characteristics of relief in the upper and smaller watersheds about their torrential water regimes.

Accurate and consistent prediction of stream flows holds paramount importance across various domains, encompassing water resource management, strategy enhancement, navigation, and maintenance operations [15]. Conducting flood risk assessments is a common practice that aims to mitigate flood-induced damages in specific locales [16]. Flood Frequency Analysis (FFA) emerges as a crucial tool in this context, facilitating the assessment and mitigation of flood risk by offering insights into the frequency and severity of flood events, pivotal for infrastructure planning and risk management strategies [17]. To ensure a precise evaluation of flood dynamics and magnitude, a considerable volume of accurately recorded peak flows is imperative for robust FFA outcomes [18]. FFA, a statistical method frequently employed to estimate flood magnitude within a designated return period, finds extensive application in studies related to water resource management [19]. The conventional approach of FFA entails extrapolating the tails of the distribution to ascertain the likelihood and severity of extreme occurrences, accomplished by fitting mathematical functions to available data [20]. Notably, the FFA serves as a linchpin in engineering practices, establishing correlations between design variables and chosen hydrological risks [21].

Therefore, the focus of this paper is to present the torrential flood phenomenon in the region of Šumadija in central Serbia based on the temporal and spatial characteristics of the recorded historical torrential flood events, as well as to determine the frequency of high waters at gauges with smaller watershed areas. Additionally, this research aims to employ, for the first time, Flood Frequency Analysis (FFA) techniques to assess the frequency and magnitude of torrential flood events in the region, thereby enhancing our understanding of flood risk dynamics. Ultimately, the goals of this study include providing valuable insights to support improved flood risk assessment and enhancing preparedness measures to mitigate the impact of future torrential flood events in the Šumadija region.

2. Materials and Methods

The dataset of registered torrential floods in the Šumadija region is a part of the Inventory of torrential floods in Serbia. The definition of the subject of the Inventory—the torrential flood (see Section 1) is applied in the process of inventorization. While the spatial scope of the documentation was defined in advance by the physical-geographical characteristics—hilly-mountainous areas south of the Sava and Danube Rivers, the temporal scope of the dataset was determined by the availability of data sources. Torrential floods in Serbia are the most frequent in spring, but all seasons are considered. The criterion related to watershed area is 500 km², taking into account small- and medium-sized watersheds, except for the Ljig River, which has a larger watershed area with a torrential water regime. Although debris flows and landslides are accompanying phenomena of torrential floods, they are not documented in this dataset (e.g., a debris flow and landslide event in a watershed of the Leva River in May 2014). The omission of phenomena similar to torrential floods should be emphasized, as many of the often-called flash flood datasets also include these events. However, information on these natural hazards may appear as accompanying in the broader description of torrential flood events. The heavy rainfall events in urban areas and resulting floods are taken as urban flash floods and are not documented in this dataset. The criterion that floods in several different watersheds or even tributaries (in the case of extreme material damage and impact on the discharge of a major river) that happened on the same date are counted as multiple events differ from some other criteria that floods in different watersheds on the same date count as one event in the dataset. The preparatory phase of torrential flood documentation focused on defining the purpose of the dataset and designing the content and attributes for each event documentation. Data collection, which took the majority of the time, was followed by the organization and analysis of the data, which enabled the distribution and publication of the obtained results. This dataset is open for data supplementation for the presented period at any time and will be updated with future torrential flood events.

The information collected was compiled using a structured documentation approach. This includes the following: name of the river with the torrential flood event and the larger river basin to which it belongs, locations and settlements that suffered material losses, date of the event (YY/MM/DD), description of the event (free text), flood deaths and the number of casualties and evacuees, descriptive information on economic losses, and the source of information. The information on the first three attributes, followed by the source of information, is a minimum and mandatory for the event documentation. The collection scheme also includes a place for hydrological and meteorological data, as well as data on warnings and alerts. A particular event has a code, e.g., 1986-02-19-VM-L-M, which is composed using the date and initials of the larger river basin to which the river with the flood event belongs, the name of the river with the torrential flood event, and the name of the affected location.

The Republic’s agencies, services, ministry centers, insurance companies, local municipalities, and belonging centers for emergencies were invited to contribute; however, many efforts were made, and only poor results were obtained in this communication. In the last decade, emergency centers have been responsible for reporting on emergency events. Still, the reports mainly focus on the assessment of material damage and vary greatly in the volume of information provided. None of them included hydrological information, and some contained meteorological data. In recent years, the Republic Hydrometeorological Service of Serbia (RHMSS) has provided public reports on extreme rainfall events. Finally, there is the challenge of storing all sources of event information with regard to the sustainability of the dataset. In particular, this refers to nternet sources, whose information can be saved as PDF files, while archival documentation of newspapers for earlier events can be stored as copies.

Therefore, the most reliable historical information was found in the scientific literature [22,23], and hydrological data for several gauged watersheds were obtained from the hydrological yearbooks (RHMSS). The maximum discharges for the recorded torrential flood events were collected, with particular attention paid to observing the discharge on the day before the flood. The specific maximum discharge, q_maxsp, which is a ratio between the peak discharge and the catchment area, was used as an indicator of the extremeness of the flood (m³s⁻¹km⁻²). If the event documentation is confirmed by one or both of these sources, then the event documentation is marked as fully reliable and validated. In addition, some of the documented events have hydrographs of torrential flood waves. Furthermore, valuable information was gathered from the archival documentation of the newspapers. Certainly, there is an uncertainty in event information from newspaper sources, which barely could be excluded, in case there are no other sources on a particular event. When a description of an event suits a torrential flood and relates to smaller watersheds and belonging and affected settlements, we used newspapers as a source of information. A good example of this case is the torrential flood of June 1969, for which an extensive description of the event and property damage is provided along with several photographs. However, when the main or only source of information is newspapers, the event documentation is marked as medium reliable. In this chronology, reliable sources of information are provided for 48% of the documented events and medium reliable for 52%.

In conducting flood frequency analysis, essential data were sourced from the RHMSS. This dataset comprises the annual maximum discharges recorded at eight gauging stations situated along the rivers in the Šumadija region (

Table 1. Characteristics for observed stations.

River	Gauging Station	Data Period	Elevation m a.s.l	Distance from the Mouth (km)	Watershed Area (km2)	Average Slope (°)
Topčiderska	Rakovica	1963–2022	83.22	5.3	138	7.05
Ljig	Bogovođa	1963–2022	110.07	7.3	679	9.2
Peštan	Zeoke	1963–2005	96.92	3.8	168.3	6.7
Dičina	Brđani	1983–2022	254.6	8.2	208	8.3
Belica	Jagodina	1963–2022	115.34	9.5	193	7.4
Jasenica	Satornja	1969–2022	262.95	65.9	83.6	12.1
Lepenica	Batočina	1975–2022	106.15	4.7	584	6.8
Lugomir	Majur	1975–2022	123.29	9.6	427	10.0

). These records span from 1963 to 2022 and are visually depicted in Figure 1.

One prominent method that emerges as a cornerstone in flood frequency analysis is the Annual Maximum Series (AMS) approach. This method stands as a stalwart in the discipline, drawing upon the foundational principles rooted in the extreme value theorem, as delineated by scholars such as [24,25,26,27]. Through the AMS method, each year forms a distinct block in which the pinnacle of the hydrological events surfaces—the annual maximum discharge. This deliberate selection process allows for a focused examination of the most extreme occurrences, shedding light on the upper echelons of flood magnitudes, a practice elucidated by the works of [28]. Moreover, the AMS method distinguishes itself through its adeptness in navigating critical challenges intrinsic to flood frequency analysis. Issues surrounding the criteria of independence, the judicious selection of thresholds, and the distribution of exceedances are deftly addressed within its framework, a testament to its robustness as underscored by the studies conducted by [29,30].

To ascertain flood frequency at specific sites within the Šumadija region in Serbia, the identification of an appropriate probability distribution is of paramount importance. In this study, we considered five widely used distributions for Flood Frequency Analysis (FFA) at eight gauging stations along the rivers in the northern part of the Šumadija region. In this study, we opted to test the Pearson Type III (P3), Log Pearson Type III (LP3), Gumbel, Generalized Extreme Value (GEV), and Log-Laplace distributions for flood frequency analysis, despite the widespread use of the GEV distribution in such analyses. The choice to explore alternative distributions stems from the known challenges associated with the GEV distribution, particularly in the context of European hydrology. The suitability of the GEV distribution in Europe has been extensively discussed by [31,32]. One crucial limitation of the GEV distribution is its sensitivity to the length of records, as highlighted by [33]. Short records can introduce a strong bias and uncertainty in estimating the shape parameter of the GEV distribution [34,35]. Even when corrected for the effect of record length, the estimate of the shape parameter tends to vary within a narrow range [33]. By examining multiple distributions, including P3, LP3, Gumbel, and Log-Laplace, we aim to provide a more robust assessment of flood frequency in the Šumadija region, considering the limitations associated with the GEV distribution.

The probability density function (pdf) and the quantile function y(F) for each distribution are presented in

Table 2. Probability density and quantiles functions of the probability distributions.

Distribution	Probability Density Function f(y)	Quantile Function y(F)
GEV	$\frac{1}{\alpha }={\left[1-k\left(\frac{y-\mu }{\alpha }\right)\right]}^{\frac{1}{k}-1}exp\left\{{-\left[1-k\left(\frac{y-\mu }{\alpha }\right)\right]}^{\frac{1}{k}}\right\}$	$\mu +\frac{\alpha }{k}\left[1-{(-logF)}^k\right]$
P3	$\frac{1}{{\beta }^{\alpha }\Gamma \alpha }{\left(y-\mu \right)}^{\alpha -1}exp\left\{-\frac{(y-\mu )}{\beta }\right\}$	Explicit analytical form is not available
LP3	$f\left(y\right)=\frac{1}{\sigma \left(y\right)\sqrt{2\pi }}exp\left[-\frac{1}{2}{\left(\frac{y-\mu \left(y\right)}{\sigma \left(y\right)}\right)}^2\right]$	$y\left(f\right)=\mu \left(y\right)+\sigma \left(y\right)x\left[Z+\frac{1}{2}\left(\frac{Z^2-1}{3}\right)\right]$
GUM	$\frac{1}{\alpha }exp\left[-\frac{y-\mu }{\alpha }-exp\left(-\frac{y-\mu }{\alpha }\right)\right]$	$\mu -\alpha \mathrm{l}\mathrm{o}\mathrm{g}(-logF)$
LL	$f\left(y\right)=\frac{1}{2b}exp\left(-\frac{\left|y-\mu \right|}{b}\right)$	y(F) = $\left\{\begin{array}{c}\mu -b\ln \left(2F\right)\\ \mu +b\ln \left(2\left(1-F\right)\right)\end{array}\begin{array}{c}ifF\le \frac{1}{2}\\ ifF>\frac{1}{2}\end{array}\right\}$

for reference. These mathematical formulations are crucial for understanding the characteristics and behaviors of the selected distributions during flood frequency analysis.

Hosking and Wallis (1997) [36] introduced L-moments as linear functions of probability-weighted moments (PWMs) to provide an alternative to conventional moments, enabling the characterization of any random variable Y with an existing mean. Several scientific sources have demonstrated the superiority of the L-moments method over the Method of Moments (MOM) [37]. The computation of L-moments involved estimating order PWMs (βr) using specific formulas, with sample estimators βr for the first four PWMs derived accordingly. Subsequently, the initial four L-moments (λ1 to λ4), depicting the mean, scale, skewness, and kurtosis of the distributions, were calculated through linear combinations of PWMs. Finally, the L-moment ratios τ2, τ3, and τ4 were computed by the definitions outlined by Hosking and Wallis (1993) [38], providing comprehensive insights into the distributional characteristics of the data.

The determination of the optimal distribution function was performed using the Kolmogorov–Smirnov (K-S) and Cramer-von Mises (CvM) tests. The K-S test assesses the consistency of the probability distribution methods by calculating the maximum unconditional deviation (D_max) between the cumulative distribution functions of the theoretical and empirical data. The Cramer–von Mises test evaluates the concordance between empirical and theoretical distributions using the Nω² statistic.

$${\mathrm{D}}_{\mathrm{m}\mathrm{a}\mathrm{x}}=\mathrm{m}\mathrm{a}\mathrm{x}\left|{\mathrm{F}}_{\mathrm{e}}\left(\mathrm{x}\right)-{\mathrm{F}}_{\mathrm{t}}\left(\mathrm{x}\right)\right|$$

$$\mathrm{N}{\mathsf{\omega }}^2=\frac{1}{12\mathrm{N}}+{\sum }_{\mathrm{i}=1}^{\mathrm{N}}{\left[{\mathrm{F}}_{\mathrm{e}}\left(\mathrm{x}\right)-{\mathrm{F}}_{\mathrm{t}}\left(\mathrm{x}\right)\right]}^2$$

These tests were employed to identify the distribution that best fits the empirical data in terms of flood frequency of the Šumadija rivers. The critical values for these tests were considered in the decision-making process [20,39]. For all of this, an algorithm was created in the Python programming language, which also enabled us to identify the most appropriate probability distribution. The statistical indicators of the observed data are presented in

Table 3. The statistical indicators of the data series for the L-moments.

River	Gauging Station	L1 (m3/s)	L2 (m3/s)	L3 (m3/s)	L4 (m3/s)	τ2	τ3	τ4
Topčiderska	Rakovica	23.6	11.8	5.1	4.4	0.500	0.431	0.374
Ljig	Bogovađa	102.5	60.5	28.7	28.6	0.590	0.473	0.474
Peštan	Zeoke	42.3	20.6	5.8	1.5	0.487	0.283	0.073
Dičina	Brđani	53.2	24.3	9.5	5.2	0.457	0.390	0.212
Belica	Jagodina	20.9	12.5	5.1	2.9	0.596	0.407	0.235
Jasenica	Šatornja	20.8	14.7	9.4	6.7	0.707	0.640	0.454
Lepenica	Batočina	38.4	15.9	4.6	4.2	0.415	0.289	0.261
Lugomir	Majur	69.1	39.2	19.7	14.6	0.567	0.502	0.373

for parameter estimation.

3. Research Area

The Šumadija region is one of the central and economically developed areas of Serbia, which covers 9.5% of the territory of the Republic of Serbia. According to morphological characteristics, in Šumadija, two parts can be distinguished—the lower northern and higher southern part. Its terrain expands at the altitude from 30 m a.s.l to 1132 m a.s.l. (Figure 2). From the regional-geographical view, it is bordered by large rivers—Sava and Danube Rivers in the north, Velika Morava in the east, and Zapadna Morava in the south. The southwestern boundary is defined by the rivers Čemernica, Dičina, and Mala Dičina to the top of Rajac Mountain, while the western border is led by the Ljig and Kolubara Rivers. In the hydrological analysis and regional chronology of torrential flood events, we excluded the large border rivers. Among river borders, only the Dičina River is included as suitable according to the subject. The Šumadija area is characterized by diverse and complex geological structures. According to the Geological Map of Serbia [40], in the Šumadija region, the most widespread are Tortonian and Messinian brackish sediments (Neogene), Pleistocene and Holocene sediments (Quaternary) in the valleys of major rivers and flysch, and other Lower Cretaceous basin sediments (Mesozoic) covering smaller mountains Rudnik and Gledićke. The Precambrian and Cambrian metamorphosed sedimentary and igneous rocks cover mountains Juhor and Crni Vrh, while those of Ordovian and Devonian are present on the Bukulja Mountain. The most widespread soils in the Šumadija region are brown forest soils (33.4%), smonitza and metamorphized smonitza (29.1%), and acid brown and podzolic soils (17.6%) [41].

The Šumadija region, likewise the largest part of Serbia, has a temperate continental climate, and the inter-annual distribution of discharges indicates that the water regime on all rivers in Šumadija is pluvial-snow (mostly dependent on rainfall and then snow melting). Though the region was named after the Serbian word for forest (“šuma”), today, the share of the area under forest is only 24.2%, according to the dataset of Corine Land Use Land Cover 2018 [42].

4. Results and Discussion

4.1. Characterization of Torrential Flood Occurrence and Trends

The spatial and temporal characterization of past torrential floods in the Šumadija region is presented by the dataset of 344 documented torrential flood events with over 43 fatalities for the period 1929–2020. The first recorded event in the Šumadija region refers to the torrential floods of the Velika Morava (Lepenica, Jasenica, Lugomir, and Lepenica) and Zapadna Morava (Gruža and Čemernica) tributaries on 16 May 1929. The dataset currently ends with the floods of 23 June 2020 which occurred on the tributaries in the Kolubara and Velika Morava River basins. At that time, the Dičina River in the Zapadna Morava River basin had the most extreme flood—specific maximal discharge of 0.64 m³s⁻¹km⁻². There are some sources of information on torrential floods in the Šumadija for the period before 1929; however, the time gap between the events is even several decades, which is a reason for their omission. Torrential floods in 2021, 2022, and 2023 in the Šumadija area were not recorded.

Out of 344 documented torrential flood events in the Šumadija region, hydrological data are available for 160 events, or 46.5%. The rest are marked as documented with the minimum attributes. In this dataset, there are 33 events with specific maximal discharges above 0.5 m³s⁻¹km⁻² and a few events with above 1 m³s⁻¹km⁻². The two most extreme torrential floods happened on:

• 10 July 1999 in the watershed of Jasenica River (profile: Donja Šatornja), when the specific maximal discharge was 2 m³s⁻¹km⁻², and the mean daily discharge on the previous day was 0.39 m³s⁻¹, and
• 19 February 1986 in the watershed of Lugomir River (profile: Majur), when the specific maximal discharge was 1.01 m³s⁻¹km⁻², and the mean daily discharge on the day before was 1.25 m³s⁻¹.

Unfortunately, meteorological data for the first several decades are scarce—there is only information on the quantity and/or duration of rainfall for just several events. The situation is not much better for the following decades and the main reason for this is the unavailability of data from rainfall stations. RHMS publishes meteorological yearbooks [43] with data from climatological stations annually and bulletins on extreme rainfall episodes sporadically in the last decade. The bulletin for the extreme rainfall episode in May 2014 [44], for example, reveals that the Bukulja, Donja Šatornja, and Rudnik stations recorded 218, 227, and 277 mm, respectively, on three days, 14–16 May.

The question arises as to whether the frequency and peaks of the occurrence of torrential floods in the Šumadija region overlap with the findings of the Inventory level. According to the Inventory of torrential floods, the frequency of their occurrence within a year shows a primary peak in May and June and a secondary peak in April and March, and within 105 years, there is a clear upward trend [4].

The distribution of recorded events per month in the Šumadija region (Figure 3) indicates that the majority of recorded torrential floods occurred in a warmer part of the year—in May and June (74 or 21.5% each), then in July (47 or 13.7%), and March (40 or 11.6%) and April (33 or 9.6%) in a colder part of the year. Consequently, most fatalities are recorded in May (19), July (18), and June (5) in several dozen flood events. This output is explained by the rainfall regime of the temperate continental climate in this part of Serbia, where the major and the highest intensity rainfall episodes occur in May and June, while the sudden snow melt occurs alone or combined with intense rainfall in early spring happens in March and April or even in February. However, the two most extreme runoff peaks in the Šumadija region mentioned above occurred in February and July.

Therefore, the retention capacity of the soils is reduced in early spring, while in the summer months of July and August, sudden flooding may occur after drought periods, disturbing the structure of soils that easily erodes in conditions of intense rainfall events. Some severe floods in July and August occurred on the right tributaries of the Kolubara and Topčiderska reka on 17 August 1956, Topčiderska reka on 30 August 1985, on the torrents throughout Šumadija region on 10 and 31 July 1999, and on the right tributaries of the Kolubara River on 11 August 2004.

The annual distribution of recorded events allows us to distinguish the peak years based on the number of documented events (Figure 4), and these are 1999 (38), 2014 (18), 1981 (16), and 1986 (16). In terms of the death toll, the peak years are 1999 (18 deaths) and 1969 (3), while 1986 and 1994 each record 1 flood death, respectively. The annual average of registered torrential flood events is 6.3 if only the years with recorded events are taken into account. This dataset is 92 years long, with 55 years of recorded torrential flood events. Though Figure 4. shows a linear growth of the number of recorded torrential flood events, if we group them into three periods: 1929–1960, 1961–1990, and 1991–2020, then we obtain the following averages per period respectively: 1.5, 5.4, 4.9. Although the largest number of torrential flood events was expected for the third period, there is a slight decline. The information from the subject literature claims that intensive rainfall showers have become more frequent in recent decades, and there is an increase in mean annual hydrological hazards in the world per decade. Hydrologists all over the world have agreed that discharges with a recurrence period of 100 years shortened their return interval to 50 or 20 years, which increases the importance of torrential flood studies. The reasons for a decline in the number of documented events may also be a lack of hydrological monitoring and other sources of information. The greatest recorded specific maximal discharge in the first period was in the Peštan watershed on 16 May 1955 with q_maxsp = 0.44 m³s⁻¹km⁻², and in the second and third periods, in addition to the two aforementioned absolute extreme specific maximal discharges for the entire dataset, the other excessive maximal discharges are documented on:

• 22 May 1967, with q_maxsp = 0.94 m³s⁻¹km⁻² and 14 June 1969, with q_maxsp = 0.96 m³s⁻¹km⁻², both in the Peštan watershed;
• 16 May 2014, with q_maxsp = 0.96 m³s⁻¹km⁻² in the Ljig watershed, and 13 August 1991, with q_maxsp = 0.91 m³s⁻¹km⁻² in the Dičina watershed. It is worth mentioning that on 15 May 2014, the maximal daily discharge for the Topčiderska River was 255 times larger than the average annual discharge for the period of hydrological monitoring.

The spatial characterization of the torrential flood phenomenon is presented in Figure 5, with mapped locations of the most extreme and destructive torrential flood events with or without fatalities, and

Table 4. Distribution of torrential flood events per river basin.

River Basin	Number of Events	Death Toll
Velika Morava	84	17
Kolubara	53	9
Zapadna Morava	35	2
Danube	13	/
Sava	39	>15
Total	334	>43

. The greatest number of torrential flood events is recorded in the watersheds of tributaries of the Velika Morava due to the highest area share, followed by the Kolubara and the Zapadna Morava. Natural conditions (steep slopes in higher lands and soil properties), but also forest exploitation and destruction and urbanization in the Šumadija region, have led to an intensification of erosion processes and frequent striking torrential floods, which cause heavy damage, including death toll and casualties.

The torrential floods of the left tributaries of the Velika Morava (17 flood deaths) and Topčiderska reka and other nearby tributaries of the Sava River (>15 flood deaths) took the most human lives. Due to the largest flood death toll in the last period, it is important to foster education and communication about risk, prevention, and preparedness and to advance the warning system.

In this sense, the thresholds for torrential flood occurrence should be defined for each watershed, which can be one of the future research directions. When the warning and monitoring systems are adequate, not a single life is lost. However, it is a question of whether politicians, citizens, farmers, and planners learn from past disastrous torrential flood events [6]. It should be borne in mind that the effects of torrential floods are increasingly intensified by anthropogenic factors such as wild dumps in river beds and uncleared river channels.

This dataset includes not only torrential flood events with major damage in urban and suburban areas but also those in the upper part of watersheds with less impact on the rural population. Villages and towns are primarily identified as affected locations, but damaged railways, motorways, and roads are sometimes recorded as the only affected locations. The material loss is not documented for all events. For the historic and extreme torrential floods of May 2014, however, the most detailed report on losses is provided by experts. The assessment of total material loss for the dataset period based on the event documentation could be a task for economic experts in further directions of dataset research.

The main limitation of the dataset is a lack of meteorological data (the data of the rainfall stations of the Republic Hydrometeorological Service are not open access). Then, the torrential flood dataset contains differently detailed event information—while some events are well documented, others provide only basic information, which depends on the available sources. The dataset includes the hydrological data for only eight gauged watersheds. For this regional chronology, further efforts can be made in the future in the direction of reconstruction of past torrential floods and estimation of maximal discharges in ungauged watersheds. The scientific papers aimed at flood simulation by employing the hydrological models, the combined method of the NRCS (Natural Resources Conservation Service of the U.S. Department of Agriculture) and synthetic unit hydrograph, the hydraulic method of flood traces, the rational method or even the dendro-geomorphological method would enhance the dataset documentation.

4.2. Flood Frequency Analysis

FFA is an essential method aimed at understanding and predicting the frequency and magnitude of flood events. This is crucial for effective water resource management and risk mitigation strategies. The fundamental expectations of the FFA and the fitting of probability distributions encompass criteria such as randomness, independence, stationarity, and skewness of the data series [38]. Thus, the initial phase involved computing basic statistics for the AMS data, including mean value, median, standard deviation, skewness, and kurtosis, for six gauging stations. The corresponding values are detailed in

Table 5. Descriptive statistics of eight gauging stations in the Šumadija region.

River	Gauging Station	N	Mean (m3/s)	Standard Error	Median (m3/s)	Standard Deviation	Kurtosis	Skewness	Minimum (m3/s)	Maximum (m3/s)
Topčiderska	Rakovica	44	23.6	3.7	15.86	24.9	0.7	2.9	1.6	138
Ljig	Bogovađa	57	102.5	12.1	85.9	91.2	0.8	4.3	24.5	651
Peštan	Zeoke	47	42.8	5.4	27.4	37.3	0.6	0.9	1.48	120
Dičina	Brđani	36	53.2	8.0	35.7	48.0	0.7	1.6	7.6	190
Belica	Jagodina	53	20.94	2.8	12.8	20.2	0.6	1.5	1.53	87
Jasenica	Šatornja	45	20.8	5.6	6.5	37.4	0.3	3.3	1.3	175
Lepenica	Batočina	41	38.4	5.1	31.7	32.4	0.8	2.9	6.4	193
Lugomir	Majur	44	69.1	12.8	41.3	84.7	0.6	2.6	3.2	431

.

The results presented in Table 3 provide a comprehensive overview of the basic statistics for the discharge data of selected gauging stations. These statistics contribute to a nuanced understanding of the hydrological characteristics at each gauging station, which is essential for effective flood frequency analysis.

The goodness-of-fit tests, consisting of the Kolmogorov–Smirnov (K-S) and the Cramér–von Mises (CvM) tests, were employed to evaluate the appropriateness of diverse distribution functions for characterizing flood frequency in the rivers of the Šumadija region (see Appendix A). At the Batočina gauging station, the Gumbel distribution emerged as the most suitable, with a K-S test statistic of 0.097 and a p-value of 0.799, reinforcing its robust fit. The CvM test also favored the Gumbel distribution, recording a low statistic of 0.053 and a high p-value of 0.861. For further visual confirmation of these results, Cumulative Distribution Function (CDF) graphs were created, providing visual validation of the statistical outcomes (see Appendix B). In Bogovađa, the Log-Laplace distribution demonstrated superior performance in both the K-S and CvM tests, showcasing minimal statistics and highly significant p-values. Brđani’s results further corroborated the excellence of the Log-Laplace distribution, revealing lower statistics and comparable p-values for both tests. Meanwhile, in Jagodina, the Log-Pearson 3 distribution prevailed in the K-S test, whereas the Generalized Extreme Value distribution excelled in the Cramér–von Mises test. Majur consistently favored the Log-Laplace distribution, exhibiting lower statistics and significant p-values in both the K-S and Cramér–von Mises tests. At Rakovica, the Log-Laplace distribution consistently outperformed other options, showcasing lower statistics and significant p-values in both tests. Šatornja’s results indicated the Generalized Extreme Value distribution as the best fit, recording lower statistics and significant p-values in both tests. Lastly, at Zeoke, the Pearson 3 distribution gained prominence, demonstrating lower statistics and significant p-values in both the K-S and CvM tests. Overall, the choice of the most suitable distribution varied across locations, emphasizing the necessity of considering both statistical measures and practical implications for comprehensive flood frequency analysis in the Šumadija region.

Estimating discharge values associated with high return periods faces a significant challenge due to uncertainties arising from both the length of the data series and the chosen distribution shape parameter [45]. Therefore, caution is essential, particularly when estimating discharges for extended return periods, such as 1000 years [20]. That is why, in our analysis, we only considered return periods up to 500 years. It is worth noting that regional FFA often combines data from multiple rivers within a region to derive generalized estimates. However, our approach focused on individual rivers in the Šumadija region, recognizing the unique hydrological characteristics of each. This strategy allowed us to capture localized variations in flow patterns, hydraulic regimes, and watershed properties that may be overlooked in regional analyses. The calculated discharges for various return periods, with 95% confidence intervals, provide crucial insights into the potential magnitudes of river flows under extreme conditions (Figure 6). The confidence intervals for each distribution were determined as depicted in [27,37]. For the Topčiderska River at the Rakovica gauging station, the estimated discharges increase substantially with longer return periods, reaching values of 48.83 m³s⁻¹ for a 10-year return period and escalating to 829.05 m³s⁻¹ for a 500-year return period. Similarly, the Ljig River at Bogovađa exhibits escalating discharges, ranging from 175.94 m³s⁻¹ for a 10-year return period to 1002.9 m³s⁻¹ for a 500-year return period. The Peštan River at Zeoke, the Dičina River at Brđani, and the Belica River at Jagodina display similar trends, with discharge values ascending notably with increasing return periods. The Jasenica River at Šatornja showcases particularly substantial discharges, with values soaring from 50.37 m³s⁻¹ for a 10-year return period to a remarkable 4038.87 m³s⁻¹ for a 500-year return period. Lastly, the Lepenica River at Batočina and the Lugomir River at Majur also manifest increasing discharges, providing essential data for well-founded decision-making in flood risk management and infrastructure planning. The findings of this study demonstrate the effectiveness of the employed technique in accurately estimating discharge values for flood peaks associated with specific return periods, particularly in watersheds characterized by incomplete data or a limited temporal record.

5. Conclusions

The first part of this research aims to streamline the process of documenting torrential floods in the Šumadija region in Serbia. This is a significant contribution to the implementation of the European Union Inspire Directive (data theme Natural risk zones), the European Union Flood Directive, and the management of torrential flood risks in Serbia. The historical context of the presented regional chronology of torrential floods in the Šumadija may be of great help to municipal authorities and insurance companies. Thus, this research serves for a better understanding of the torrential phenomenon and fosters further torrential flood research. The dataset of torrential flood events in Šumadija, with 344 recorded torrential flood events and more than 43 fatalities, is presented for the period 1929–2020. The intra-annual peaks are distinguished, and the spatial pattern of this natural phenomenon is characterized.

Furthermore, this study undertakes a comprehensive flood frequency analysis for the rivers in the Šumadija region. Utilizing a diverse array of probability distributions, such as Pearson 3, Log Pearson 3, Gumbel, Generalized Extreme Value, and Log-Laplace. The analysis incorporated goodness-of-fit tests, such as the Kolmogorov–Smirnov and Cramér–von Mises tests, to identify the most suitable distribution for each monitoring station. The results revealed that, based on both tests, the Log-Laplace distribution consistently outperformed other distributions, providing the best fit for several stations, including Bogovađa, Brđani, Majur, and Zeoke. These outcomes were further confirmed through Cumulative Distribution Function graphs, detailed in the appendices. Additionally, the calculated discharges for different return periods were presented, offering valuable insights into extreme flow events. The study’s contribution lies in its rigorous examination of distribution fitting methods and their application in flood frequency analysis for the region, providing crucial information for effective water resource management and infrastructure planning in the face of changing hydrological patterns.