Assessing Climate Change Impacts on Combined Sewer Overflows: A Modelling Perspective

1. Introduction

A number of latest studies conducted by the Intergovernmental Panel on Climate Change (IPCC, 2014) reported evidence for changes in the occurrence of certain meteorological events during the 21st century. Regarding precipitation extremes in Europe, increased precipitation intensity in a future climate was one of the model results of most recent studies [1,2,3,4,5,6]. The IPCC assessments relied heavily upon the use of Global Circulation Models (GCMs), namely large mechanistic models for the global atmosphere, which are indeed coupled atmosphere-ocean models. Regional Climate Models (RCMs) are forced by the output of the GCMs as initial or boundary conditions and simulate regional climate with a finer resolution, enhancing the spatial resolution of climatic projections. During the previous years, different RCMs have been used to produce high-resolution climate scenario calculations on the European scale (e.g EU-projects PRUDENCE, ENSEMBLES and CORDEX), including the entire Mediterranean basin [7,8,9].

Future climate projections produced by climate models (GCMs and RCMs) are still subject to large uncertainties possibly originated from three major sources such as the emission uncertainty, the model uncertainty reflecting the limited understanding of the different atmospheric processes and their representation in the climatic models, as well as the uncertainty due to natural variability [10]. Especially when using RCM data, the drivers of uncertainty increase, as the model resolution, the forcing lateral conditions and physical parameterizations [11,12] induce more bias and variance in the results. Therefore, the aforementioned uncertainty sources should be considered when assessing the impacts of climate change on different domains, for the estimation process to be more reliable and robust. Model uncertainty can be considered by creating a multi-model ensemble, consisting of climate models that are ideally structurally independent. Kendon et al. [13] combined climate models through a multi-model ensemble to produce more reliable and consistent climate projections, suggesting the use of minimum three components to estimate changes in precipitation extremes.

Different scales (spatial and temporal) between RCMs and local measuring gauges or local urban drainage systems can have significant effects and produce inaccuracies on the precipitation process or the different sewer processes [14,15]. The aforementioned arguments may necessitate the use of dynamic or statistical downscaling methods to reduce the effects of the coarse scales of climate models to precipitation and urban drainage studies. Statistical downscaling includes methods that are quite easy to apply, are not time consuming and are considered rather efficient. This category of methods includes statistical models to transfer the “predictor” variables, which are available on a coarse grid, to the “predictand” variable, corresponding to finer spatial scales. Statistical downscaling methods include the construction of empirical transfer functions, weather typing models or stochastic modelling methods [16]. The former methods create and utilize transfer functions between the coarse scale variable of the climatic model and the smaller scale variable. Bias correction methods include the development of transfer functions between the cumulative distribution functions of the climate model and observed data. Bias represents the error component of the climate model that is independent of time [17]. Among the bias correction techniques, recent studies mainly using precipitation and temperature data, indicated the quantile mapping methods as the most efficient, even for the most extreme part of the distribution of the studied variables [18,19,20,21]. Gudmundsson et al. [22] reviewed and classified bias correction methods using observations of a large number of precipitation stations in Norway achieving to systematically reduce biases in RCM precipitation projections. Teutschbein and Seibert [23] reviewed and evaluated different bias correction methods for hydrological climate change impact studies. Tani and Gobiet [20] developed an extended version of the quantile mapping approach to provide more reliable bias-corrected precipitation extremes by combining the advantages of nonparametric methods and extreme value theory. Chandra and Papalexiou [24] introduced a semi-parametric approach for correcting daily precipitation for bias, providing good results for extreme quantiles. Holthuijzen et al. [21] proposed a hybrid bias-correction method for precipitation, including a robust linear correction for extremes.

Another significant issue that should be considered refers to the fact that hydro-meteorological variables, such as precipitation, resulting from climate models, are not in general available in time series with fine temporal resolution (e.g., 10 min to 1 h time steps). However, these time scales are essential for studying and estimating the effects of climate change especially on medium or small catchment areas. Temporal downscaling and temporal disaggregation are used to overcome such a difficulty. Temporal downscaling usually refers to the generation of data of high temporal resolution by means of statistical techniques, most commonly stochastic models, calibrated using information on the statics of data from lower resolution temporal scales. The most commonly used downscaling approaches are statistical downscaling methods. Nguyen et al. [25] proposed a statistical temporal downscaling method which is based on the scaling Generalized Extreme Value (GEV) distribution to describe the relationships between daily and subdaily extreme precipitation. Terti et al. [26] evaluated the suitability of the abovementioned approach for determining the distribution of short-duration extreme rainfall in Thessaloniki, Greece. Ghanmi et al. [27] used the simple scale invariance concept to statistically downscale daily rainfall in northern Tunisia to subdaily amounts. Yeo et al. [28] used GEV scaling to estimate subdaily annual maximum precipitation from daily values and applied their methodology to available data from Quebec, Canada and Seoul, South Korea. Galiatsatou and Iliadis [29] applied a simple scaling approach based on the GEV distribution to assess short-duration rainfall extremes at ungauged sites transferring information from gauged ones with a relatively homogeneous extreme rainfall climate. Iliadis et al. [30] used a GEV scaling approach to daily rainfall extremes assessed from a long time series, to estimate fine temporal scale extreme events in Thessaloniki, Greece.

Intensity-Duration-Frequency (IDF) curves based on local precipitation measurements are currently utilized for engineering design and management applications, such as flood risk protection structures and infrastructures or flood mitigation projects. The IDF curves are constructed for different return periods representing the variation of rainfall intensity with duration. Theoretical probability distribution functions are fitted to annual maximum rainfall intensities of particular durations ranging from shorter ones e.g., 5 min to daily precipitation events. These curves are currently designed under the assumption of stationarity [31,32,33], namely the hypothesis that the occurrence probability of precipitation events will remain unaltered or will not exhibit significant changes in time, mainly due to the lack of predictability on details of possible future changes. But the fact is that such hydro-meteorological signals, especially at their extreme levels, exhibit phenomena of nonstationarity. Natural climatic variability, human interventions in the hydrologic cycle of different catchments areas and climate change are some of the prominent causes of such nonstationarities [34,35]. Yan et al. [36] reviewed literature on updating IDF curves under climate change, mainly emphasizing on covariate-based distributions for modelling extremes. Martel et al. [37] studied changes in extreme precipitation and related sources of uncertainty. Schlef et al. [38] presented a review on current knowledge and research on nonstationary IDF curves. Nonstationary IDF curves are usually constructed under the assumption of covariate-based modelling or simulated precipitation from RCMs or GCMs. In the latter case, precipitation extremes are obtained over historical and future periods (projections) from a number of climate models (RCMs or GCMs). However, there is currently an ongoing discussion on whether and to what degree nonstationarity should be addressed in design for hydroclimatic extremes, mainly focusing on uncertainty, model complexity, and data limitation issues [39,40].

Considering that there is an increasing evidence of climate change associated with destructive extreme events of higher intensity and frequency, combined with the deterioration process of the sewer network and the rapid growth of cities, an increase is expected in the number of people and properties in urban areas affected by the harmful effects of urban stormwater. In this context, urban drainage infrastructures are unlikely to be adequate in the future in the majority of cities worldwide, because their design relied upon the historical precipitation data ignoring the climatic variability. Kourtis and Tsihrintzis [41] reviewed the state-of-the-art of scientific approaches for adaptation of urban drainage networks to climate change, discussing climate change impacts on precipitation, IDF curves, flooding and urban drainage, and defined a novel approach for assessing urban drainage network adaptation to climate change and other drivers. Galiatsatou and Iliadis [29] introduced a methodological framework for designing urban drainage networks under climate change conditions, considering nonstationarity of the extreme precipitation in the study area. Galiatsatou et al. [42] conducted capacity assessment of a combined sewer network under the impact of rainfall storm events of different return periods, and proposed mitigation strategies, such as the application of nature-based solutions (NBSs) or low-impact developments (LIDs) for controlling combined sewer overflows. In this framework, the effects of climate change on the performance and capacity of the combined sewer system in Thessaloniki’s historic center, Greece.

2. Materials and Methods

2.1. Study Area

This study focuses on the historic center of Thessaloniki, a city in northern Greece. Thessaloniki is the country’s second-largest city, with a population exceeding one million and a particularly high population density in the city center, estimated at 200 residents per hectare. The city’s annual average rainfall was recorded at 445.7 mm for the period 1960–2020, based on monthly data from the Hellenic National Meteorological Service (HNMS). During this period, the average monthly rainfall varied from a low of 21.1 mm in August to a peak of 54.5 mm in December. The study area is defined by Aristotelous Street to the west and Aggelaki Street to the east, with its northern boundary along Egnatia and Svolou streets, and its southern boundary along Leoforos Nikis on the city’s coastal front. The selected area, covering approximately 40 hectares, along with the layout of Thessaloniki’s examined sewer network, is presented in Figure 1.

2.2. Climate Data Assessment and Downscaling Methodology

Historical precipitation measurements were obtained from the Meteorological Station (MS) of the Aristotle University of Thessaloniki (AUTH) which is located near the study area. These historical records covered the period from October 1965 to September 2005. The AUTH station has been widely recognized for its long-term, high-quality data collection, making it a reliable source for climatological analysis and model evaluation.

Future precipitation projections used in this study were sourced from the EURO-CORDEX Coordinated Downscaling Experiment [43], an initiative under the CORDEX framework. This dataset provides high-resolution regional climate simulations, ensuring the robustness of climate projections across Europe. For this analysis, monthly precipitation projections were derived from three distinct Regional Climate Models (RCMs) as presented in Table 1.

Each RCM provides a spatial resolution of 0.11o (approximately 12.5 km), offering a fine spatial grid suitable for local-scale climate impact studies. The future projections cover the period from October 2020 to September 2100 and were performed under the Representative Concentration Pathway (RCP) 4.5 scenario. This scenario represents an intermediate greenhouse gas emission pathway, where emissions peak around mid-century and then stabilize.

To ensure consistency and comparability, historical precipitation observations (1965–2005) were directly compared to the corresponding historical simulation outputs from the RCMs for the same period. Given the inherent biases present in raw precipitation outputs from RCMs, bias correction was necessary. The Bias Correction Quantile Mapping (BCQM) method was employed for this purpose, utilizing the qmap package [22] in R. This statistical technique adjusts the cumulative distribution of modeled precipitation to align with that of observed precipitation, thereby addressing systematic discrepancies. The BCQM approach uses a non-linear transformation to model the quantile–quantile relationship between observed and simulated precipitation values. Specifically, the exponential tendency to an asymptote transformation function was applied, as described by Gudmundsson et. al. [22]:

$${\widehat{X}}_o=(a+b\times X_m)(1-e^{-{\displaystyle \frac{\left(X_m-x\right)}{\tau }}})$$

where ${\widehat{X}}_o$ is the observed variable, X_m is the modelled variable, and a,b,x and τ are parameters.

The bias-corrected transformation function was subsequently applied to the future climatic projections to adjust for biases in the raw outputs of the climate models. To enable a more detailed analysis of future climatic trends, the projected period spanning from 2010 to 2100 was divided into two distinct 40-year intervals: 2020-2060, representing the short-term projection, and 2060-2100, representing the long-term projection. The mean annual precipitation and its associated standard deviation were computed for both the historical period and the future time intervals for each of the three climate models after bias correction.

2.3. The Scaling GEV Model

The GEV distribution has been used in several scientific studies for modelling annual series of extreme rainfall, as it has been proven that it matches well the observed extreme values. The cumulative distribution function (CDF), G(x), for the GEV distribution is given as:

$$G\left(x;\mu ,\sigma ,\xi \right)=\left\{\begin{array}{c}\exp \left\{-{\left[1+\xi {\displaystyle \frac{x-\mu }{\sigma }}\right]}^{-{\displaystyle \raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$\xi $}\right.}}\right\}\xi \ne 0\\ \exp \left[-\exp \left({\displaystyle \frac{x-\mu }{\sigma }}\right)\right]\xi =0\end{array}\right.$$

(2)

where μ, σ and ξ are respectively the location, scale and shape parameters, the latter also determining its limiting behaviour. When the shape parameter, ξ, is equal to zero the GEV distribution is equivalent to the Gumbel distribution, while when ξ is higher or lower than zero, the GEV corresponds to the Fréchet and the (reversed) Weibull distribution families, respectively. In this work, due to the small sample size of the annual maximum rainfall series (40-year maxima for all study periods and climate models) the three parameters are estimated by the method of L-moments [44,45]. For the random variable X with distribution function F, the theoretical L-moments λr+1 with r=0, 1,.. are expressed as linear functions of the specific probability weighted moments (PWM):

$$a_r=M_{\mathrm{1,0},r}=\mathrm{E}\left\{X{\left[1-F\left(X\right)\right]}^r\right\}$$

(3)

with the dimensionless L-moment ratios L-skewness, τ₃=λ₃/λ₂, and L-kurtosis, τ₄=λ₄/λ₂, calculated as functions of the L-scale, λ₂, and the third, λ₃, and fourth, λ₄, L-moments. Let a_r be the unbiased estimator of a_r for an ordered sample x_1:n≤…≤x_n:n [45]:

$${\mathrm{a}}_{\mathrm{r}}={\displaystyle \frac{1}{n}}{\sum }_{i=1}^n\left(\begin{array}{c}n-i\\ r\end{array}\right)x_{i:n}{\left(\begin{array}{c}n-1\\ r\end{array}\right)}^{-1},r=0,1,\dots ,n-1$$

(4)

with the unbiased sample L-skewness, t₃=l₃/l₂, and L-kurtosis, t₄=l₄/l₂, calculated as functions of the sample L-scale, l₂, and the third, l₃, and fourth, l₄, sample L-moments, respectively. If the L-moments estimation method results in a quite low estimate of the shape parameter (estimate close to -0.5) of the extreme precipitation process for a certain time interval and climate realization scenario, the maximum likelihood estimation (MLE) procedure is also applied to estimate the GEV parameters [46,47].

After the parameters’ estimation, the quantiles XΤ for each return period T can be calculated as follows:

$$X_T=\left\{\begin{array}{c}\mathsf{\mu }-{\displaystyle \frac{\sigma }{\xi }}\left\{1-{\left[-\mathrm{l}\mathrm{n}\left(1-{\displaystyle \frac{1}{T}}\right)\right]}^{-\xi }\right\}\xi \ne 0\\ \mathsf{\mu }-\mathsf{\sigma }\ln \left[-\mathrm{l}\mathrm{n}\left(1-{\displaystyle \frac{1}{T}}\right)\right]\xi =0\end{array}\right.$$

(5)

It can be shown that the k-th order of Non-Central Moments (NCM), m_k, can be expressed as a function of the three GEV parameters as:

$$m_k={\sum }_{i=0}^k\left(\begin{array}{c}k\\ i\end{array}\right){\left(-1\right)}^i{\left(-{\displaystyle \frac{\sigma }{\xi }}\right)}^{k-i}\Gamma \left(1-i\xi \right)$$

(6)

For k=1, 2, 3 equation (5) can be written as follows:

$$m_k={\left(\mu -{\displaystyle \frac{\sigma }{\xi }}\right)}^k+{\left(-1\right)}^k{\left(-{\displaystyle \frac{\sigma }{\xi }}\right)}^k\Gamma \left(1-k\xi \right)+k{\sum }_{i=1}^{k-1}{\left(-1\right)}^i{\left(-{\displaystyle \frac{\sigma }{\xi }}\right)}^{\iota }{\left(\mu -{\displaystyle \frac{\sigma }{\xi }}\right)}^{k-i}\Gamma \left(1-k\xi \right)$$

(7)

where Γ(.) is the gamma function. In this work, the temporal downscaling method proposed by Nguyen et al. [25] and based on the “scale-invariance” concept, is utilized to statistically downscale daily rainfall hindcasts and predictions into finer temporal scales to assist the construction of IDF/DDF curves under different climate realizations. By definition, a function is scaling (or scale-invariant) if f(x) is proportional to the scaled function f(λx) for all positive values of the scale factor λ. Then, there exists a function C(λ) such that:

$$f\left(\lambda x\right)=C\left(\lambda \right)f\left(x\right)$$

Nguyen et al. [25] proved that:

$$C\left(\lambda \right)={\lambda }^{\beta }\mathrm{a}\mathrm{n}\mathrm{d}f\left(x\right)=x^{\beta }f\left(1\right)$$

(8)

in which β is the scaling parameter and is considered constant. Therefore, the relationship between the NCM of order k, m_k, and the variable x, can be expressed as [48]:

$$m_k=E\left\{f^k\left(x\right)\right\}=a\left(k\right)x^{\beta \left(k\right)}$$

(9)

in which:

$$a\left(k\right)=E\left\{f^k\left(1\right)\right\}\mathrm{a}\mathrm{n}\mathrm{d}\beta \left(k\right)=\beta k$$

(10)

If equation (9) is log-transformed the NCM of order k can be written as:

$$Ln\left(m_k\left(t\right)\right)=Ln\left(E\left\{f^k\left(1\right)\right\}\right)+\beta kLn\left(t\right)$$

(11)

Therefore, the three GEV parameters for subdaily rainfall annual maxima can be assessed from the respective annual maximum daily GEV parameters. More specifically, the GEV parameters for the two different time scales t and λt (λ≤1) are assessed as follows: the shape parameters of durations t and λt are assumed equal, and the location and scale parameters of durations t and λt are related by λ^β. In the present work, a more elaborate approach is used. The shape parameter of the process, ξ, is assumed to change with duration and not to remain constant, in order to achieve better accuracy of the method. Precipitation can present different scaling regimes, therefore during distinct time intervals, the process can present different scaling parameters, β. Therefore, if a process is characterized by scaling regimes in the intervals [t₁, t₂] and [t₂, t₃], then it can be described as:

$$f\left(t\right)=\left\{\begin{array}{c}{\left({\displaystyle \frac{t}{t_2}}\right)}^{{\beta }_1}f\left(t_2\right),t\in \left[t_1,t_2\right]\\ {\left({\displaystyle \frac{t}{t_3}}\right)}^{{\beta }_2}f\left(t_3\right),t\in \left[t_2,t_3\right]\end{array}\right.$$

(12)

After assessing the GEV parameters of finer temporal scales from daily rainfall, IDF and DDF curves can be fully defined for all time periods and climate realizations considered. The general IDF and DDF relationship used in this work, considering that dependence in rainfall duration, d, and return period, T, can be modelled separately is:

$$i_{{\rm T}}\left(t\right)={\displaystyle \frac{K{T}^c}{d^{1-b}}}\mathrm{a}\mathrm{n}\mathrm{d}{d}_{{\rm T}}\left(t\right)=K{T}^c{d}^b$$

(13)

where K, b, and c are coefficients to be determined by the available data for each period of interest.

2.4. Hydrologic and Hydraulic Model Configuration

In this study, the InfoWorks ICM (Integrated Catchment Modelling) software was selected to evaluate and assess the drainage system within the study area [49]. Developed by Innovyze, InfoWorks ICM is an advanced tool designed for managing combined sewer and stormwater systems. It facilitates detailed hydrologic and hydraulic modeling, enabling the analysis of sewer networks, assessment of existing infrastructure capacity, and planning of network upgrades. The software has been extensively applied in urban catchment modeling, as well as in studies related to urban flood hazards and flood risk assessment [50,51,52,53,54,55]. InfoWorks ICM supports integrated 1D and 2D hydrological and hydraulic simulations. To model surface flow, the nonlinear shallow water equations were employed, while the continuity and momentum equations were solved using a finite volume scheme with a Riemann solver.

Detailed and comprehensive data on the sewer system within the study area were obtained from the Geographic Information System (GIS) database of Thessaloniki’s Water Supply and Sewerage Company S.A. (EYATH S.A.). The collected data included analytical information on system links (conduits), such as their lengths (m), diameters (mm), and materials, as well as details regarding system nodes (e.g., inlets, manholes), including ground and invert elevations (m). Shapefiles containing all network elements were imported into InfoWorks ICM using the software’s open data import center. Additionally, IDF (Intensity-Duration-Frequency) curves were developed based on historical rainfall events and data from the Regional Climate Models: CCLM, RACMO, and REMO. Rainfall amounts corresponding to return periods of 2, 10, 50, and 100 years with a duration of 1 hour were initially determined. Subsequently, these total rainfall amounts were distributed over the duration of the storm event using the alternating block method.

3. Results

3.1. Projected Outcomes of Climate Models

The results, presented in Table 2, Table 3 and Table 4, illustrate the variability and potential shifts in precipitation dynamics at different temporal scales, highlighting inter-model differences and uncertainties inherent in future projections, and emphasizing the necessity of incorporating variability and uncertainty considerations into climate adaptation strategies.

The analysis of annual maximum daily precipitation reveals significant variability in its temporal distribution across the three climate models, as shown in Figure 2. These variations highlight the distinct responses of each climate model to the projected climatic conditions under the RCP4.5 scenario.

For the historical period (1965–2005), the average annual maximum daily precipitation was estimated at 41.04 mm/day, serving as the baseline for comparison. Projections for the short-term future period (2020–2060) indicate that the average annual maximum daily precipitation is expected to decrease to 38.71 mm/day (a reduction of -5.68%) according to the CCLM4-8-17.v1 - CLMcom.ICHEC-EC-EARTH model. In contrast, an increase to 50.61 mm/day (a rise of +23.32%) is projected by the RACMO22E.v1 - KNMI.ICHEC-EC-EARTH model, while the REMO2009.v1 - MPI-CSC.MPI-M-MPI-ESM-LR model forecasts a decline to 37.33 mm/day (a reduction of -9.04%).

In the long-term future period (2060–2100), the projected trends exhibit also notable differences among the models. The CCLM4-8-17.v1 - CLMcom.ICHEC-EC-EARTH model projects a substantial increase in the average annual maximum daily precipitation to 52.95 mm/day, representing a +29.02% change relative to the historical baseline. Similarly, the RACMO22E.v1 - KNMI.ICHEC-EC-EARTH model predicts an elevated value of 48.16 mm/day, corresponding to a +17.35% increase. Conversely, the REMO2009.v1 - MPI-CSC.MPI-M-MPI-ESM-LR model estimates a moderate increase to 39.02 mm/day, equating to a -4.92% change compared to the historical period.

Furthermore, as illustrated in Figure 2, the annual maximum daily precipitation consistently exhibits average values that are higher than their respective medians. This characteristic underscore the right-skewness of the precipitation dataset, a statistical trait that reflects the occurrence of infrequent but intense precipitation events that significantly influence the mean. Right-skewed distributions are a well-documented feature of precipitation datasets, particularly in regions subject to pronounced seasonal variability or periodic extreme weather phenomena. This skewness is indicative of a climate regime where extreme events dominate precipitation trends, and it has important implications for hydrological modeling and risk assessment.

Further analysis of the box plots presented in Figure 2 reveals notable differences between the historical and future periods. Specifically, the whiskers of the box plots for the future periods (2020–2060 and 2060–2100) concerning the CCLM4-8-17.v1 - CLMcom.ICHEC-EC-EARTH model are significantly longer than those observed for the historical baseline (1965–2005). This observation suggests a wider range of variability in annual maximum daily precipitation in the future, indicating the potential for more extreme values. The lengthened whiskers in the future projections highlight the increased likelihood of extremely high values, emphasizing the importance of accounting for climate extremes in future planning scenarios. This is particularly relevant for infrastructure design, and flood risk assessment, where planning for variability and uncertainty is critical.

3.2. Development of IDF Curves

The scaling behaviour of rainfall extremes in the study area is assessed using available annual maximum rainfall data at the Thessaloniki airport (station Mikra) provided by the HNMS covering a period of 25 years (1963-1987), corresponding to durations of 5min, 10min, 15min, 30min, 1h, 2h, 3h, 6h, 12h, and 24h. L-moments approach is initially used to assess the GEV parameters for all rainfall durations of the available fine scale dataset. Equation (6) is then used to assess the three NCMs for all rainfall durations and a logarithmic plot of the first three NCMs of rainfall at Thessaloniki station with rainfall duration is then constructed (Figure 3). Figure 3 also includes the power laws (Eq.9) resulting for all three NCMs.

The log-linearity of the GEV NCMs of annual maximum rainfall is illustrated in the plot, indicating the power law dependency (i.e., scaling) of the statistical moments with duration (Eqs. 9 and 11). The slopes in the plot are proportional to the scaling exponent β. Two different scaling regimes are evident from Figure 1, one in the interval 5min to 30min and one in 30min to 24h, indicating a difference in extreme rainfall dynamic at the finest temporal scales. This difference is confirmed from linearity of the scaling exponent β(k) with moment order, k, shown in Figure 4(a, b) for the intervals of [5,30] min and [0.5, 24] h.

After assessing the relationships of GEV NCMs with extreme rainfall duration, the three NCMs of the daily rainfall predictions are used combined with the scaling exponent β(k), which is kept constant for all time intervals and climate realizations, to assess the coefficients a(k) (see Eq. 9) of the annual maximum precipitation process. These coefficients can be used to determine the newly defined NCMs m_k for each subdaily rainfall duration (using Eq. 9). The three GEV NCMs for each rainfall duration are then utilized to assess the respective GEV parameters by means of Eq. (6). After assessing the GEV parameters for all rainfall durations and for all 40-year intervals considered, the extreme quantiles (return levels) corresponding to return periods of 2, 5, 10, 20, 50, 100, 200, 500 and 1000 years are estimated using Eq. (4). IDF and DDF curves are then created for all 40-year periods considered (historical and future periods from RCMs) based on Eq. (13). Figure 5 presents IDF and DDF curves in Thessaloniki, Greece for the historical/reference climate period (1965-2005). Figure 6 presents IDF curves for all climate models considered in this work (CCLM, RACMO, REMO) and for both future periods 2020-2060 and 2060-2100.

Table 5 presents IDF and DDF curve equations (Eq. 13) for the historical/reference climate period (1965-2005), as well as for the two future periods (2020-2060 and 2060-2100) for all three RCMs studied in this work.

Rainfall extremes for the future period 2020-2060 are assessed higher than those of the historical climate only for RACMO RCM. More specifically 100-year rainfall events of 1h duration are almost 39% higher than those of the historical climate. CCLM and REMO RCMs produce lower 100-year rainfall extremes (1 h duration) up to 17.3% and 20%, respectively. The future period 2060-2100 appears more energetic with respect to extreme precipitation events, compared to the future period 2020-2060. In the future period 2060-2100, 100-year rainfall events of 1 h duration are almost 150%, 7%, and 51.5% higher than those of the future period 2020-2060 for CCLM, RACMO and REMO RCM, respectively.

3.3. Hydraulic Results

The rainfall storm events used in this study were derived from the IDF curves developed in Section 3.2, which were constructed using advanced statistical methods and techniques. These IDF curves were specifically designed for Thessaloniki’s city center. Total rainfall amounts corresponding to return periods of 2, 10, 50, and 100 years with a duration of 1 hour were utilized. The model that was calibrated and validated in Galiatsatou et al. [42] is also employed in this study.

Figure 7 illustrates the hydraulic results for current conditions under the worst-case scenario of a 100-year return period storm event. Meanwhile, Figure 8, Figure 9 and Figure 10 show the results derived from rainfall distributions based on the CCLM, RACMO, and REMO climate models, respectively, for the periods 2020-2060 and 2060-2100. Flooding correlation is assessed by examining conduit surcharge conditions and the hydraulic grade line (HGL) at nodes relative to a predefined theoretical basement elevation of 2 meters below ground level. To simplify result visualization, the color-coded figures display only the HGL. Black dots represent manholes, while dark blue lines along the coastal front indicate CSO points in the network. Sewer segments are color-coded as follows: green when the HGL remains more than 2 meters below the surface, purple when basement flooding occurs (HGL is less than 2 meters from the surface but remains underground), and red when surface flow is present. For simplicity, results for events with a return period of less than 100 years are not presented, as the issue is more effectively highlighted under the worst-case scenario.

Table 6 provides a summary of the simulation results and the impact of future climate projections on combined sewer overflow volumes (CSO volumes).

As shown in Figure 7, Figure 8, Figure 9 and Figure 10 and Table 6, in the near-term period from 2020 to 2060, the projected CSO volumes vary significantly across the models. The CCLM (10,735 m³) and REMO (10,012 m³) models predict lower overflow volumes compared to the existing conditions (12,273 m³), suggesting that these models anticipate either a decrease in extreme rainfall intensity or a redistribution of rainfall events. However, the RACMO model (17,117 m³) shows an increase in CSO volumes, indicating a potential rise in extreme precipitation events under this scenario. For the long-term period from 2060 to 2100, all models predict an increase in CSO volumes compared to existing conditions, highlighting a worsening overflow problem in the future. Among them, CCLM (25,514 m³) projects the highest overflow volume, more than doubling the existing conditions. This suggests that under this scenario, intense rainfall events will become significantly more frequent or severe, putting additional stress on urban drainage systems. The RACMO model (18,605 m³) also indicates a continuous increase in CSO volumes, though at a lower magnitude than CCLM, while REMO (15,799 m³) suggests a moderate increase over time. These projections carry important implications for future flood management strategies. The substantial increase in overflow volumes predicted by the CCLM model indicates that the most robust adaptation measures may be necessary in this scenario. These could include significant infrastructure enhancements such as increasing storage capacity, upgrading conveyance systems, and implementing nature-based solutions like green infrastructure. The variations among the models highlight the uncertainty in climate projections, underscoring the need for adaptable and flexible flood management strategies capable of addressing a broad spectrum of potential future conditions.

4. Conclusions

This study emphasizes the substantial effects of climate change on combined sewer overflows, particularly within the Thessaloniki City Centre network. The analysis of current and future IDF curves, derived from climate models CCLM, RACMO, and REMO, indicates that the system will encounter significant challenges due to increasing rainfall intensity. Modeling with InfoWorks shows that the expected rise in rainfall intensity may exceed the system’s capacity to handle peak flows, potentially leading to flooding risks and higher overflow volumes. To mitigate these issues, strategies such as separating stormwater and sanitary flows, disconnecting downspouts, and implementing nature-based solutions should be considered. These measures would improve the system’s resilience, reduce flood risks, and support long-term urban sustainability in the face of climate change. Future research will focus on addressing capacity issues using offline storage and low-impact development techniques.