Modeling of flood regulation for ecosystem accounting: a case study of Ogosta river basin

The System of Environmental-Economic Accounting – Ecosystem Accounting (SEEA-EA) provides the methodological framework for organizing data about habitats and landscapes, measuring the ecosystem services (ES), tracking changes in ecosystem assets, and linking this information to economic and other human activity. It is proposed and developed by the United Nations Statistical Division (UNSD) through a series of guidance documents (UN et al. 2014; UN 2017; UN et al. 2021). SEEA-EA is a spatially-based, integrated statistical framework for organizing biophysical information about ecosystems, measuring ecosystem services, tracking changes in ecosystem extent and condition, valuing ecosystem services and assets, and linking this information to measures of economic and human activity (UN et al. 2021). It is considered a way to integrate biophysical and economic information on the changes in the stock of natural capital and the value of ES to provide regular information to decision-makers (Vardon et al. 2019). Water flow regulation ES and biophysical modeling are among the main topics in the individual ES part of the SEEA-EA framework and flood regulation ES is one of the important services. Journal of the Bulgarian Geographical Society, Volume 46 (2022) 3–10

The System of Environmental-Economic Accounting -Ecosystem Accounting (SEEA-EA) is a spatially-based, integrated statistical framework for organizing biophysical information about ecosystems, measuring ecosystem services (ES). Water flow regulation ES and biophysical modeling are among the main topics in the individual ES part of the SEEA-EA framework and flood regulation ES is one of the important services. Characterizing and assessing flood regulation is a challenging task as both assessment and accounts of this ES need various data which are usually not available through direct or indirect measurements, therefore modeling approaches of water regulation are much needed. Despite growing attention and studies using hydrologic models to assess and/or map flood regulation ES, the accounting of this service is still not well developed. In this paper, we present an approach for accounting flood regulation at a local scale using ArcSWAT modeling. It is based on the results of flood regulation ES assessment, where modeling results are used to quantify the ES indicators and delineate the service providing areas (SPA) and service demand areas (SDA). The actual flow of flood regulation is calculated as a ratio between ES demand and ES potential and it represents the area of SPA which corresponds to the demand for flood regulation represented by SDA. The results show that predominant flood regulations ES supply is provided by the forest ecosystem as well as the actual flow. The accounting of flood regulation is strongly determined by ecosystem extent mapping. The CORINE Land Cover (CLC) provides the most appropriate and available data for mapping ecosystem extent at smaller scales. However, at a larger scale, it is too coarse and the combination of Mapping and Assessment of Ecosystems and their Services (MAES) national ecosystem mapping gives better results.
Although, there is some progress in the accounting of waterrelated regulating ES and flood regulation in particular, further development in this area is needed (Vardon 2014;Vardon et al. 2019). Characterizing and assessing flood regulation is a challenging task as both assessment and accounts of this ES need various data which are usually not available through direct or indirect measurements, therefore modeling approaches of water regulation are much needed (Crossman et al. 2019). Modeling water flow and flood regulation is often data-intensive and also analytically complex and generally requires the use of hydrological models (UN 2017). Popular solution is the usage of specialized hydrologic models such as KINEROS (KINematic Runoff and EROSion model), SWAT (Soil and Water Assessment Tool), STREAM and hydraulic models such as HEC-RAS (Hydrologic Engineering Center) in combination with GIS based interfaces such as AGWA (Automated Geospatial Watershed Assessment) and ArcSWAT (Nikolov and Nedkov 2020). The model's application requires a large amount of data and relatively long time for adjustment, calibration and validation. However, they could generate the reliable parameter values for indicators such as soil infiltration, surface runoff or water yield, which could be used as indicators for quantification of the service (Nedkov and Burkhard 2012) and further for the needs of ecosystem accounting.
The hydrological model SWAT is the most frequently used tool for the assessment of water-related ES (Jujnovsky et al. 2012;Francesconi et al. 2016;Xu et al. 2017;Cheng et al. 2017;Lee et al. 2018;Netzer et al. 2019;Carvalho-Santos et al. 2019;Izydorczyk et al. 2019). The application effort for the SWAT model is particularly higher, which is due to the high time requirement, the training effort to apply, and the necessary post-processing (Lüke and Hack 2017). However, it is available for application in a GIS environment through the ArcSWAT tool, which enables it for application to a wide audience of ArcGIS users. It delivers detailed results on an HRU (Hydrological Response Units) basis, which enables reflection on the spatial variability of the studied service throughout the basin. The SWAT model has been successfully tested for application in Bulgaria in the area of the Ogosta river basin). Studies on flood regulation ES assessment (Nikolov and Nedkov 2020), water supply and water purification ES assessment (Boyanova 2015), integrated geodatabase for surface water, soil and groundwater pollution (Cherkezova et al. 2019), land cover changes (Stoyanova et al. 2020) ensure appropriate research basis for further development of SWAT modeling of the basin and generation of parameters which could be used in the ecosystem accounting.
Despite growing attention and studies using hydrologic models to assess and/or map flood regulation ES, the accounting of this service is still not well developed. The approach proposed by Vallecillo et al. (2020) provides some methodological clues, but it is applicable at the supranational (European Union, EU) level. Flood regulation ES accounting at a national and local level is still a research gap that needs to be filled. This paper aims to present an approach for accounting flood regulation at a local scale using ArcSWAT modeling. Klisurski manastir  The Neozoic and Mesozoic underlying rocks of the Ogosta River floodplain in the study areas are covered by alluvial deposits with a two-layer structure: a lower layer built of gravels and boulders with sand, and an upper layer of sandy-clayey deposits (Stoyanova et al. 2020). The climate is temperate-continental characterized by relatively warm summers and cold winters. The mean annual temperatures are between 10.1 and 11.1°С. The annual precipitation varies from 600 and 800 mm. The soil types in the valley are determined according to Food and Agriculture Organization (FAO) classification. The eutric and dystric Fluvisols, and Gleysols are the most widespread soil types within the bottom of the valley. Usually, gleyic and dystric Colluvisols, rendzic Leptosols, dystric Cambisols, and luvic Phaeozems are found in adjacent areas (Tcherkezova et al. 2019).

Accounting of flood regulation
Flood regulation services are the ecosystem contributions to the regulation of river flows by prevention and mitigation functions. They are derived from the ability of ecosystems to absorb and store water and hence mitigate the effects of the flood. The regulating role of wetlands, floodplains, and coastal ecosystems is usually emphasized but it is also important to pay attention to the functions of other ecosystems throughout river basins that control the processes of water balance. Hence, it is necessary to separate between a service production (supply) area and a service benefit (demand) area. The ecosystems affect the water balance mainly through the processes of interception and infiltration. The first one depends on the aboveground structure of the ecosystem (land cover) while the second is strongly determined by the soil properties. The surface runoff, which is the main factor for flood formation, depends also on abiotic factors such as rocks and topography. Regulating ecosystem services can have preventing or mitigating functions. In the first case, the ecosystems (i.e. forests) redirect or absorb parts of the incoming water (from rainfall), reducing in such a way the surface runoff and consequently the amount of rivers discharge. This ecosystem service plays its role before flood occurrence and in some cases, it can even prevent it. However, the flood mitigation function comes into effect when the flood is already formed. The ecosystems (i.e. flood plains and wetlands) provide retention space for the water surplus to spill, thus reducing the flood's destructive power (Nedkov and Burkhard 2012).
For the needs of ES accounting, it is necessary to quantify three different components that are essential to understand and properly assess the amount of service used, and therefore, the benefit generated. These key components are: 1) ES potential, which is the amount of ES that can be delivered sustainably; 2) ES demand which is the need for a specific ES by society; and 3) ES use or actual flow which is the amount of service that is mobilized (used) in a specific place and time (Vallecillo et al. 2020). ES potential is usually mapped based on the ecosystem's capacity to retain water. Mapping ES demand is based on the flood risk assessment and the outline of the areas under threat. The actual flow can be calculated as a ratio between ES demand and ES potential. The spatial relationship between ES supply and demand can be represented by the so-called Service Providing Areas (SPA) and the Service Demanding Areas (SDA) (Syrbe et al. 2017). For accounting purposes, the SPA corresponds to the service supply, while the SDA to the service use. The spatial interaction between SPA and SDA can be represented by the actual flow of flood regulation. For the needs of flood control accounting, Vallecillo et al. (2019) define the actual flow as the "extend of the demand with upstream protection from the upstream ecosystems with high runoff retention potential".
For this study, we adapted the approach proposed by  which is based on the results of flood regulation ES assessment. The main difference is the use of the ArcSWAT tool to quantify the ES indicators and delineate the SPA. Thus, the ES supply is assessed on the results of biophysical modeling and the results are in the form of a flood regulation supply capacity map presented in six categories ranging from 0 (no relevant capacity) to 5 (very high relevant capacity). The SPA is defined from the upper three categories of the assessment scale (medium, high, and very high capacity).

Hydrological modeling and flood regulation ES assessment
For this study, we utilized the ArcSWAT, which is a river basin/ watershed scale model developed to predict the impact of land management practices on water, sediment, or agricultural yields in large complex watersheds with varying soils, land, and management conditions over a long period. It requires specific information about the topography, land management practices, vegetation, soil properties, and weather in the watershed. The physical processes associated with water movement, sediment movement, crop growth, nutrient cycling, etc. are directly modeled by SWAT using input data. The watershed may be partitioned into several subbasins (Nikolov and Nedkov 2020). The input information for each subbasin is grouped or organized into the following categories: climate, hydrologic response units (GRUs), ponds/wetlands, groundwater, and the main channel draining the subbasin. HRUs are lumped land areas within the subbasin that are comprised of unique land cover, soil, and management combinations (Winchell et al. 2010).
The required spatial data for ArcSWAT includes a digital elevation model (DEM), and soil and land cover data. Free DEM data of 12 m resolution is obtained from ALOS PALSAR. The soil data are obtained from the Bulgarian Soil Resources agency. The land cover data are obtained from the Ministry of Agriculture and Food. More details about these data sources and their utilization for the model requirements are given in Nikolov and Nedkov (2020).
For the quantification of flood regulation ES supply, we selected three indicators: actual evapotranspiration, surface runoff, and water yield. The parameters of the indicators are generated through ArcSWAT simulations and obtained from the model result tables. Actual evapotranspiration in SWAT is calculated as evaporated rainfall intercepted by the plant canopy therefore higher amount indicates a lower capacity for flood regulation and vice versa. Surface runoff, or overland flow, is the flow that occurs along a sloping surface. A high amount of surface runoff indicates low capacity. Water yield is the net amount of water that leaves the subbasin and contributes to streamflow in the reach. A high amount indicated low capacity for flood regulation and vice versa.
The assessment of flood regulation supply is based on a 0 to 6 relative scale where: 0 is no capacity; 1 indicates low relevant capacity, 2 relevant capacity, 3 is medium capacity, 4 for high relevant capacity, and 5 for very high relevant capacity. Based on these assessments, estimation values representing the capacity of ecosystem services are calculated. For this purpose, we use the ArcSWAT result table, regarding the HRU. For each HRU there is information for several different indicators including the specified ones. The calculation of the assessment scores is made by a tool, which extracts data for the specified indicators from ArcSWAT results for a user-preferred date (Nikolov and Nedkov 2020).

Delineation of SPA and SDA
The delineation of SPA should be based on a calculation of the potential runoff retention provided by spatially explicit data to identify key areas for flood control (Vallecillo et al. 2020). In our case, the results of the flood regulation ES assessment ensure differentiation of the basin into areas with different capacities i.e., potential runoff retention. The areas with low or no capacity (categories 0 to 2) have limited retention potential and consequently are excluded from the service-providing areas. Therefore, the delineation of SPA is made by integration of the areas with higher capacity for flood regulation i.e., categories 3 to 5.
The SDA is a function of economic assets in floodplains and it is used to calculate the demand for flood control (Vallecillo et al. 2020). The delineation of service demand areas is based on the assumption that the most vulnerable areas would have the highest demand for flood regulation (Nedkov and Burkhard 2012). The vulnerability has different dimensions (e.g., social, economic, environmental, institutional) and the most vulnerable places are using different sources of demographic, statistical, topographic, and economic data (Nikolova et al. 2009). Such areas have the highest (5-value) demand for flood regulation. The analyses of their spatial extent showed that they are located within the floodplains and occupy predominantly urban land cover classes. The agriculture areas within the floodplains are also under direct hazard in cases of flood events so they are also marked by high relevant demand. Thus, the areas with high demand for flood regulation are delineated and form the SDA of the case study.

Actual flow and accounting tables
For the assessment of the actual flow of flood control by ecosystems, Vallecillo et al. (2020) propose an integration of the spatial dimension between the ES potential and the ES demand represented by SPA and SDA. The map of the actual ES flow of flood control is thus expressed as the number of hectares of demand (SDA) protected by upstream ecosystems (SPA) in a given year (Vallecillo et al. 2020). This approach enables quantifying the role of ecosystems to control floods in relative terms and the actual flow, defined in this way, is dependent on changes in ecosystems situated upstream. For this study, the extent of the SPA and SDA is defined using a 1x1 km grid. The grid is intersected with the SPA and SDA polygons and each cell of the grid is assigned to a particular category. The actual flow of flood regulation is calculated as a ratio between ES demand and ES potential and it represents the area of SPA which corresponds to the demand for flood regulation represented by SDA.
The accounting tables are structured to record the flow of goods and services among economic units, between economic units and the environment, and among ecosystem assets (UN et al. 2021). They represent predominantly the supply and use of ecosystem services. In the case of flood regulation, the core of ES accounts is focused on the amount of ES used (the actual flow), which refers to the transaction between ecosystems and socio-economic systems (Vallecillo et al. 2020). For this study, we calculated ES potential (based on SPA), ES demand (based on SDA), and ES actual flow. The ecosystem's extent is defined using MAES ecosystem classification and CORINE Land Cover (CLC) data. CLC data has been used in this study as the most appropriate data source for mapping ecosystems at a national level, and for consistency with the MAES mapping at the EU level. The CLC classes were correlated to the ecosystem subtypes to develop a relevance table (Hristova and Stoycheva 2021). The correlated classes were incorporated into the CLC GIS data and the resulting dataset is appropriate for mapping ecosystems. CLC data is available for 2000, 2006, 2012, and 2018, therefore the extent of ecosystem subtypes is defined for these four periods. The ecosystem subtypes are intersected with the data for SPA and SDA, which enables us to calculate the areas of SPA and SDA for each ecosystem subtype. The actual flow of flood regulation is calculated as a ratio between ES demand and ES potential and it represents the area of SPA which corresponds to the demand for flood regulation represented by SDA.

Hydrologic modeling and Flood regulation ES assessment
The ArcSWAT hydrologic model generated an HRU shapefile with 25 subbasins and 3007 unique HRUs. The subbasins have an estimated elevation between 352 m and 1260 m. The smallest subbasin (№ 7) has an area of 0.01 km 2 with just 3 HRUs. The largest subbasin (№ 6) has an area of 55.5 km 2 , which is 9.9% of the study area and contains 314 HRUs. The average area of the subbasins is 22.4 км 2 . About 332.9 km 2 (59.4%) of the case study is covered by HRUs with an area larger than 1 km 2 . Nevertheless, a significant part of the HRUs (29909) has an area smaller than 1 km 2 , and 2319 of them are even below 0.1 km 2 . The model was run to simulate the elements of the water balance for the period 2000 -2005. The calibration of the simulated results was performed using monthly runoff data from the hydrometric station of Gavril Genovo which is located at the outlet of the basin. The statistic parameters calculated to validate the results are coefficient of determination and Nash-Sitcliffe efficiency. The coefficient of determination is 0.61 which is satisfactory as the recommended values for successful calibration are between 0.5 and 1. The Nash-Sitcliffe efficiency is 0.52.
Flood regulation assessment necessitates the estimation of the indicators for a period of time with a flood event. The most significant flood event during the simulated period occurred on 07.08.2005. The increase in the runoff started on 05.08.2005 and the decrease to the normal level finished on 09.08.2005. Therefore, the period between 07.08.2005 and 09.08.2005 was chosen to quantify the flood regulation ES through the indicators of evapotranspiration, surface runoff, and water yield. The values of the selected indicators for this period were extracted from the simulated results using the tool, which extracts data for the specified indicators from ArcSWAT results for a user-preferred date. The values of the indicators for the selected period were normalized to according the assessment score scheme from 0 to 5. The results for the three indicators are given in Table 1. The overall score was calculated as an average of the three indicators.

SPA and SDA
The delineation of SPA was made using the results of the flood regulation ES assessment. The differentiation of the basin into areas with different capacities was used as a basis for the delineation. The resulting GIS layer with overall capacity was reclassified into two classes. The first covers the areas with low or no capacity (categories 0 to 2) which have limited retention potential and are not part of the service providing area. The second covers the areas with higher capacity for flood regulation (categories 3 to 5) and together they represent the SPA of the Ogosta river basin. The delineation of SDA was made using the results of the flood risk analyses. They cover the areas with high and very high demand for flood regulation (4 and 5-value).
The SDAs comprise 420,1 km 2 which is about 80% of the whole basin (Fig. 2). SPAs are located evenly throughout the basin. They have a compact extent and form a continuous area from the south to the north in the higher part of the basin. More heterogeneous is their distribution in the lower parts of the basin especially the downstream of the Dalgodelska Ogosta and the Chiprovska rivers. There are also SPAs in the floodplain especially downstream of the main river around the villages of Gorna Kovatchitsa and Belimel.
The SDAs comprise 6.3 km 2 which is 1.2% of the whole basin. There are 25 separate areas with an average size of 0.57 km 2 located predominantly downstream of the main river and the Dalgodelska Ogosta tributary.

Actual flow and accounting table
The extent account is performed using a 1x1 km grid. The grid is intersected with the SPA and SDA polygons and each cell of the grid is assigned to a particular category. Then, the grid is intersected with the CLC data (available for 2000, 2006, 2012, and 2018) and the results are distributed into ecosystem subtypes following the MAES typology and its implementation in the mapping of ecosystems in Bulgaria. The actual flow of flood regulation is calculated as a ratio between ES demand and ES potential and it represents the area of SPA which corresponds to the demand for flood regulation represented by SDA.
The accounting table contains the ES potential (calculated by SPA), ES demand (calculated as SDA), and ES actual flow for four periods corresponding to the time series of CLC data ( Table 2). The predominant part of the ES potential (76%) is provided by Woodland and forest ecosystems. Grassland (13%) and cropland (8%) have a small contribution, while the urban and sparsely vegetated areas have almost no contribution to the ES potential. The changes in woodland and forest ecosystems during the four periods are quite small. There  ES demand areas are located predominantly in urban (27%) and cropland (69%) ecosystems. They also do not change significantly over time, but it should be mentioned that both ecosystems show a slow decrease during the whole period.

Discussion
The results show that hydrologic modeling ensures appropriate data for quantification of water balance parameters which can be used for indicators of flood regulation ES. This confirms the conclusions made in previous studies on the modeling of flood regulation ES (Nedkov and Burkhard 2012;Stürck et al. 2014;Boyanova et al. 2016a). The specific contribution of this study is the application of the SWAT model which is the most popular hydrologic model. Its GIS applications (ArcSWAT and QSWAT) make it easy to work with spatial data and generate spatially explicit outputs appropriate for accounting purposes. A combination of model results with further data from hydrological measurements or flood risk assessment can add further opportunities to obtain more precise accounting-related data.
The results presented here seem to deliver a good representation of the prevention flood regulation function. Better estimations of the mitigation function could be achieved by incorporating a hydraulic model that ensures the calculation of the amount of water absorbed by floodplains and wetlands. This would give the opportunity for a more precise estimation of ecosystems that have predominantly prevention flood regulation function.
Besides the advantages of the approach described above, there are also shortcomings. One of the model's limitations is that it is not applicable to floods caused by snow melting. The application effort for the SWAT model is particularly higher, which is due to the high time requirement for the preprocessing of input data, and the training effort to apply the model (Lüke and Hack 2017). This makes the application of the approach limited to several river basins but not to the whole country. Therefore, in the accounting of flood regulation on a national scale, another approach is needed.
The accounting tables produced in this study are based on calculations of spatial units derived from land cover data correlated to the MAES ecosystem types. However, the flood regulation function is strongly dependent also on the soil properties and the topography (especially slope). The landscape studies which incorporate such characteristics could be an appropriate solution for this problem. One particular example is the ecosystem services assessment based on multi-level landscape classification (Prodanova, 2021). The integration of the traditional landscape approach and the ES assessment into ecosystem accounting is a promising perspective.

Conclusions
The proposed approach employs a combination of hydrologic modeling, GIS-based techniques, and ES assessment to produce flood regulation accounting in relatively small river basins. The application of the SWAT model in the ArcSWAT GIS environment en-sures automated and precise delineation of the SPA. The flood regulation ES demand analyses based on flood risk assessment data ensure appropriate information for the delineation of SDA. The calculation of actual flow and development of accounting tables is based on the concept provided by Vallecillo et al. (2020) adapted to the specifics of the local-scale study. The critical point of the approach is the identification of the SPA based on the hydrological modeling results . The main question that needs to be solved is how to determine the threshold value which outlines the SPA.
The accounting of flood regulation is strongly determined by ecosystem extent mapping. The CLC provides the most appropriate and available data for mapping ecosystem extent at smaller scales. However, at a larger scale, it is too coarse and the combination of MAES national ecosystem mapping gives better results. The correlation between CLC and the MAES ecosystem subtypes at the national level in Bulgaria (Hristowa and Stoycheva 2021) gives an appropriate basis for the elaboration of more precise accounting tables. The main challenge for the near future is to develop the approach for application and effective integration into accounting and reporting systems at the national level.

Funding program
This study was funded by the MAIA (Mapping and Assessment for Integrated ecosystem Accounting) project, that receives funding from the European Union's Horizon 2020 research and innovation programme under grant agreement No 817527 and INES (INtegrated assessment and mapping of water-related Ecosystem Services for nature-based solutions in river basin management) funded by the Bulgarian National Science Fund grant agreement КП-06-Н54/4.