Urban Flood Modelling under Extreme Rainfall Conditions for Building-Level Flood Exposure Analysis

: The expansion of urban areas and the increasing frequency and magnitude of intense rainfall events are anticipated to contribute to the widespread escalation of urban ﬂood risk across the globe. To effectively mitigate future ﬂood risks, it is crucial to combine a comprehensive examination of intense rainfall events in urban areas with the utilization of detailed hydrodynamic models. This study combines extreme value analysis techniques applied to rainfall data ranging from sub-hourly to daily durations with a high-resolution ﬂood modelling analysis at the building level in the centre of Thessaloniki, Greece. A scaling procedure is employed to rainfall return levels assessed by applying the generalised extreme value (GEV) distribution to annual maximum ﬁne-temporal-scale data, and these scaling laws are then applied to more reliable daily rainfall return levels estimated by means of the generalised Pareto distribution (GPD), in order to develop storm proﬁles with durations of 1 h and 2 h. The advanced ﬂood model, CityCAT, is then used for the simulation of pluvial ﬂooding, providing reliable assessments of building-level exposure to ﬂooding hazards. The results of the analysis conducted provide insights into ﬂood depths and water ﬂowpaths in the city centre of Thessaloniki, identifying major ﬂowpaths along certain main streets resulting in localised ﬂooding, and identifying around 165 and 186 buildings highly exposed to inundation risk in the study area for 50-year storm events with durations of 1 h and 2 h, respectively. For the ﬁrst time in this study area, a detailed analysis of extreme rainfall events is combined with a high-resolution Digital Terrain Model (DTM), used as an input into the advanced and fully featured CityCAT hydrodynamic model, to assess critical ﬂowpaths and buildings at high ﬂood risk. The results of this study can aid in the planning and design of resilient solutions to combat urban ﬂash ﬂoods, as well as contribute to targeted ﬂood damage mitigation and ﬂood risk reduction.


Introduction
Urban surface water floods are amongst the most widely distributed natural hazards, endangering lives and causing damage to properties worldwide.The extent and severity of the damage is a product of both the intensity and duration of extreme rainfall events (variable in space and time), and their interaction with the complex flowpaths in a city on the ground surface and below it [1].Impermeable surfaces in urban areas, impeding infiltration and creating overland flow which exceeds the drainage capacity of the existing infrastructure, renders cities vulnerable to flash floods.Climate change and urbanisation are expected to increase urban flood risk, contributing to different components of the flooding system.Climate change, associated with global warming and an increase in the frequency and severity of extreme weather events, is anticipated to intensify flooding hazards.Urbanisation will contribute to the increase in flooding hazards, caused by a decrease in infiltration, baseflow and lag times and an increase in runoff volumes and peak discharges [2,3], but also to the increasing impacts of urban floods (increase in potential flood damages) caused by the growth of settlements and assets in flood-prone areas followed by a rise in property value in such areas [4,5].
Nowadays, urban floods affect both developing and developed countries, with the impacts of pluvial flooding being a major problem due to frequent natural disasters in urban areas [6].Flood damage mitigation in urban areas includes both structural and non-structural measures.Structural measures involve urban flood defences and flood control structures or designing or upgrading stormwater and drainage networks.Non-structural measures mainly focus on flood early warning systems and preventive actions [5,7].Urban flood modelling combined with exposure assessment of the buildings and population affected represents a principal non-structural measure to effectively manage urban flooding events and their adverse effects, as well as a prerequisite for disaster prevention and mitigation.
Flood modelling is a powerful tool in understanding the hydrodynamics of historic flood events, and in some cases, the construction of accurate IDF (intensity-durationfrequency)/DDF (depth-duration-frequency) curves can be used to predict future events that will cause damage to the urban fabric [8].There exist several studies combining flood inundation modelling with hydrological modelling and the unit hydrograph theory, mainly focusing on fluvial flooding in small or larger catchments [9][10][11][12][13][14], using DEMs (Digital Elevation Models) or DTMs (Digital Terrain Models) with a computational grid resolution ranging between 5 m and 100 m.However, it should be noted that the combination of an incorrect representation of urban features in the flood model such as buildings, bridges, infrastructure, etc., and the resampling of the DEM/DTM multiple times might cause large inconsistencies [1,15,16] and overestimation of the flooding hazard in areas with minor inundation issues and vice versa.Unlike studies modelling fluvial flood inundation in urban areas, studies focusing on the exposure assessment of urban areas to pluvial floods are rather limited, and started receiving significant interest quite recently.Zhu et al. [17] used the LISFLOOD-FP hydrodynamic model to simulate flooding in Lishui City, China, and employed the building-scale population distribution map to assess the population affected.Park et al. [5] evaluated flood risk for different building types, conducting vulnerability and exposure analysis in five regions of Ulsan City, South Korea.Their analysis resulted in a classification of each building type into five risk-related classes.Stefanidis et al. [18] presented a coupling of hydrological and hydraulic modelling on a national scale to produce flood hazard maps regarding flooding exposure in residential areas and infrastructure in Greece.Bertsch et al. [19] presented a sensitivity analysis and validation of a generic flood exposure analysis following a large storm event in Newcastle upon Tyne, UK, where more than 70% of the inundated buildings in the area were correctly identified during the storm event.
Pluvial flood risk assessment in urban areas, associated with estimating the hazard, exposure, and vulnerability components for the affected system, is therefore a major challenge for future societies.There is a great need to combine hydrological and hydrodynamic modelling to understand the impacts of urban floods, the water flowpaths in a city, and the urban features exposed to high flood risk.This study combines a detailed analysis and modelling of extreme rainfall events in the centre of Thessaloniki, Greece, with an advanced hydrodynamic model, CityCAT, to simulate pluvial flooding, significantly assisting a reliable assessment of exposure to flooding.The results of this study can aid in the planning and design of resilient solutions against urban flash floods, as well as contributing to targeted flood damage mitigation and flood risk reduction.
CityCAT has previously been applied in studies in the UK [1,[19][20][21][22] and the USA [23] (Environmental Justice of Urban Flood Risk and Green Infrastructure Solutions-Urban Systems Lab, Urban Flooding, Equity, and Green Infrastructure (arcgis.com,accessed on 10 July 2023)), where detailed and reliable spatial datasets were available, such as DTMs, building footprints, green spaces, and roads.This study also aims to demonstrate and assess the universal applicability of CityCAT, even in regions where comparable datasets may not be readily available, emphasising the value of integrating the model with extreme rainfall data to enhance flood resilience in urban areas.The practical implementation of the model in this study will assist local authorities and engineers in improving their future flood adaptation strategies.This study also marks the first published implementation of a flood exposure analysis calculator at the building level in a large Greek city, as opposed to conventional assessments limited to flood zoning.

Study Area and Available Datasets
The historic centre of Thessaloniki city in Greece was the study area of the present work, located in the northern part of Greece.Thessaloniki is part of the municipality of Central Macedonia, and it is the second largest city in Greece with a population of around 814,000 (Thessaloniki Population 2023: worldpopulationreview.com,accessed on 5 June 2023).The dense city centre facing the coastal front is characterised by historic buildings, residential properties, marketplaces, and a few green open spaces (see Figure 1).This part of Thessaloniki has suffered from severe storms and flash floods in the last decade, causing significant damage to roads, basements, local stores, etc.It should be noted that, during storm events, the roads are seen to become the main flowpaths for floodwater.

Extreme Rainfall Assessment
Extreme value theory includes two main approaches for identifying and modelling the extreme values of a random process, namely the block-maxima approach where the extremes follow a generalised extreme value (GEV) distribution, and the peaks-overthreshold (POT) approach that fits the extremes using a generalized Pareto distribution (GPD).In the former approach, the observation period is divided into nonoverlapping equal intervals (of length usually equal to one year) and block maxima are selected to be fitted according to a GEV distribution [24]: Two different datasets are available in the study area for the analysis of extreme rainfall events.The first dataset consists of daily rainfall data at AUTh (Aristotle University of Thessaloniki) station located in the centre of the city and covering 64 years (1958-2021) of measurements (no missing data are present), obtained from the database of the School of Geology, AUTh.The second dataset includes monthly maximum rainfall depths for rainfall durations of 5 min, 10 min, 15 min, 30 min, 1 h, 2 h, 6 h, 12 h and 24 h at Mikra station, located in the eastern part of the city.This dataset was made available by the Hellenic National Meteorological Service (HNMS) and covers a period of 25 years (1963-1987).It is therefore evident that the second dataset includes rainfall measurements of finer temporal scales than the first one, but it contains only monthly maximum values, and its length is significantly shorter than that of the daily rainfall series available at AUTh.It should also be noted that the second dataset contains missing values.To proceed with the extreme value analysis of all available datasets, annual maxima were first extracted for both the daily and the sub-daily series, and each dataset was tested for stationarity and trends [24].The datasets examined satisfy the hypothesis of stationarity, while no statistically significant trends were detected.

Extreme Rainfall Assessment
Extreme value theory includes two main approaches for identifying and modelling the extreme values of a random process, namely the block-maxima approach where the extremes follow a generalised extreme value (GEV) distribution, and the peaks-overthreshold (POT) approach that fits the extremes using a generalized Pareto distribution (GPD).In the former approach, the observation period is divided into nonoverlapping equal intervals (of length usually equal to one year) and block maxima are selected to be fitted according to a GEV distribution [24]: where µ, σ and ξ are the location, scale and shape parameters of the distribution, respectively.The POT approach employs two probability distribution functions: one for the intensity of exceedances over an appropriately defined threshold, typically a GPD, and another for the number of events per year (typically a Poisson distribution, or alternatively a constant number is used).The cumulative distribution function of the GPD is given by [24]: where σ and ξ are the scale and shape parameters of the GPD, respectively, while u is the defined threshold value.The scale parameter of the GPD, sometimes referred to as the modified scale parameter, is expressed as a function of the respective GEV scale parameter as: considering that the shape parameter, ξ, of the GPD is equal to that of the corresponding GEV.The GPD-Poisson, which employs the Poisson distribution to model the number of exceedances over the threshold value per year, is characterized by three parameters, the exceedance rate, λ, the scale, σ, and the shape, ξ, parameters.
Considering the threshold of the POT approach, a high value improves the validity of the asymptotic approximation of the GPD, but at the same time increases the variance of parameter estimates because of the reducing dimensions of the excess sample.In contrast, a very low threshold may increase bias from model misspecification [25,26].Finding a trade-off between these two issues is critical in fitting an extreme value distribution and producing reliable estimates of extremes.The parameters of the extreme value distributions were assessed using both maximum likelihood estimation (MLE) and the L-moments [27] approach.

POT Threshold Selection
When the POT approach is used to model rainfall extremes, an appropriate threshold should be selected to detect exceedances and define the extreme sample.Various threshold selection methods have been proposed in the literature, such as empirical methods, distance measure approaches, or diagnostic plots such as the mean residual life plot and GPD parameter estimate stability plots [28].However, these approaches have a significant level of subjectivity in the threshold selection process.This work used two methodologies, aiding in creating a more automatic threshold selection process [29].These threshold selection methods were proposed by Bader et al. [30] (henceforward referred to as threshold selection method (a)) and Silva Lomba and Fraga Alves [31] (henceforward referred to as threshold selection method (b)) facilitating a less ambiguous and more objective selection of daily extreme rainfall events.
Threshold selection method (a) [30] considers a set of candidate thresholds u 1 < . . .< u l , each having n i exceedances, i = 1, . .., l.Let H 0 (i) denote the null hypothesis that the distribution of n i exceedances above the threshold u i follows the GPD.Following the forward stop rule of G'Sell et al. [32], a rejection rule was constructed by returning a cutoff level k, such that H 1 to H k are rejected: where a is a prespecified significance level and p i I = 1, . .., l are the corresponding p-values of the l hypotheses.If there is no k ∈ [1, . .., l], there is no rejection of the null hypothesis.
The p-values in Equation ( 4) are assessed using the Anderson-Darling (AD) test for each candidate threshold, with the respective statistic assessed as: where z (i) = F y (i) θn is the probability integral transformation of the order statistics of the exceedances y (1) ≤ . . .≤ y (n) , based on the maximum likelihood estimator of θ, θn , under the null hypothesis H 0 .F denotes the cumulative distribution function of the GPD for each candidate threshold.Threshold selection method (b) [31] is based on L-moments.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): with the dimensionless L-moment ratios L-skewness, τ 3 = λ 3 /λ 2 , and L-kurtosis, τ 4 = λ 4 /λ 2 , calculated as functions of the L-scale, λ 2 , and the third, λ 3 , and fourth, λ 4 , L-moments, respectively.Let a r be the unbiased estimator of a r for an ordered sample x 1:n ≤ . . .≤ x n:n : with the unbiased sample L-skewness, t 3 = l 3 /l 2 , and L-kurtosis, t 4 = l 4 /l 2 , calculated as functions of the sample L-scale, l 2 , and the third, l 3 , and fourth, l 4 , sample L-moments, respectively.A set of candidate thresholds {u i } I i=1 are first defined with I = 10 or 20 sample quantiles.The minimum Euclidean distance is then defined between the sample L-skewness, τ 3,u i , and L-kurtosis, τ 4,u i , for each threshold value and the respective quantities of the theoretical GPD curve: The best candidate threshold is then defined as: characterized by L-moment statistics that fall closer to the respective values of the theoretical L-moment ratio curve.

Scaling Rainfall Extremes
Rainfall features of different temporal scales can be linked using scaling models, mainly based on the multifractal behaviour of rainfall.Temporal downscaling and temporal disaggregation methods are used to produce finer-temporal-scale rainfall data from coarser resolution observations.Temporal downscaling usually refers to the generation of data of a high temporal resolution by means of statistical techniques, most commonly stochastic models, calibrated using information on the statistics of data from lower resolution temporal scales.Temporal disaggregation indicates the generation of high-resolution temporal data based on coarser time scales, so that the former add up to the totals of the second scale.This can be performed by means of a temporal partitioning of low-temporal-resolution amounts using a recursive rule or by repeatedly adjusting stochastic models to the generated high-resolution data.Within the general framework of temporal disaggregation, different methodologies have been developed in the literature.Some quite simple techniques, based on assumptions regarding the association of specific characteristics of the probability distribution functions of rainfall amounts at different time scales, have been developed by Liu and Wang [33] and Chen et al. [34,35], among others.However, these techniques do not represent the basic statistics of the fine temporal scales of precipitation in a satisfactory way, nor the intermittency of precipitation events.Precipitation stochastic generators are also utilized for temporal disaggregation purposes of rainfall amounts.Methods based on point-process models for the temporal disaggregation of hydro-meteorological data are quite widespread, producing satisfactory results [36][37][38][39][40]. Lee et al. [41], Salas and Lee [42] and Lee and Jeong [43] introduced a nonparametric model for the temporal disaggregation of hydro-meteorological variables, which incorporates k-nearest-neighbour resampling and a genetic algorithm.Temporal disaggregation of hydro-meteorological data is also performed by means of machine learning techniques (i.e., Kumara et al. [44]).
The hypothesis of scale invariance [45] is usually applied to link rainfall intensities of different temporal scales.More specifically, the hypothesis of scale invariance states that annual maximum rainfall intensities, I d and I λd , corresponding to durations d and λd, respectively, can be related by the following equation [46][47][48]: where equality corresponds to the similarity of probability distributions.The coefficient λ is the ratio of scale invariance between the known duration D and the duration to be assessed, d, and β is the self-similarity index of the studied rainfall process.The qth moments of rainfall intensity are obtained from Equation (10) as follows [46]: where β(q) is the scale exponent of order q, estimated by log-transforming Equation ( 11): Therefore, the exponent β can be assessed as the slope of the linear relationship described by Equation (12).The abovementioned scaling behaviour can also be detected in quantiles of rainfall intensities corresponding to durations d and λd, considering that their cumulative distribution function (CDF) has a standardized form independent of the rainfall duration [48].In this work, a scaling procedure is applied to rainfall intensity quantiles corresponding to different durations, considering that their CDF has a standardized form independent of the rainfall duration.The scaling laws are assessed for all return periods for rainfall durations from 5 min to 30 min and from 30 min to 24 h, considering that rainfall dynamics change quite significantly in convective events.
Intensity-duration-frequency (IDF) and depth-duration-frequency (DDF) curves based on local precipitation measurements summarise the relationships between rainfall dynamics, namely rainfall intensity or depth, duration and frequency (return period), and are currently utilized for engineering design and management applications, such as flood risk protection structures and infrastructure or flood-mitigation projects.The IDF (DDF) curves are constructed for different return periods, representing the variation of rainfall intensity (depth) with duration.In this work, IDF and DDF curves were created following the methodology described in Galiatsatou and Iliadis [48].Theoretical probability distribution functions were fitted to annual maximum or POT rainfall intensities of particular durations ranging from shorter periods, e.g., 5 min, to daily events.When annual maximum rainfall intensities (or depths) were available, the GEV distribution (Equation ( 2)) was fitted to the samples of different durations, and rainfall return levels were assessed as: for a defined return period, T. When the GPD (Equation ( 2)) is fitted to rainfall POTs, the return levels are given by [24]: where n y is the number of observations per year and ζ u is the total exceedance rate of the threshold u.Confidence intervals for return levels estimated using both the GEV and GPD are assessed using the delta method [49].

Flood Exposure
An efficient flood exposure tool, developed by Bertsch et al. [19], was used to calculate the flood exposure likelihood for buildings.Figure 2 presents the schematic workflow of the flood exposure analysis tool.Each building was assessed for flood risk using simulated flood depths in a 3 m buffer zone around its perimeter, where the mean and the 90th percentile values were calculated, with the proposed value for the buffer zone (see Figure 3) depending on the resolution of the computational grid (2 m in this study).A simple classification scheme shown in Table 1 was used to categorise buildings at low, medium and high flood risk.3) depending on the resolution of the computational grid (2 m in this study).A simple classification scheme shown in Table 1 was used to categorise buildings at low, medium and high flood risk.

Modelling System and Model Set Up
Over the last decade, many studies have reviewed hydraulic and hydrodynam models which have been developed to simulate surface flows [6,[51][52][53], and one of t most advanced and fully featured is the City Catchment Analysis Tool-CityCAT (fo full description, see Glenis et al. [20]) developed at Newcastle University.CityCAT i unique hydrodynamic model able to simulate fully coupled surface and pipe netwo flows, and it can represent natural drainage systems and built-up areas, with the expl  Over the last decade, many studies have reviewed hydraulic and hydrodynamic models which have been developed to simulate surface flows [6,[51][52][53], and one of the most advanced and fully featured is the City Catchment Analysis Tool-CityCAT (for a full description, see Glenis et al. [20]) developed at Newcastle University.CityCAT is a unique hydrodynamic model able to simulate fully coupled surface and pipe network flows, and it can represent natural drainage systems and built-up areas, with the explicit representation of buildings, where the buildings' footprint is excluded from the computational grid, and different types of blue-green infrastructure (BGI) [20][21][22]54] (such as blue/green roofs, water butts, swales, etc.), thus enabling the assessment of different alleviation measures.The model outputs include maps and time series of water depth, flow velocity, and the volume in and out of manholes, gully drains, buildings, etc.The required inputs are: (a) a highresolution Digital Terrain Model (DTM) to unlock the full potential of the model, although depending on urban layouts, a lower resolution model may still be very functional; (b) the buildings' footprint; (c) green spaces to calculate the infiltration with the Green-Ampt method [55]; and (d) the IDF and DDF curves to generate storm profiles with the rainfallrunoff method [56], to be applied over the whole domain with a uniform assumption.
The flood domain studied here was modelled using CityCAT for two different design storm events, with a magnitude of 1 in 50 years and a duration of 1 h and 2 h.The buildings' footprint was extracted from the ONEGEO data (https://onegeo.co/data,accessed on 10 May 2023) and the permeable areas from OpenStreetMap (https://www.openstreetmap.org, accessed on 10 May 2023).The computational grid was constructed using the DTM provided by the Hellenic Cadastre (http://www.ktimatologio.gr/en,accessed on 10 May 2023) at a resolution of 2 m (each cell with an area of 4 m 2 ), so the total number of computational cells in the flow domain was 125,192, covering an area of 0.78 km 2 .The representation of the buildings in the model was performed following the 'Building Hole' approach where a non-flow boundary is generated around buildings to redistribute the rainfall to the nearest grid square (for a full description and performance relative to other methods, see Iliadis et al. [1]).Following the previously described approach, this study explored the likelihood of flood exposure for 1165 buildings.The roughness coefficient (Manning's n) was defined as 0.02 for impermeable areas and 0.035 for permeable areas.Due to the limitations in the Hellenic National regulations in urban flood modelling, and the intended design of the combined sewer system for storms with a return period of 10 years in the city centre of Thessaloniki (based on design cross-sections of existing combined sewers), a simple assumption was made in this work, that 20% of the rainfall enters the drainage system.In other countries, i.e., the UK, there is an instruction to flood modellers, when they do not combine the drainage system with the surface, to exclude specific rainfall from the model (e.g., 6 mm-15 mm).

Extreme Rainfall Assessment
Threshold selection for the daily rainfall data was performed using both well-known threshold selection techniques, such as the mean residual life plot (MRL) and parameter stability plots, while also accounting for site-specific characteristics of extreme rainfall, together with the two new threshold selection techniques presented in Section 2.2.1.The MRL plot of the daily rainfall sample at AUTh station is presented in Figure 4.The MRL plot identifies a range of possible threshold values in the interval 10 mm ≤ u ≤ 30 mm. Figure 5 presents the GPD parameter stability plots for the modified scale (scale parameter of the GPD), and shape parameter, ξ, for thresholds in [10,30] mm.Based on stability characteristics of the GPD stability plots, while also considering the uncertainty of the parameters represented using 95% confidence intervals shown as vertical lines for each threshold, a threshold between 14 mm ≤ u ≤ 27 mm is considered to be a good candidate.The MRL plot identifies a range of possible threshold values in the interval 10 mm ≤ u ≤ 30 mm. Figure 5 presents the GPD parameter stability plots for the modified scale (scale parameter of the GPD), and shape parameter, ξ, for thresholds in [10,30] mm.Based on stability characteristics of the GPD stability plots, while also considering the uncertainty of the parameters represented using 95% confidence intervals shown as vertical lines for each threshold, a threshold between 14 mm ≤ u ≤ 27 mm is considered to be a good candidate.
The threshold selection method (a) was applied for threshold values in the entire range, 10 mm ≤ u ≤ 30 mm.The significance level selected was set at 5%.The forward stop rule applied did not explicitly identify a cutoff level.However, it has been observed that the threshold u = 22 mm was the only one giving a p-value of the AD statistic lower than 5%.The threshold selection method (b) provided clearer results, indicating u = 22 mm as the threshold providing a local minimum to the Euclidean distance criterion.The threshold u = 28 mm provided the global minimum of the distance d u i in the studied interval.However, this threshold level ended up with only 134 POT samples, corresponding to just λ = 2.09 exceedances per year.Therefore, the threshold u = 22 mm was selected to perform the extreme value analysis of the daily rainfall data, corresponding to 221 POT samples, with around λ = 3.45 exceedances per year.
(scale parameter of the GPD), and shape parameter, ξ, for thresholds in [10,30] mm.Ba on stability characteristics of the GPD stability plots, while also considering the unc tainty of the parameters represented using 95% confidence intervals shown as vert lines for each threshold, a threshold between 14 mm ≤ u ≤ 27 mm is considered to b good candidate.The threshold selection method (a) was applied for threshold values in the en range, 10 mm ≤ u ≤ 30 mm.The significance level selected was set at 5%.The forward s rule applied did not explicitly identify a cutoff level.However, it has been observed t the threshold u = 22 mm was the only one giving a p-value of the AD statistic lower th 5%.The threshold selection method (b) provided clearer results, indicating u = 22 mm the threshold providing a local minimum to the Euclidean distance criterion.The thre old u = 28 mm provided the global minimum of the distance  in the studied interv However, this threshold level ended up with only 134 POT samples, corresponding to j Using a threshold u = 22 mm, the GPD was fitted to daily rainfall maxima (Equation ( 2)) using both MLE and the L-moments approach, with both approaches providing consistent results, with higher return level estimates assessed using MLE. Figure 6 presents rainfall return level estimates assessed using Equation ( 14), with the parameters of the GPD calculated using the MLE approach.The black line represents maximum likelihood rainfall return level estimates, while the blue lines represent the upper 97.5% and the lower 2.5% confidence limits (95% confidence interval).Round marks in the return level plot correspond to measured data from the available extreme rainfall sample.It should be noted that the most extreme 24-hourly rainfall measurement was 98 mm, and was observed in 1985 and 2014.The maximum likelihood return level estimate significantly underestimates this value, while the upper 97.5% confidence limit seems to better fit the most extreme part of the observed sample.More specifically, for a return period of 64 years, equal to the daily rainfall sample length, the maximum likelihood estimate of the rainfall return level is about 84 mm, and the respective 97.5% upper confidence limit is 100.5 mm.Based on this finding, the upper 97.5% confidence limit for daily rainfall return levels was used in the scaling methodology presented in Section 2.2.2 to extract rainfall return levels corresponding to finer temporal scales.
The GEV distribution (Equation ( 1)) is then fitted to annual maximum rainfall intensities for time periods of 5 min, 10 min, 15 min, 30 min, 1 h, 2 h, 6 h, 12 h and 24 h, available at Mikra station for the period 1963-1987 using L-moments (due to the small sample size of this dataset).Rainfall return levels for the different durations are assessed for return periods of 2, 5, 10, 20, 50, 100, 200 and 500 years using Equation (13).For each return period, plots of Log(i) and Log(λ) are created, and linear functions are then fitted, dividing the plots into two parts, the first one corresponding to rainfall durations from 5 min to 30 min, and the second one from 30 min to 24 h.To end up with two different rainfall duration groups, a number of trials were performed considering different duration groups, and finally we selected those providing the highest coefficients of determination, R 2 , for all return periods.Figure 7 presents the linear relationships between the log-transformed quantiles (log-transformed return levels) of rainfall intensity and log-transformed scale factors of different durations, for return periods of 5 years (left panel) and 50 years (right panel).The plots include the linear function equations for rainfall durations in the intervals of 5 min to 30 min and 30 min to 24 h, and the respective coefficients of determination.Table 2 presents estimates of the self-similarity index (estimates of −β) assessed for all return periods and for the two groups of rainfall duration.lower 2.5% confidence limits (95% confidence interval).Round marks in the return level plot correspond to measured data from the available extreme rainfall sample.It should be noted that the most extreme 24-hourly rainfall measurement was 98 mm, and was observed in 1985 and 2014.The maximum likelihood return level estimate significantly underestimates this value, while the upper 97.5% confidence limit seems to better fit the most extreme part of the observed sample.More specifically, for a return period of 64 years, equal to the daily rainfall sample length, the maximum likelihood estimate of the rainfall return level is about 84 mm, and the respective 97.5% upper confidence limit is 100.5 mm.Based on this finding, the upper 97.5% confidence limit for daily rainfall return levels was used in the scaling methodology presented in Section 2.2.2 to extract rainfall return levels corresponding to finer temporal scales.The GEV distribution (Equation ( 1)) is then fitted to annual maximum rainfall intensities for time periods of 5 min, 10 min, 15 min, 30 min, 1 h, 2 h, 6 h, 12 h and 24 h, available at Mikra station for the period 1963-1987 using L-moments (due to the small sample size of this dataset).Rainfall return levels for the different durations are assessed for return periods of 2, 5, 10, 20, 50, 100, 200 and 500 years using Equation (13).For each return period, plots of Log(i) and Log(λ) are created, and linear functions are then fitted, dividing the plots into two parts, the first one corresponding to rainfall durations from 5 min to 30 min, and the second one from 30 min to 24 h.To end up with two different rainfall duration groups, a number of trials were performed considering different duration groups, and finally we selected those providing the highest coefficients of determination, R 2 , for all return periods.Figure 7 presents the linear relationships between the log-transformed quantiles (log-transformed return levels) of rainfall intensity and log-transformed scale factors of different durations, for return periods of 5 years (left panel) and 50 years (right panel).The plots include the linear function equations for rainfall durations in the intervals of 5 min to 30 min and 30 min to 24 h, and the respective coefficients of determination.Table 2 presents estimates of the self-similarity index (estimates of −β) assessed for all turn periods and for the two groups of rainfall duration.The self-similarity indices presented in Table 2 are then used in Equation ( 12) to te  The self-similarity indices presented in Table 2 are then used in Equation ( 12) to temporally downscale daily rainfall return levels assessed from fitting the GPD to data from AUTh station (1958-2021).More specifically, daily rainfall return level estimates corresponding to the 97.5% upper confidence limit (see Figure 6) are used in the scaling process.Rainfall return level estimates extracted for durations of 5 min, 10 min, 15 min, 30 min, 1 h, 2 h, 6 h, and 12 h are used to construct IDF and DDF curves for the study site.The formulas extracted to describe the IDF and DDF curves are given below: 16.63T 0.2152 t 0.7116 and p(mm) = 16.63T 0.2152 t 0.2884 (15) where T is the return period (years) and t is the rainfall duration (h).

Modelled Flow Depth
A detailed analysis of the areas where maximum water depths are highlighted and analysis identifying the critical roads during heavy rains will be presented in this section.The rainfall depth for 50-year events, specifically for durations of 1 h and 2 h, is evaluated.The findings indicate that the assessed rainfall depth is approximately 70% higher than the quantiles derived from the DDF curves extracted from the shorter-duration dataset spanning 25 years (1963-1987).This high difference is attributed to: (i) using the upper 97.5% confidence limit to assess daily rainfall return levels of the longer time series (1958-2021), (ii) fitting a POT model to the 64-year daily series to assess extreme quantiles, and (iii) missing observations in the shorter series perhaps leading to an underestimation of the extreme sample.To simulate these two storm events, the CityCAT model is employed.It should be noted that when selecting the duration for modelling purposes, it is essential

Modelled Flow Depth
A detailed analysis of the areas where maximum water depths are highlighted and analysis identifying the critical roads during heavy rains will be presented in this section.The rainfall depth for 50-year events, specifically for durations of 1 h and 2 h, is evaluated.The findings indicate that the assessed rainfall depth is approximately 70% higher than the quantiles derived from the DDF curves extracted from the shorter-duration dataset spanning 25 years (1963-1987).This high difference is attributed to: (i) using the upper 97.5% confidence limit to assess daily rainfall return levels of the longer time series (1958-2021), (ii) fitting a POT model to the 64-year daily series to assess extreme quantiles, and (iii) missing observations in the shorter series perhaps leading to an underestimation of the extreme sample.To simulate these two storm events, the CityCAT model is employed.It should be noted that when selecting the duration for modelling purposes, it is essential to consider the critical duration that triggers the most significant flood response, taking into account factors such as time-to-peak and other relevant characteristics.In the case of a catchment area spanning only a few square kilometres, a duration of 1 or 2 h is often sufficient to adequately represent the hydrological processes and capture the flood dynamics effectively.These durations are typically suitable for encompassing the key rainfall patterns and associated runoff generation within the catchment, enabling accurate flood modelling and analysis.The application of the CityCAT model in simulating the two storm events (50-year events with durations of 1 h and 2 h, see Figure 9) provides valuable insights into flood depths and water flowpaths within the study area.Note that the simulated storm events here exhibit similarities to previously observed storms as reported by the Hellenic National Meteorological Service (HNMS).
Hydrology 2023, 10, x FOR PEER REVIEW 14 of 20 simulated storm events here exhibit similarities to previously observed storms as reported by the Hellenic National Meteorological Service (HNMS).The flood maps produced by the model outputs, depicted in Figure 10 (maximum flood depths), clearly identify a major water flowpath along Agias Sofias street (see Figure 1 to locate the street), where the darker blue illustrates water depths exceeding 30 cm.This indicates that the street is highly susceptible to flooding during intense rainfall events.Furthermore, the presence of small ponds in various parts of the catchment, attributed to the complex and dense topography of the area, highlights the potential for localised flooding.Identifying these ponding areas is crucial for understanding flood risk and implementing measures to minimise the impacts, such as sacrificial zones, the creation of retention ponds, the improvement of surface drainage in specific locations, or converting the impermeable pavements to permeable pavements.
The study area's locations and roads, discussed below, have experienced substantial water buildup during intense rainfall events in the past, as reported by local authorities and residents.However, additional efforts are required to compare and confirm these observed occurrences with the results obtained through modelling.The modelled water depths of this work were calculated to estimate the maximum levels on the following roads (see Figure 1 to locate the streets): (a) Palaion Patron Germanou and Pavlou Mela.This particular area demonstrates a significant propensity for water pooling, with estimated water depths exceeding 30 cm.(b) Notably, Proxenou Koromila experiences frequent ponding, with estimated flood depths ranging from 25 cm to 41 cm.This road is particularly susceptible to water accumulation during storm events, which can lead to hazardous conditions.(c) In certain parts of Mitropoleos street, water depths exceeding 25 cm have been estimated.This poses a risk of localised flooding which would result in traffic disruption.
The estimated water depth and the flow direction for the two storms with 1 h and 2 h durations can be seen in Figures 11-13, where we zoom in on these areas.Overall, the contribution of a detailed flood model, such as CityCAT, is crucial to developing a better understanding of the flood dynamics, quantifying water depths with high accuracy, and The flood maps produced by the model outputs, depicted in Figure 10 (maximum flood depths), clearly identify a major water flowpath along Agias Sofias street (see Figure 1 to locate the street), where the darker blue illustrates water depths exceeding 30 cm.This indicates that the street is highly susceptible to flooding during intense rainfall events.Furthermore, the presence of small ponds in various parts of the catchment, attributed to the complex and dense topography of the area, highlights the potential for localised flooding.Identifying these ponding areas is crucial for understanding flood risk and implementing measures to minimise the impacts, such as sacrificial zones, the creation of retention ponds, the improvement of surface drainage in specific locations, or converting the impermeable pavements to permeable pavements.
The study area's locations and roads, discussed below, have experienced substantial water buildup during intense rainfall events in the past, as reported by local authorities and residents.However, additional efforts are required to compare and confirm these observed occurrences with the results obtained through modelling.The modelled water depths of this work were calculated to estimate the maximum levels on the following roads (see Figure 1 to locate the streets): (a) Palaion Patron Germanou and Pavlou Mela.This particular area demonstrates a significant propensity for water pooling, with estimated water depths exceeding 30 cm.(b) Notably, Proxenou Koromila experiences frequent ponding, with estimated flood depths ranging from 25 cm to 41 cm.This road is particularly susceptible to water accumulation during storm events, which can lead to hazardous conditions.(c) In certain parts of Mitropoleos street, water depths exceeding 25 cm have been estimated.This poses a risk of localised flooding which would result in traffic disruption.
The estimated water depth and the flow direction for the two storms with 1 h and 2 h durations can be seen in Figures 11-13, where we zoom in on these areas.Overall, the contribution of a detailed flood model, such as CityCAT, is crucial to developing a better understanding of the flood dynamics, quantifying water depths with high accuracy, and locating areas at high flood risk to improve inundation resilience in dense cities.

Exposure Likelihood to Buildings
In order to identify the urban features exposed to flood risk, an innovative tool was used, as described in Section 2.3.The analysis of flood exposure to buildings in the study area provides valuable insights into the vulnerability of urban features to flood risk.Note that this area has faced inundation issues from extreme events in the past, for which no formal reports exist, but are well known by local people.
Table 3 provides the total number of inundated buildings per scenario in the study area.The number of buildings classified as being at high risk for the first storm event (1 h duration) is 165, and that for the second storm event (2 h duration) is 186.These values are nearly twice as high as for the buildings with medium flood exposure.Most of the high-risk buildings are located on the streets mentioned in Section 3.2, where the flood depth is more than 30 cm.Furthermore, in the studied area of the city centre, many buildings house businesses, particularly on their ground floor, often containing vulnerable assets, while there also exist numerous buildings of historical value.

Exposure Likelihood to Buildings
In order to identify the urban features exposed to flood risk, an innovative tool was used, as described in Section 2.3.The analysis of flood exposure to buildings in the study area provides valuable insights into the vulnerability of urban features to flood risk.Note that this area has faced inundation issues from extreme events in the past, for which no formal reports exist, but are well known by local people.
Table 3 provides the total number of inundated buildings per scenario in the study area.The number of buildings classified as being at high risk for the first storm event (1 h duration) is 165, and that for the second storm event (2 h duration) is 186.These values are nearly twice as high as for the buildings with medium flood exposure.Most of the high-risk buildings are located on the streets mentioned in Section 3.2, where the flood depth is more than 30 cm.Furthermore, in the studied area of the city centre, many buildings house businesses, particularly on their ground floor, often containing vulnerable assets, while there also exist numerous buildings of historical value.

Exposure Likelihood to Buildings
In order to identify the urban features exposed to flood risk, an innovative tool was used, as described in Section 2.3.The analysis of flood exposure to buildings in the study area provides valuable insights into the vulnerability of urban features to flood risk.Note that this area has faced inundation issues from extreme events in the past, for which no formal reports exist, but are well known by local people.
Table 3 provides the total number of inundated buildings per scenario in the study area.The number of buildings classified as being at high risk for the first storm event (1 h duration) is 165, and that for the second storm event (2 h duration) is 186.These values are nearly twice as high as for the buildings with medium flood exposure.Most of the high-risk buildings are located on the streets mentioned in Section 3.2, where the flood depth is more than 30 cm.Furthermore, in the studied area of the city centre, many buildings house businesses, particularly on their ground floor, often containing vulnerable assets, while there also exist numerous buildings of historical value.Figure 14 illustrates the flood depths and the resulting flood exposure of buildings during the two generated storm events (50-year storm events with durations of 1 h and 2 h).The use of colour-coded zones helps to categorize buildings based on the flood depths in the buffer zone.The red-coloured buildings indicate high-risk, where the flood depth exceeds the 30 cm threshold.These buildings are estimated to be more vulnerable to damage from flooding, and it is crucial to prioritize them for adaptation measures and enhance their resilience to future flooding.Buildings depicted in orange indicate a medium risk of flooding, where damage from flooding is still significant.Lastly, a grey colour highlights the buildings at low risk, with minimal flood depths and lower vulnerability to flooding.

Storm Scenarios
Medium High 50-year event with a duration of 1 h 90 165 50-year event with a duration of 2 h 99 186 Figure 14 illustrates the flood depths and the resulting flood exposure of buildings during the two generated storm events (50-year storm events with durations of 1 h and 2 h).The use of colour-coded zones helps to categorize buildings based on the flood depths in the buffer zone.The red-coloured buildings indicate high-risk, where the flood depth exceeds the 30 cm threshold.These buildings are estimated to be more vulnerable to damage from flooding, and it is crucial to prioritize them for adaptation measures and enhance their resilience to future flooding.Buildings depicted in orange indicate a medium risk of flooding, where damage from flooding is still significant.Lastly, a grey colour highlights the buildings at low risk, with minimal flood depths and lower vulnerability to flooding.It should be noted that further investigation is needed into the 30 cm threshold to categorise buildings according to their flood risk, in order to provide more accurate estimations to understand the risk and the vulnerability profile of the city's buildings.

Conclusions
This study combines a detailed contemporary analysis of extreme rainfall events in Thessaloniki, Greece, with an advanced hydrodynamic model to simulate pluvial flooding, assisting in the reliable assessment of building exposure to flooding risks.A dual scheme is employed to assess extreme rainfall: (i) extreme daily rainfall, resulting from a It should be noted that further investigation is needed into the 30 cm threshold to categorise buildings according to their flood risk, in order to provide more accurate estimations to understand the risk and the vulnerability profile of the city's buildings.

Conclusions
This study combines a detailed contemporary analysis of extreme rainfall events in Thessaloniki, Greece, with an advanced hydrodynamic model to simulate pluvial flooding, assisting in the reliable assessment of building exposure to flooding risks.A dual scheme is employed to assess extreme rainfall: (i) extreme daily rainfall, resulting from a long daily series, is analysed using a GPD.Two threshold detection methods are applied, to assist a less ambiguous selection of daily extreme rainfall events.(ii) Extreme rainfall of shorter annual maximum series ranging from sub-hourly to sub-daily durations is analysed using the GEV distribution.A scaling procedure is applied to rainfall return level estimates assessed from (ii), and the resulting scaling laws are applied to the more reliable daily rainfall return levels of (i), in order to finally derive storm profiles with durations of 1 h and 2 h.The resulting storm profiles are used to drive the hydrodynamic model CityCAT to simulate flooding, estimate the water depths, identify the critical water flowpaths and finally assess the total number of inundated buildings through a novel exposure analysis calculator per extreme rainfall scenario in the historic centre of Thessaloniki.Furthermore: 1.
Typical storm events have durations spanning 1 h to 2 h, so both durations have been used here to see how sensitive the damages are to storm duration.For storms of the same return period, a modest increase is found for the 2 h storm relative to the 1 h storm.

2.
The CityCAT model provides valuable insights into flood depths and water flowpaths, identifying a major water flowpath along Agias Sofias street, which is highly susceptible to flooding during intense rainfall events.The presence of small ponds in various parts of the studied catchment further highlights the potential for localised flooding.

3.
The estimated likelihood of flood exposure to buildings reveals the vulnerability of urban features to flood risk.Due to the previous flood events in the area, the number of buildings at high risk for both storm events underscores the importance of addressing flood impacts on the built environment.4.
The modelling system is suitable for assessing the performance of flood-resilience strategies such as retention ponds, surface drainage improvements, and permeable pavements.
This study showcases the unique capabilities of CityCAT in its application to a country like Greece, which faces challenges of limited data availability.By leveraging globally accessible datasets, a high-resolution Digital Terrain Model (DTM, provided by the Hellenic Cadastre), and a detailed analysis of extreme rainfall events, this model facilitates a better understanding of the dynamics of urban flooding.It is noteworthy that in Greece, flood exposure analysis is conducted here for the first time at the level of individual buildings, moving away from the conventional approach of assessing flood risk in predefined zones.The identification of critical flow paths and the assessment of buildings at high flood risk serve as key considerations for future work.This includes expanding the catchment area, adding the current sub-surface drainage system or developing new synthetic methods to represent the system, implementing the model, and validating against historical storm events.These efforts are aimed at making informed decisions to develop flood-resilience solutions that safeguard people, assets, and infrastructure from future flood events.

20 Figure 1 .
Figure 1.An overview of the study area in Thessaloniki, Greece: (1) the city centre and the computational domain (with red colour); (2) the urban features, where grey denotes the buildings, green denotes the permeable areas and yellow to brown shading indicates the surface elevation of the area.

Figure 1 .
Figure 1.An overview of the study area in Thessaloniki, Greece: (1) the city centre and the computational domain (with red colour); (2) the urban features, where grey denotes the buildings, green denotes the permeable areas and yellow to brown shading indicates the surface elevation of the area.

Figure 2 .
Figure 2. Schematic workflow of the flood exposure analysis tool for the classification of buildings according to the water depth in the buffer zone[50].

Figure 2 .
Figure 2. Schematic workflow of the flood exposure analysis tool for the classification of buildings according to the water depth in the buffer zone[50].

Figure 2 .
Figure 2. Schematic workflow of the flood exposure analysis tool for the classification of buildin according to the water depth in the buffer zone[50].

Figure 3 .
Figure 3. Example of the buffer zone to calculate inundation depth from grid squares [50].

Figure 3 .
Figure 3. Example of the buffer zone to calculate inundation depth from grid squares [50].

Figure 4 .
Figure 4. MRL plot of daily rainfall in Thessaloniki for the interval 1958-2021.

Figure 4 .
Figure 4. MRL plot of daily rainfall in Thessaloniki for the interval 1958-2021.

Figure 5 .
Figure 5. GPD parameter stability plots for daily rainfall in Thessaloniki for the interv 1958-2021.

Figure 5 .
Figure 5. GPD parameter stability plots for daily rainfall in Thessaloniki for the interval 1958-2021.

Figure 6 .
Figure 6.Rainfall ML return levels and 95% confidence interval (mm) assessed by fitting the GPD to daily rainfall in Thessaloniki for the period 1958-2021.

Figure 6 .
Figure 6.Rainfall ML return levels and 95% confidence interval (mm) assessed by fitting the GPD to daily rainfall in Thessaloniki for the period 1958-2021.

Figure 7 .
Figure 7. Scaling of rainfall return level estimates at the station of Mikra, in Thessaloniki, for ret periods of 5 (left panel) and 50 (right panel) years.

Figure 7 .
Figure 7. Scaling of rainfall return estimates at the station of Mikra, in Thessaloniki, for return periods of 5 (left panel) and 50 (right panel) years.

Figure 9 .
Figure 9. Storm profiles corresponding to the constructed DDF curves for Thessaloniki, Greece: (a) 39 mm rainfall with a 1 h duration; and (b) 46 mm rainfall with a 2 h duration.

Figure 9 .
Figure 9. Storm profiles corresponding to the constructed DDF curves for Thessaloniki, Greece: (a) 39 mm rainfall with a 1 h duration; and (b) 46 mm rainfall with a 2 h duration.

Hydrology 2023 , 20 Figure 10 .
Figure 10.Example of maximum flood depths from a CityCAT simulation for a 50-year storm event with durations of (a) 1 h and (b) 2 h, for the centre of Thessaloniki.

Figure 11 .
Figure 11.Flood depths and flow direction (black arrows) for a 50-year storm event with durations of (a) 1 h and (b) 2 h at Palaion Patron Germanou and Pavlou Mela streets (marked as (b) in Figure 1).

Figure 10 . 20 Figure 10 .
Figure 10.Example of maximum flood depths from a CityCAT simulation for a 50-year storm event with durations of (a) 1 h and (b) 2 h, for the centre of Thessaloniki.

Figure 11 .
Figure 11.Flood depths and flow direction (black arrows) for a 50-year storm event with durations of (a) 1 h and (b) 2 h at Palaion Patron Germanou and Pavlou Mela streets (marked as (b) in Figure 1).

Figure 11 .
Figure 11.Flood depths and flow direction (black arrows) for a 50-year storm event with durations of (a) 1 h and (b) 2 h at Palaion Patron Germanou and Pavlou Mela streets (marked as (b) in Figure 1).

Figure 12 .
Figure 12.Flood depths and flow direction (black arrows) for a 50-year storm event with durations of (a) 1 h and (b) 2 h at Proxenou Koromila (marked as (c) in Figure 1).

Figure 13 .
Figure 13.Flood depths and flow direction (black arrows) for a 50-year storm event with durations of (a) 1 h and (b) 2 h at a part of Mitropoleos street (marked as (d) in Figure 1).

Figure 13 .
Figure 13.Flood depths and flow direction (black arrows) for a 50-year storm event with durations of (a) 1 h and (b) 2 h at a part of Mitropoleos street (marked as (d) in Figure 1).

Figure 13 .
Figure 13.Flood depths and flow direction (black arrows) for a 50-year storm event with durations of (a) 1 h and (b) 2 h at a part of Mitropoleos street (marked as (d) in Figure 1).

Figure 14 .
Figure 14.Maximum flood depths and flood exposure of buildings for a 50-year storm event with durations of (a) 1 h and (b) 2 h for the centre of Thessaloniki.

Figure 14 .
Figure 14.Maximum flood depths and flood exposure of buildings for a 50-year storm event with durations of (a) 1 h and (b) 2 h for the centre of Thessaloniki.

Table 1 .
Classification scheme to calculate flood exposure likelihood for buildings.

Table 1 .
Classification scheme to calculate flood exposure likelihood for buildings.

Table 3 .
Total number of inundated buildings per scenario for the centre of Thessaloniki.

Table 3 .
Total number of inundated buildings per scenario for the centre of Thessaloniki.