The two models—InVEST UFRM (ver. 3.13.0, (Natural Capital Project (2023))) and IMECOGIP WRRM (Alpha Release, Ruhr-Universität Bochum, Geogr. Institut (2021))—were selected for this study because they can be run with open access data and their processing time for each case/scenario is relatively low, provided that their results are applicable for different uses, including urban planning (Assumma et al. 2022). They have relatively similar modelling approaches using the SCS-CN method and similar outputs that can represent the FRS-runoff retention index (RRI) (Natural Capital Project 2023) and discharge coefficient (CW) (Ruhr-Universität Bochum, Geogr. Institut 2021). However, there are certain differences which render them appropriate for comparison, both in terms of the quality of the input data and the precision of the results for the flood regulation ES.

We developed a specific five-stage approach to assess the performance of the two models, based on different inputs, which enables the execution of models with different input data and comparative analyses of the results (Fig. 1). The first stage consisted of collecting and analysing all available and applicable data for the case study area, as well as studying the two models and their characteristics. In the second stage, the models were executed using the initial datasets. The next stage involved updating the initial datasets to improve their quality. The improvement in quality and accuracy was performed for the model input data—soil types, ecosystem types datasets and Urban Atlas classes by updating the already existing datasets and applying an approach for soil types update. Stage 4 included new modelling with already updated datasets and the final stage was data analysis and FRS assessment (Fig. 1). Each step is described in more detail in Subsections Input data, Flood regultion modelling and Comparison of model results.

Figure 1.

Conceptual scheme of the proposed approach.

UFRM—InVEST Urban Flood Risk Mitigation;

WRRM—IMECOGIP CICES 2.2.1.3: Waterflow Regulation (Retention Model);

CW—Discharge coefficient;

MAES BG—Mapping and Assessment of Ecosystems and their Services in Bulgaria.

Sofia is the capital city and largest urban area in Bulgaria. The case study area of Sofia is located in western Bulgaria and covers an area of 221.53 km² (Fig. 2). It represents urban ecosystems and their subtypes, as well as other ecosystem types, such as croplands, grasslands, woodlands and forests, heathlands and freshwater. The same case study has been studied as part of the INES project (Petkova et al. 2022).

Figure 2.

Case study area. A) Ecosystem subtypes (Level 3) (see Table S1, Suppl. material 1, Apostolova et al. (2017a), Apostolova et al. (2017b), Kostov et al. (2017), Velev et al. (2017), Yordanov et al. (2017), Zhiyanski et al. (2017)). B) Urban Atlas classes (2018) (see Table S2, Suppl. material 1, European Environment Agency 2020).

The boundaries of the area do not follow the administrative division of the Sofia Municipality, but reflect the spatial extent of the impact of the densely urbanised area. The case study was delineated following the extent of the urban heat island intensity in the 20–22 h interval (August 2019) (Sofiaplan 2020) and the boundaries were corrected according to the extent of the polygons of the ecosystem types from the national methodologies for Mapping and Assessment of Ecosystems and their Services in Bulgaria (MAES BG) (Maes et al. 2013, Apostolova et al. 2017a, Apostolova et al. 2017b, Karamfilov et al. 2017, Kostov et al. 2017, Sopotlieva et al. 2017, Uzunov et al. 2017, Velev et al. 2017, Yordanov et al. 2017, Zhiyanski et al. 2017).

Six of eight terrestrial ecosystem types from the MAES classification (Maes et al. 2013) were presented in the case study (Fig. 2A). These include urban, cropland, grassland, woodland and forest, heathland and shrub and freshwater ecosystems. Urban ecosystems were the predominant type presented by all ten subtypes developed for Bulgaria according to the methodology for urban ecosystems (Zhiyanski et al. 2017). The largest proportion was in J1: residential and public areas; J5: urban green areas (including sports and leisure facilities); J6: industrial sites (including commercial sites); and J7: transport networks and other constructed hard surface sites.

We also used classes from the Urban Atlas dataset (European Environment Agency 2020) to evaluate the performance of the models using different input data. There were 22 of 27 classes in the case study area, which corresponded to a highly diverse urban environment (Fig. 2B). Class 12100, ‘Industrial, commercial, public, military and private units’, was the most common as it covers over 1/5 of the case study area (21.74%). The other highly represented areas were classes 11210 ‘Discontinuous Dense Urban Fabric (S.L.: 50– 80%)’ (17.72%) and 14100 ‘Green urban areas’ (15.64%).

Sofia has most recently experienced flooding in 2005 (Council of Ministers of the Republic of Bulgaria 2022), July 2020 (Council of Ministers of the Republic of Bulgaria 2022) and June 2023 (National Institute of Meteorology and Hydrology of Bulgaria 2023). In the 2020 flood (around 50 mm/24 h in different locations), the most affected were the southwest districts and certain areas alongside the urban rivers Vladayska, Perloska and Slatinska (National Institute of Meteorology and Hydrology of Bulgaria 2020). The most recent flood, in June 2023 (61.30 mm on 15 June), had a slightly different spatial extent because the affected areas were in the south and southeast of Sofia. The three flooding events were preceded by periods of rainfall that exceeded the rainfall norms for each month.

The evaluation of model performance using different input data was the primary objective of this study. To perform the comparison, we have selected a few approaches to identify how the two models react to different input datasets and rainfall scenarios. We conducted a comparative analysis of FRS results in terms of two primary aspects. The assessment of FRS necessitates the identification of appropriate indicators to quantify the ecosystem regulation function. First, FRS maps derived from different models were compared and the results were analysed. The FRS values were normalised so they could be compared in each simulation. The final score for the FRS varied between zero (no retention capacity) and one (full retention capacity). Second, the FRS maps derived using different input data were compared and analysed. The analyses were based on the type and area coverage, as well as distribution of each LULC class and FRS capacity. FRS for each ecosystem type (subtype) and Urban Atlas class; the initial and updated datasets, were derived, as well as overall statistics (average, median, quartiles) for these datasets.

Another crucial aspect of input data analysis is modelling using different rainfall scenarios. We applied three rainfall scenarios (Table 2) using different combinations of LULC classes. After thoroughly reviewing and comparing the results of each simulation, we determined the performance of the models under different rainfall scenarios calculating the FRS for each LULC scenario and overall statistics (average, median, quartiles). For the analysis of different scenarios, a set of data was selected to better represent the differences in the model outputs. For this analysis, the simulation with initial and updated ecosystem subtypes, Urban Atlas datasets (for the LULC) and the 69 mm scenario were selected.

The geospatial data visualisation was performed using ArcMap, where each indicator provided on the assessment value had a different colour scheme. Data processing was performed using statistical operations within ArcMap and analysed using data processing software (MS Excel).

The InVEST UFRM outputs include three runoff retention parameters representing the amount of precipitation captured per grid cell (Table 4). Runoff values varied between 0 (no retention) and 69 mm (the entire amount of precipitation captured by the ecosystem). For the runoff retention parameters, these values are transformed into quantities (m³) and relative retention capacities ranging from zero to one. The latter was defined as the most appropriate to depict the spatial distribution of runoff retention (Fig. 3) and was used for further analyses of the FRS. These outputs provide the option for estimating the FRS of each ecosystem subtype and Urban Atlas classes.

Figure 3.

InVEST UFRM (A) and IMECOGIP WRRM (B) outputs that can be used for flood regulation supply mapping (based on the initial Urban Atlas dataset).

The only IMECOGIP WRRM output that has similarity to the UFRM output is the CW value, which serves as an avoided discharge (or avoided runoff) captured per unit (in our case, Urban Atlas classes). The CW varied from zero (fully avoided runoff) to one (no retention and full runoff) (Table 4). CW is the most appropriate for analysing the spatial distribution of runoff retention and it is a non-parametric coefficient which can be used for comparison with the RRI from InVEST UFRM.

The analyses showed that the RRI (InVEST UFRM) and CW values (IMECOGIP WRRM) were the most appropriate parameters to be used as indicators for quantifying FRS. Both parameters were measured in the same dimensions (zero–one) but in reverse, which resulted in inverse representations in the maps (Fig. 3). Therefore, the CW values were normalised for the FRS, whereas the RRI derived from the InVEST UFRM was used without changes. The runoff retention values vary from zero to one, where zero is full retention (full runoff, for example, roads) and one is avoided runoff (full infiltration, for example, forests).

The differences between the two models are due to the different approaches used to calculate the water flow components in them. The FRS from the UFRM model varies between zero and one based on the type of LULC datasets. Therefore, the results (FRS) for UFRM showed extensive variability within the LULC classes (Fig. 3 and see FRS using different data sources). In contrast, the FRS from the WRRM consisted of fixed values for each LULC class and changes in the class area did not reflect the final FRS values.

The outputs of the two models showed significant differences in the FRS values for each LULC class. The UFRM assesses different impervious urban classes with different values, as for the case study they vary between 0.1 to 0.5, whereas the WRRM assesses the FRS in these densely sealed areas with one fixed value (0.2) (where nearly all rainfall is transformed into runoff). The most significant difference was the lower FRS in the WRRM modelling results. The lower results (0.2–0.3) were derived from the lower assessment values assigned for impervious urban classes (primarily residential units).

Urban green areas received relatively similar values, varying between 0.8 and 1.0, where the UFRM values exceeded those from the WRRM. Neither model assessed green areas, urban green areas or urban forests, with large differences in flood regulation values. UFRM assessed green areas with similar values, varying between 0.9 and 1.0, whereas the WRRM provided slightly lower values (0.8–0.9) for the same areas.

In contrast, sub-urban areas within the case study received fluctuating values. The UFRM indicated relatively medium values (0.3–0.5), while the WRRM assessed the FRS as lower, between 0.3 and 0.4. Moreover, the variability of the flood regulation values was a notable difference between the two models. The WRRM provided lower assessment values to sub-urban areas, substantially lower than those in the same classes in the UFRM.

Artificial (impervious) areas, including industrial sites and transportation networks, also received varying FRS values in both models. UFRM assessed certain industrial sites with values varying between 0 and 0.3, whereas the WRRM values were closely clustered between 0.3 and 0.4, indicating a propensity for less rainfall to be transformed into direct runoff.

The two models have different results when used for a series of rainfall events for the same case study area. The UFRM shows changes in its output, especially the FRS values for different classes, whereas the WRRM did not adapt its FRS values for each rainfall scenario, which results in preserving the same values for all modelled events, regarding different precipitation depths. The FRS values of the two models were highly diverse (Fig. 9). First, their comparison showed a gradual increase in the FRS values for the UFRM results, whereas the WRRM results did not change in all modelled rainfall scenarios. This is a result of the input data, which included only LULC and CW values, whereas the use of soil datasets in UFRM can be interpreted as a second adjustment of the results.

Figure 9.

FRS based on different precipitation depth scenarios.

The FRS values derived from the UFRM had distinctive differences based on rainfall scenarios. Green urban areas (14100; Fig. 2) preserved or lowered (decreased by 7.5%) their FRS in all scenarios. The pattern of preserving the same FRS values was also observed in open waterbodies.

The changes in FRS are primarily within several urban classes, which correspond to a decrease in the efficiency of rainfall-runoff reduction. The most remarkable changes in FRS values can be found in the urban fabric-related classes (primarily 11100 and 11210) and industrial, commercial, public, military and private units (12100) which had decreased FRS by 18% (for the 20-years probability scenario) and by 40% (for the 1000-years probability scenario; from 0.44 to 0.26), based on the percentage changes.

The current study compared the performance of two ecosystem service assessment models, InVEST and IMECOGIP, while using two LULC classifications: ecosystem subtypes and Urban Atlas classes. The study added more insights to the current research as a comparative model assessment study was performed. The proposed soil data updating approach could be applied to different case studies so the used data matches the needs in the local context.

UFRM and WRRM are applicable to stormwater floods because they estimate the retention capacity of urban ecosystems to store part or the entire amount of incoming water. Both models have various outputs that can be used as indicators to quantify the FRS. In the UFRM, these are the runoff retention parameters that estimate the amount of captured water quantitatively (mm, m³ and RRI). The RRI is most appropriate for complex ES assessments when flood regulation is assessed together with other services (e.g. Cortinovis et al. (2022), Marino et al. (2023)). These quantitative parameters can be useful in flood risk assessment studies that need to calculate the amount of water captured to propose specific management measures. The flood risk service parameters of the UFRM (Table 3), in actuality, normalise the runoff retention values per spatial unit. This is useful for further analyses of the FRS per ecosystem subtype or LULC class. In the WRRM, the most appropriate output for the assessment of FRS is the CW, which is relatively similar to the RRI (UFRM).

The soil characteristics (including their initial water conditions) could affect the flood regulation potential of the ecosystems (Zepp et al. 2023). The UFRM requirement of soil data and their incorporation in the modelling process ensure more precise results of flood regulation assessment, in comparison to the WRRM.

After comparing results from the two models and different modelling scenarios, the InVEST UFRM using the Urban Atlas dataset (but connected with ecosystem subtypes) could be recommended for usage in FRS assessment in urban environments, as it is more comprehensive and precise, based on the input data and internal modelling processes. The model has been used in different environments and the results have been found to have different applicability (Cortinovis et al. 2022, Marino et al. 2023). However, the IMECOGIP WRRM could also be applied, at the initial stages of flood regulation ES assessment.

In this study, the models were tested to assess FRS, based on the available data. However, the analyses of the model outputs showed that the values built by the UFRM service can potentially be used to assess also flood regulation demand, but after including other sets of indicators representing the human activity concentration area (Dai et al. 2021) and socio-demographic indicators (Fang et al. 2023). However, they cannot fully cover the entire range of flood exposure factors. Our experience with the WRRM does not provide evidence for possible applications for the assessment of flood regulation demands. The least favourable scenario for runoff retention can be interpreted as a step towards generating a quantitative measure for assessing demand (as part of the supply scenarios).

The InVEST UFRM is sensitive to the presence of specific data types. For example, a biophysical table (.csv) must contain all LULC classes in the raster geospatial file; otherwise, the model would not complete the iteration. The same applies to the presence/absence of the soil type in the raster file. The modelling process has certain uncertainties related to the different stages of data processing. The spatial units, as well as the minimum mapping unit (for vector data) and resolution (for raster data), can be sources of uncertainty, as they both can cause uncertainty in the final results, especially when converting from one format to another. The other source of uncertainty for this input dataset is the LULC classification used. This can be an uncertainty in the interpretation of the final results, as different researchers use different classifications and cross-walking might be challenging. For the IMECOGIP WRRM, a predefined LULC was prepared and used, but its linkage (within the model literature) with other commonly used LULC classes was not fully described. This can lead to possible mismatches in the results. Another possible source of uncertainty is the lack of soil data in the WRRM input data. The soil data (infiltration) could be used as an adjustment point (one of the major influencing factors for flood regulation) of the other input datasets which are similar between both models.

With the growing number of models and tools for flood regulation assessment, the number of studies has increased; however, there is a lack of comparative work on model results with other similar models that can reveal their actual advantages and disadvantages (Veerkamp et al. 2023). Both models used in this study have strengths and limitations. The primary advantage of using the InVEST UFRM is that it includes the required data for the curve number (incorporated into the biophysical table) and soil data. This enables users to obtain comparatively more precise outputs. The soil hydrological groups are essential for the second verification (along with the curve number) and for improving the final results by connecting them to the local characteristics of the case study area.

A common limitation of both models is the lack of inclusion of the specific topography of the area (e.g. elevation, slopes), referring to the lack of relief input data in both models. The inclusion of such elements in the model can serve as a second round of verification of the model outputs (if we interpret the soils as the first). The approach suggested by Vallecillo et al. (2020) to adjust the curve number by slope is a good example of integrating topography into the modelling process and improving the results. We have tested applying the adjustment; therefore, the results were not prominent, as the urban territory of the case study does not exceed the 5° average slope. Regarding this, this approach could apply to areas with more than a 5° slope. The other suggested by Vallecillo et al. (2020) improvement that we applied includes the adjustment of the CN values by the percentage of imperviousness for each LULC class. This adjustment enables the differentiation of the artificial classes (sealed) and refining the CN number.

Another possible method of improving the input data and overall model performance could be to use the vector file format of the input and output data, which would contribute to preserving the original scale and accuracy of the input datasets rather than converting them to raster format (especially applicable for the UFRM model). A possible method to increase the performance of the UFRM model, including its usage with other models, is by preserving the attribute data from the initial AOI such that the output data can be used for comparison between different spatial scales.

The possibilities for improving the IMECOGIP WRRM can be categorised into two groups. The first includes all improvements in modelling data preparation. There are several mismatches between the LULC classes in the English and German versions of the sample data table, including the names and CW data. This could be solved by the possible inclusion of a clearer explanation of the urban classes, as there are two versions of them in the predefined table and sample data. The second group is related to possible output enhancements, which include suggestions for calculating the output layer of the ES balance output (diagram and vector file) within the spatial dimensions of the AOI (spatial units). Thus, the outputs (spatial units) should remain stable for comparison across scenarios following the changes in their time and spatial evolution. This can help in the following analysis, including the possibility of further applying box-plot analysis for the spatial units, in which the performance in the current version of the WRRM does not present differences in the LULC units.

The WRRM has the same FRS results when used for a series of rainfall events for the same case study area. The model did not adapt its FRS values for each rainfall scenario, which results in preserving the same values for all modelled events, regarding different precipitation depths. The reason behind this outcome is the lack of usage of soil dataset, resulting in fixed FRS values, without further adjustment of the LULC and CW values, which in the UFRM are included in the CN values (soil hydrological groups).

As both models assess the capacity of the ecosystem to provide specific flood regulation services due to extreme rainfall events and not river flood events, uncertainties regarding validation could be identified. The validation is done when there are enough field measurements, which in our study is not applicable, as the models use modelling variables which are not studied in detail and further research is needed for validation.