How Frequent Is an Extraordinary Episode of Precipitation? Spatially-Integrated Frequency in the Júcar-Turia System (Spain)

1. Introduction

1.1. Context

Extreme precipitation constitutes one of the primary natural hazards in the coastal regions of the western Mediterranean. In particular, the Júcar-Turia system (Spain) is among the most exposed areas in Europe to torrential rainfall episodes, with historical records frequently approaching or exceeding 500 mm in less than 24 hours [1,2]. This hazard arises from quasi--stationary mesoscale convective systems (MCS) as a result of the interaction of: i) a massive supply of warm, moist air from the Mediterranean Sea, ii) an upper-level flux divergency due to the proximity of an isolated depression (DANA or cut-off low), and iii) a complex orography that facilitates an efficient humidity convergence. The confluence of all these factors have driven recurrent catastrophic events in the Júcar-Turia system, mainly in October, including those of 1957 (València), 1982 (Tous), 1987 (Oliva) and the episode of 29 October 2024, which reached 771.8 mm in 24 hours at Turís [3], the second highest 24-hour rainfall total recorded for Spain, after 817.0 mm observed in Oliva on November 3, 1987. Such phenomena lead to devastating flash floods, substantial economic losses, and high mortality rates, underscoring the increasing systemic risk in the region.

In the context of climate change, the intensification of these events has been widely documented [3,4,5,6,7]. Both regional scenarios and observational studies indicate an increasing intensity of rainfall up to +20% in areas that are already highly vulnerable, such as the Valencian coastline [3,8]. Therefore, the study of the temporal structure and statistical recurrence of these episodes is essential for risk management, hydrological planning, and territorial adaptation.

1.2. Episode Recurrence

Numerous studies have addressed the characterization of these phenomena from multiple perspectives. At the synoptic scale, atmospheric patterns associated with extreme events are classified through statistical circulation analysis, confirming the decisive role of blocking configurations and retrograde flows [9,10,11]. On the regional scale, studies such as González-Hidalgo et al. [2] have cataloged a very large number of extraordinary precipitation episodes in the Iberian Peninsula since the beginning of the twentieth century, while Ruiz et al. [12] investigated palaeofloods in the Júcar–Turia system through geomorphological analyzes of alluvial fans and fluvial plains, with at least 10 major floods in five centuries (1571, 1590, 1632, 1776, 1779, 1805, 1864, 1897, 1923, 1949 and 1957). These reconstructions reveal that the actual frequency of catastrophic events far exceeds estimates based solely on instrumental meteorological records [13,14].

An important limitation repeatedly identified in the hydrological literature is the excessive reliance on the classical micro-catchment or point-based approach of Intensity-Duration-Frequency (IDF) curves [15]. Here, micro-catchments are understood as basins of small spatial extent (typically tens to <100 km²), with short concentration times (hours) and rapid hydrological responses strongly controlled by local orography and drainage structure. These units exhibit pronounced spatial and temporal variability in rainfall distribution and runoff generation: the same convective event may produce markedly different maxima at stations separated by only a few kilometers. This spatial heterogeneity is a pattern widely documented in recent studies [16,17].

The point-based (micro-catchment) approach—summarizes hazards using a single representative value per basin—fails to capture this behavior and, critically, the spatial probability of exceedance, i.e., the likelihood that a given threshold is exceeded in at least one subregion of the domain. Recent studies show that extreme precipitation events often display complex spatial dependence structures whose extent and intensity evolve with climate forcing [18], and that neglecting spatial variability leads to systematic underestimation of flood hazard [17]. Moreover, the spatial dependence of extreme rainfall can substantially modify recurrence intervals in a basin, making the analysis of a single station insufficient for regional hazard assessment [19].

Therefore, since basin- or regional-scale hazard depends fundamentally on the joint spatial behavior of extremes, it is necessary to complement (or replace) purely point-based estimates (e.g., Gumbel fitted to a single station or micro-catchment) with explicitly spatial methodologies that integrate information from multiple stations and model the spatial dependence of extremes. The continued use of point-based representations assumes a homogeneous, station-equivalent behavior that is incompatible with the spatial complexity of extreme precipitation processes under current and future climate conditions [20]. As a consequence, the erroneous idea has become widespread that extreme events such as those mentioned above have return periods of hundreds or thousands of years, whereas empirical, geological, and statistical records indicate that they recur several times per century [12,14].

1.3. Rainfall Concentration

In addition, recent studies highlight the need to use new statistical indicators that capture the temporal variability of precipitation beyond simple accumulations. Among them, the fractal intermittency n-index proposed by Monjo [21] stands out as a dimensionless measure that describes the structure of precipitation events according to their degree of concentration and intensity across time scales. Specifically, the maximum precipitation amount ($P_t$) expected in a duration of t hours is given by

$$P_t=P_{1h}{t}^{1-n}⟹P_{\mathrm{tot}}=P_{1h}{t_{\mathrm{ef}}}^{1-n}$$

(1)

where $P_{1h}$ is the maximum precipitation recorded in 1 hour, and the effective duration $t_{\mathrm{ef}}$ is defined when the maximum amount is equal to the total amount. Eq. 1 can be derived from the Maximun Average Intensity (MAI) relationship implemented in empirical IDF curves with power-law behavior [15]. The n-index is a dimensionless indicator that describes precipitation variability on different temporal scales. Its value ranges between 0 and 1, where values close to 0 indicate more uniform–continuous rainfall, whereas values close to 1 indicate highly concentrated–intermittent rainfall, and $n=0.5$ denotes a balance between intensity and duration [21].

The n-index approach has been applied in both global and regional analyzes and has proven useful for complementing the hydrological modeling of intense episodes [16,22,23,24,25,26]. Its application is key in climate modeling, allowing a better understanding of rainfall patterns and their impact on water resource management and the prediction of extreme events.

Against this background, the present work proposes an advanced statistical approach to characterize extreme precipitation episodes in the Júcar-Turia system, using a combination of high-resolution data and long-term historical series. The analysis begins with a detailed examination of the event that occurred on 29 October 2024, one of the most devastating in recent decades, which particularly affected the municipality of Turís (Valencia Province).

Through the use of the fractal–intermittency n-index and the estimation of empirical and theoretical return periods through fitting to different probability distributions, the representativeness and frequency of this episode are evaluated within the context of the last decades. Unlike previous studies, the approach adopted here integrates multiple stations within each water basin, allowing the analysis of annual maxima and their aggregated spatial behavior. In addition, a sensitivity analysis is introduced as a function of the number of stations to assess how return periods and fitting errors vary as greater spatial coverage is incorporated.

The article is structured as follows: Section 2 presents the data used and describes the methodological framework, including the calculation of spatially-integrated return periods with sensitivity analysis and the n-index; Section 3 presents the results obtained; while Section 4 discusses their climatic, hydrological, and territorial interpretation, and draws general conclusions to outline future research directions.

2. Materials and Methods

2.1. Observed Data

A total of 98 time series of daily precipitation (1981–2024), recorded in stations located in the Júcar-Turia system, were collected from the Spanish Meteorology Agency (AEMET) database: 58 stations correspond to the Júcar basin and 40 are located in the Turia basin (Figure 1). To analyze the extreme precipitation episode that occurred on 29 October 2024 and severely affected Valencian territory, 10-minute precipitation data recorded by an automatic weather station located in the municipality of Turís were also used (records from 28/10/2024 00:00 h to 31/10/2024 00:00 h, local time).

2.2. Point and Spatially-Integrated Frequency

Point and spatially-integrated return periods were estimated for several reference thresholds of extraordinary episodes, including the 2024 case. Point return periods ($\rho$) represent the expected average time between precipitation (p) events that exceed a certain threshold of pluvial amount (P, our case), intensity, or fluvial avenue, according to the cumulative probability $\Pi (p\le P)$ of a unique time series ($N=1$); that is,

$$\rho \equiv {\displaystyle \frac{1}{\Pi (p>P)}}={\displaystyle \frac{1}{1-\Pi (p\le P)}}.$$

(2)

In contrast, spatially-integrated periods also represent return periods $\rho$ but using $N>1$ time series (from stations located nearby), so the probability of exceeding thresholds is significantly higher. Both indicators are used to assess the probability of occurrence of extreme events over a given time span. The longer the return period, the less frequent and more intense the associated event.

The time series used for all subsequent calculations consist of the annual maximum values extracted from a set of N time series (subsequently adding stations); that is, the maximum of the annual maxima is considered to build a synthetic annual time series that integrates all the N stations considered. Distribution fitting was performed using families commonly employed in extreme value theory; Weibull, Gumbel, Generalized Extreme Value (GEV), Gamma, and Pareto distributions were considered [27], as these families are specifically designed to model distribution tails and extreme phenomena. In contrast, central-limit distributions such as the normal distribution are not suitable for describing the probability of occurrence of extreme values due to their rapidly decaying tails. The resulting areal IDF curves were compared to the empirical n-index power law found by Moncho and Caselles [15]:

$$P_{t,\rho ,N}=P_{t_0,{\rho }_0,N_0}{\left({\displaystyle \frac{\rho N}{{\rho }_0{N}_0}}\right)}^m{\left({\displaystyle \frac{t}{t_0}}\right)}^{1-n}⟹\rho \sim 1/N,$$

(3)

where m is the scaling parameter, while $P_{t,\rho ,N}$ and $P_{t_0,{\rho }_0,N_0}$ are, respectively, the maximum precipitation amounts recorded in a small area by N and $N_0$ stations, as expected for a return period $\rho$ and ${\rho }_0$, within an episode of duration t and $t_0$. If $N=N_0=1$, Eq. 3 provides a classical empirical IDF curve as $I_{t,\rho }=P_{t,\rho }/t$. On the other hand, for $N\ge N_0>1$ statistically-independent stations, Eq. 3 represents an empirical relationship between the spatially-integrated return period $\rho$ and the maximum precipitation amount $P_{t,\rho ,N}$. To analyze the effect of N on the return period $\rho$, we set the curves of Eq. 3 for constant precipitation thresholds $P_{t,\rho ,N}=P_{t_0,{\rho }_0,N_0}$, for example, 200 or 500 mm. Here, we assumed that the return period changes as follows:

$$\rho ={\rho }_{min}+({\rho }_{max}-{\rho }_{min}){\left({\displaystyle \frac{1}{N}}\right)}^{\lambda \left(N\right)},$$

(4)

where ${\rho }_{max}$ is the return period considering only the extremest station, while ${\rho }_{min}$ is the lowest return period, found when all the stations are used, and $\lambda \left(N\right)$ is an independence parameter: Statistically, $\lambda =1$ means that the stations are purely independent (i.e., recovering the $\rho \sim 1/N$ relationship of Eq. 3 with ${\rho }_{min}=0$), while $\lambda =0$ implies that all stations are purely dependent and do not contribute to new extremes. For this study, we considered a reduction of $\lambda$ from 1 to 0 as $N/N_0$ increases with respect to a certain $N_0$, specifically $\lambda \left(N\right)\approx exp(-N/N_0)$.

2.3. Sensitivity Analysis of Episode Recurrence

The fit of the theoretical cumulative distributions was compared with the empirical cumulative function to evaluate which model best reproduces the observed behavior across different ranges (e.g., fitted versus extrapolated). To differentiate between well-sampled and extreme cases, two disjunct ranges of return periods were considered: the fitting range (40 values) and the validation range (3 extremest values). The fitting range includes return periods where multiple observations are available, allowing for robust statistical fitting. In contrast, the validation range corresponds to very large return periods, often represented by only one or two observations. These extreme cases are exceptionally rare, leading to return periods that are extremely high and potentially distant from the true values due to the lack of additional extreme observations. Evaluating the validation range separately helps quantify the uncertainty associated with extrapolating beyond the observed data.

The accuracy of the fitted distributions was quantified using the mean absolute error (MAE) and the root mean square error (RMSE), calculated between the modeled and empirical values of the spatially-integrated return periods. Given a set of K observations (i.e. K annual maxima and their return periods), the errors are defined as:

$$\mathrm{MAE}={\displaystyle \frac{1}{K}}\sum _{i=1}^K\left|{\rho }_i^{\mathrm{model}}-{\rho }_i^{\mathrm{empirical}}\right|$$

and

$$\mathrm{RMSE}=\sqrt{{\displaystyle \frac{1}{K}}\sum _{i=1}^K{\left({\rho }_i^{\mathrm{model}}-{\rho }_i^{\mathrm{empirical}}\right)}^2},$$

where ${\rho }_i^{\mathrm{model}}$ and ${\rho }_i^{\mathrm{empirical}}$ are the modeled and observed return periods, respectively. These metrics were calculated separately for the fitting range and the validation range, allowing assessment of model performance in well-sampled versus extreme, poorly-sampled conditions.

The spatial-integration analysis was then conducted by computing the evolution of the return period as stations are progressively added to the analysis, in order to capture a spatial evolution rather than a purely temporal one. Calculations were performed on maximum series derived from different station subsets, with blocks of five stations. For each subset, the different fitting methods were compared in order to select the one with the lowest error.

2.4. Rainfall-Concentration Comparison

To complete the time structure of the obtained IDF curves, a comparison of different extreme precipitation episodes was made together with the data for the 29 October 2024 episode using the fractal–intermittency n-index (Eq. 1). It is calculated from the multifractal analysis of the precipitation series, using a multiplicative cascade approach to assess variability at different scales [21]. The analysis of rainfall concentration was carried out mainly using the n-index.r script and functions_spell.r scripts [28] to compute the n-index and display the results of the case study. For the estimation of return periods, both fitdistrplus and evd libraries were applied [29,30].

For general data processing and analysis, RStudio was used, an integrated development environment for the R programming language that facilitates data manipulation, statistical model development, and result visualization, providing advanced tools for the exploration and analysis of precipitation time series,

3. Results

3.1. Episode Recurrence

For the Júcar basin (Figure 2 a), the best fit is obtained using either the Gumbel or the GEV distribution, although it should be noted that for very high precipitation values none of the methods is able to adequately reproduce the observed maximum. In contrast, for the Turia basin (Figure 2 b), the Pareto distribution appears to provide the best fit to the empirical data, except for the highest precipitation values (validation range), where Weibull or GEV show a better performance. To more precisely assess which fit performs best, fitting errors were calculated with respect to the empirical distribution (Table 1).

After computing the errors, Pareto is identified as the best fit for the Turia basin, except for the values used for validation, for which Weibull provides the best performance. For the Júcar basin, Gumbel exhibits the lowest error across most of the distribution, except for the most extreme values, where GEV better defines the return periods.

Although the analysis was carried out using the best-performing method for each individual case (and even for each subset of stations), overall, GEV produced the lowest average error. Finally, for stations located within the Júcar basin, a GEV fit was applied to all subsets due to the small differences in error in the estimation of the return period compared to the large differences observed in return periods for extreme values. In contrast, for the Turia basin, GEV was the dominant method for the first two station subsets, whereas Weibull was applied to the remaining subsets.

Figure 3. Evolution of the spatially-integrated return periods for Júcar basin (a) and Turia basin (b) according to distinct severity thresholds (four colors). Intervals represent the MAE error of the fitted return periods. Dashed lines represent fitted curves (Eq. 4) with three parameters (Table 2).

Figure 3. Evolution of the spatially-integrated return periods for Júcar basin (a) and Turia basin (b) according to distinct severity thresholds (four colors). Intervals represent the MAE error of the fitted return periods. Dashed lines represent fitted curves (Eq. 4) with three parameters (Table 2).

The results showed that the evolution of the return period follows a power-law decay curve with a final stabilization (Figure 4). This structure is more stable for the Turia basin, as there is less variation between extreme values, while for the Júcar basin the evolution for the highest thresholds (e.g. 400 and 500 mm) is more irregular, as the incorporation of many stations ($>30$) with less extreme rainfall reduces the severity. It should also be noted that the evolution of the fitting errors is proportional to the number of stations considered: the greater the number of stations, the larger the uncertainty in the return periods.

For both basins, the five or ten stations that exhibit the highest maxima are those that logically introduce the largest reduction in the estimates of the spatially-integrated return period, with the station that records the highest precipitation maximum in each basin standing out most clearly. Table 2 confirms that a critical number of stations is $N_0\approx 9$ for Júcar and $N_0\approx 3$ for Turia, which provides most of the reduction in the return period (from ${\rho }_{max}$ to ${\rho }_{min}$).

Moreover, if $N=N_0$ is fixed, Eq. 3 provides a representative areal IDF curve such as

$$P_{1day,\rho }=P_{1day,1year}{\left({\displaystyle \frac{{\rho }_{min}}{1\mathrm{year}}}\right)}^m$$

(5)

with scaling parameter $m=0.31\pm 0.02$ for Júcar and $m=0.30\pm 0.09$ for Turia basin, which reproduces the thresholds ($P_{1day,\rho }$, Table 2) from $P_{1day,1year}=171\pm 8$ mm and $P_{1day,1year}=134\pm 10$ mm, respectively.

3.2. Rainfall Concentration

For the 2024 catastrophic event, two major precipitation peaks were observed, with the second peak clearly standing out from the first; therefore, the main analysis focuses on this latter segment. An n-index of 0.27 and an effective duration of 4.94 h were calculated, considering that the maximum hourly precipitation reached 170.2 mm and that the maximum accumulated precipitation within the segment corresponding to the largest peak amounted to 620.6 mm (Fig).

The effective duration was also computed considering the entire precipitation time series, while maintaining an n-index of 0.27, as this value best represents the entire series. In this case, the expected duration is 7.3 h. If n is replaced by 0.5 to evaluate the duration that would be required in an “efficient” case for the same total amount of precipitation, an effective duration of 13.3 h is obtained.

To compare with other cases, the n-index and effective duration were estimated according to Eq. 1 and shown in Table 3. The n-index of Turís during the 2024 DANA episode is the second lowest record for extreme precipitation after the Sueca event of 23 September 2008 (326 mm in 3 hours and 144 mm in 1 hour), with $n=0.14\pm 0.01$ for an effective duration of 2.9 hours. The difference can be characterized by the spatial scale of the storms. The anchored convective cell of the 2008 episode (horizontal scales of 5–10 km) was four times smaller compared to the quasi-stationary MCS of the 2024 DANA episode (storm focus size of about 20-40 km).

Assuming stationarity (i.e. omitting climate change effects) and taking into account the parameters of Eq. 5 and Table 3 in Eq. 3, total amount (772 mm) of the 29 October 2024 episode in Turís is given by the following curve:

$$772\mathrm{mm}=(171\pm 8)\mathrm{mm}{\left({\displaystyle \frac{{\rho }_{min}}{1\mathrm{year}}}\right)}^{0.31}{\left({\displaystyle \frac{t_{ef}}{7.3\mathrm{hours}}}\right)}^{1-0.27}$$

(6)

where the effective duration of the 1-day accumulation was $t_{ef}=7.3$ hours. Therefore, the return period of this catastrophe is obtained from Eq. 6 as ${\rho }_{min}=130\pm 40$ years for the Júcar basin, while the threshold of 500 mm has only $33\pm 2$ years of return period (Table 2).

4. Discussions and Conclusions

4.1. Spatially-Integrated Return Periods

This work was conceived as a feasibility assessment of an explicitly spatial approach to estimate return periods of extreme precipitation, progressively incorporating stations within each basin (Júcar and Turia). The results confirm that the approach is operational and consistent with probabilistic reasoning: as the number of stations increases, the return period (when at least one station exceeds a given threshold within the domain) decreases up to a certain spatial sample size of about 15-20 stations. Moreover, beyond that sample, the spatially-integrated return period stabilizes into a plateau (i.e., a minimum return period ${\rho }_{min}$), because newly added stations no longer contribute annual maxima higher than those previously observed. This decreasing power law (e.g., Eq. 4) is particularly clear in the Turia basin and somewhat more irregular in the Júcar basin for high thresholds, suggesting differences in internal heterogeneity between domains. Overall, the evidence supports the idea that a highly localized anomaly should not inflate regional return periods when the objective is basin-scale hazard and when the metric is correctly defined in spatial terms, so with areal IDF curves. Moreover, the value of the plateau ${\rho }_{min}\sim$ 33-130 years for catastrophic thresholds (e.g. 500–700 mm/day) is consistent with the observed palaeofloods in the Júcar–Turia system: at least 10 major floods in five centuries [12].

The observed uncertainty of the estimates is also consistent with the spatial variability and the aggregation effects. As stations are added, estimation errors initially increase: expanded spatial coverage introduces heterogeneity (different local regimes), and spatial dependence among series reduces the effective sample size, so uncertainty does not decrease mechanically with the number of stations. Once the regional tail becomes “well sampled”, both return periods and errors tend to stabilize. This interpretation is consistent with the reported pattern that “the probability is proportional to the number of stations”, as well as with the plateau observed in the return period–number of stations relationship for different precipitation thresholds.

4.2. Selection of Theoretical Distributions

Regarding tail modeling, the study documents differences by basin and range. In the Júcar basin, Gumbel exhibits the lowest mean error across the distribution, although GEV better represents the extremes, which is why GEV is used for the spatial station sets in this basin. In the Turia basin, Pareto provides the best overall fit except in the validation range, where Weibull improves performance; for block aggregation, GEV dominates in the first two station sets, while Weibull is applied to the subsequent ones. On the aggregated scale, the GEV yields the lowest mean error, but a case-by-case criterion is adopted, prioritizing the best-performing fit for each station set. This result reinforces two key ideas: (i) the need for adaptive selection of the distribution family depending on the domain and the range considered (observed versus extrapolated), and (ii) the importance of reporting uncertainty bands specific to the fitted model used in each range.

Finally, several limitations and avenues for improvement, already implicit in the study design, should be acknowledged: (a) the use of annual maxima facilitates comparability but may underutilise POT (peaks-over-threshold) information and seasonality; (b) spatial dependence among stations suggests the use of resampling techniques or spatial extreme-value frameworks to better quantify uncertainty; and (c) the order in which stations are incorporated could be explored through sensitivity analyses (e.g., starting with stations exhibiting the highest maxima or based on climatic proximity criteria). Nevertheless, the primary objective is fulfilled: the proposed spatial method is viable, behaves as expected, and produces stable and traceable return period estimates once the relevant climatological heterogeneity of each basin is adequately covered. In any case, the annual maximum approach is enough for the areal IDF analysis. More sophisticated options such as three- and four-parametric distributions would be more adequate for daily precipitation modeling from entire time series [27].

4.3. Episode Time Structure

The use of the fractal–intermittency n-index provides physical–statistical context to the analysed episodes (e.g., the 2024 DANA event in Turís), allowing comparisons of temporal efficiency and intermittency against historical cases (such as the 2008 event), and linking temporal structure to accumulated severity. The IDF combination of n-index, annual maxima, and tail fitting offers a consistent interpretation of hazard: episodes with values close to $n\approx 0.3$–$0.5$ (more efficient rainfall) can concentrate large amounts within relatively short effective windows, which is highly relevant for threshold definition and for the interpretation of regional maxima employed in spatial return period estimation.

4.4. Final Recommendation

From an operational perspective, the study supports three recommendations: (i) Explicitly define the spatially-integrated return period (e.g., “at least one station”, k-of-N, or areal average) according to the management question, as values are not interchangeable across definitions. (ii) When aggregating stations in blocks (here, groups of five), jointly assess return periods and errors, and stop aggregation when changes in return periods fall below a practical threshold (e.g., $<5$ %) and the error reaches its plateau. (iii) Accompany each return period estimate with its associated uncertainty and with the fitting method selected for that specific station set and range, in line with the distribution-specific error evidence for each basin. (iv) Analyze the dependence of the spatially-integrated return period and the number of stations considered using an empirical relationship (e.g. Eq. 4), and finally obtain a better estimate of the annual maxima ($P_{1day,1year}$, Eq. 5) and its corresponding IDF curve (e.g. Eq. 6).