Pre/Post-Fire Soil Erosion and Evaluation of Check-Dams Effectiveness in Mediterranean Suburban Catchments Based on Field Measurements and Modeling

: The present study was conducted in the suburban forest of Thessaloniki city (Seich Sou), which constitutes one of the most signiﬁcant suburban forests in Greece. In 1997, more than the half of the forest area was destroyed by a wildﬁre, after which soil erosion and ﬂood control works (check-dams) were constructed in the burned areas. The aim of the study is to estimate the annual soil erosion rate for the last 30 years (pre- and post-ﬁre periods) applying the WaTEM/SEDEM model, in order to investigate the impact of this wildﬁre on soil erosion, the effectiveness of the ﬂood- and erosion-control works and the level of forest regeneration. It is the ﬁrst time that WaTEM/SEDEM was calibrated and validated in Greece, taking into account soil erosion records from the 18 check-dams that were constructed in the study area in 2001. The mean annual erosion rate was 0.0419 t/ha/year, 0.998 t/ha/year and 0.08 t/ha/year for the pre-ﬁre period, the ﬁrst 3 years and 20 years after the ﬁre, respectively. The results showed a very low erosion rate for the pre-ﬁre period, an expected signiﬁcant increase 3 years after the wildﬁre and a gradual decrease in the subsequent years until 2021. However, it seems that the post-ﬁre regeneration of the forest has not been fully achieved, since the annual soil erosion rate at the long-term post-ﬁre period is double compared with the pre-ﬁre period. Concerning the check-dams’ effectiveness, it was observed that after 20 years of operation, they were non-silted, and most of them retained a small amount of sediments. This fact could be attributed to multiple factors such as the very thin soil depth, ﬁre severity and catchment geomorphology, though the main reason seems to be the time elapsed between ﬁre occurrence and the check-dams’ construction. The results of this study advance/strengthen the knowledge concerning the pre/post-ﬁre soil erosion processes in already degraded ecosystems, while the calibrated model could serve as a useful tool able to be applied in other Mediterranean catchments of similar characteristics.


Introduction
Soil erosion constitutes a serious environmental problem, appearing on a global scale and posing an environmental threat of different levels of intensity among different regions [1]. It involves two processes: the erosive power of rain, which detaches soil particles, and the subsequent transportation through the surface runoff and wind [2]. Soil erosion rates in the Mediterranean region are expected to accelerate because of the frequent appearance of intense rainfall episodes, which contribute to the intensification of the erosion process and, therefore, the watersheds' degradation [3,4]. Significant environmental consequences of soil erosion include water quality deterioration, soil degradation, floods and land disasters [5], as well as increased sediment generation, transportation and deposition in natural lakes and water reservoirs [6][7][8]. In the natural environment, the balance of soil loss is achieved by soil regeneration through rock weathering and pedogenesis. However, the rate of erosion may increase because of the destruction of the protective plant cover, potentially caused by biotic or abiotic factors. Changes in land uses and non-sustainable hydrological management of watersheds could substantially influence the sediment generation and transport, as well as the surface runoff. This fact could affect the geomorphology of the stream bed and potentially give rise to undesirable, catastrophic results with potential disastrous consequences [9,10]. The main factors that inevitably entail significant changes in land uses are the abandonment of farms and pastures, the expansion of touristic activities, urbanization, as well as deforestation (forest fires, unsustainable logging, insect outbreaks) [11][12][13][14][15][16].
Density, structure, type and diversifications of vegetation in combination with landuse changes form a critical complex of factors that influence the production and transportation of sediments and soil erosion rates [17][18][19]. The increase of vegetation density contributes to the reduction of the raindrop impact force and prevents the detachment of soil particles and the transport of sediments, resulting in the stabilization of sediment production slopes [10,20,21]. The introduction of meadow vegetation on slopes of former arable lands significantly reduces the runoff and sediment transport [22], while the same effect is observed after the reforestation of watersheds [10,23,24]. Additionally, decrease of connectivity between sediment production areas and the hydrographic network could potentially prevent the discharge of sediments out of the watersheds [25,26].
In Mediterranean watersheds, two well-established practices of watershed management are reforestation and the construction of check-dams in first-and second-order streams, which are applied with the aim of decreasing the soil loss rates [27][28][29][30]. It is widely accepted that the construction of gravity check-dams is an effective management practice for the controlling of erosion and flooding [31][32][33]. Check-dams reduce water velocity and water drag capacity [34][35][36][37], while contributing to the stabilization of the riverbed, as the accumulation of sediment materials behind the check-dams prevents landslides and soil subsidence [34,[38][39][40]. A series of check-dams could potentially achieve a compensation gradient, and subsequently an equilibrium gradient, to fix the longitudinal profile, to reduce channel gradients and slope failures or sidewall erosion [41][42][43][44][45]. In parallel, check-dams reduce the amount of sediment losses, retaining large amounts of ground materials [46][47][48].
It is widely known that soil loss equations and the respective modeling software could be efficiently used in the estimation of soil erosion rates [49][50][51]. Applying both GIS techniques and mathematical models such as universal soil loss equation (USLE) and revised universal soil loss equation (RUSLE) has been proven to be an effective method for the estimation of both the quantity and spatial distribution of soil erosion [52][53][54]. Widely applied erosion modeling applications are the Unit Stream Power-based Erosion Deposition (USPED), Integrated Valuation of Environmental Services and Tradeoffs (InVEST), Lisem Integrated Spatial Earth Modeller (LISEM), Spatially Distributed Scoring Model (SPADS), Soil and Water Assessment Tool (SWAT) and Water and Tillage Erosion Model and Sediment Delivery Model (WaTEM/SEDEM) [55,56]. In Greece, USLE and RUSLE have been mainly used so far in combination with GIS for modeling erosion processes [57][58][59][60], as well as the SWAT model [61][62][63][64][65].
The WaTEM/SEDEM is a spatially distributed model, which estimates the long-term annual erosion rate and sediment transport and deposition at the catchment scale [66,67]. The model has been calibrated and validated in the past using data with a 20 × 20 m spatial resolution [66][67][68]. Recently, the model has been used by the European Union to create a database of expected erosion rates in the context of the goals towards zero soil degradation [69,70]. WaTEM/SEDEM has been previously used in the Mediterranean region [10,[70][71][72][73][74], though without extended application in Greek territory [75,76]. The two previous studies that were conducted in Greece did not calibrate the model and used the default model parameters, a fact that constitutes a limited factor for the accuracy of model predictions.
The aim of the current study is to estimate the annual soil erosion rate for the last 30 years using the WaTEM/SEDEM model in two catchments of the suburban forest of Thessaloniki (Seich Sou), which suffered a forest wildfire in 1997. Seich Sou is one of the Land 2022, 11, 1705 3 of 18 most significant suburban forests of Greece, located at the northeast side of Thessaloniki city. After the forest fire, soil erosion and flood control works were constructed in the burned areas. The model was calibrated and validated employing sediment records derived from eighteen (18) check-dams, which have been constructed in the hydrographic network since 2001. According to the authors' knowledge, this is the first time that WaTEM/SEDEM was calibrated in Greek catchments under post-fire conditions. The main goal is to estimate the sediment production and transportation rates in the years subsequent to the wildfire and to investigate the potential effect of wildfire on vegetation regeneration and recovery. The calibrated model would undoubtedly serve as a tool to be applied in other Mediterranean watersheds of similar characteristics.

Study Area
The study area consists of two typical Mediterranean catchments located in the suburban forest of Thessaloniki (Seich Sou), the "Konstadinidi" and "Ano Toumba" catchments ( Figure 1, Table 1). The motivation for the selection of these two catchments was the presence and operation of many technical works (check-dams) in the area, the significant changes of vegetation cover caused by the wildfire spread on 7 July 1997 and the easy accessibility of the locations where the technical works have been constructed.
default model parameters, a fact that constitutes a limited factor for the accuracy of model predictions.
The aim of the current study is to estimate the annual soil erosion rate for the last 30 years using the WaTEM/SEDEM model in two catchments of the suburban forest of Thessaloniki (Seich Sou), which suffered a forest wildfire in 1997. Seich Sou is one of the most significant suburban forests of Greece, located at the northeast side of Thessaloniki city. After the forest fire, soil erosion and flood control works were constructed in the burned areas. The model was calibrated and validated employing sediment records derived from eighteen (18) check-dams, which have been constructed in the hydrographic network since 2001. According to the authors' knowledge, this is the first time that WaTEM/SEDEM was calibrated in Greek catchments under post-fire conditions. The main goal is to estimate the sediment production and transportation rates in the years subsequent to the wildfire and to investigate the potential effect of wildfire on vegetation regeneration and recovery. The calibrated model would undoubtedly serve as a tool to be applied in other Mediterranean watersheds of similar characteristics.

Study Area
The study area consists of two typical Mediterranean catchments located in the suburban forest of Thessaloniki (Seich Sou), the "Konstadinidi" and "Ano Toumba" catchments ( Figure 1, Table 1). The motivation for the selection of these two catchments was the presence and operation of many technical works (check-dams) in the area, the significant changes of vegetation cover caused by the wildfire spread on 7 July 1997 and the easy accessibility of the locations where the technical works have been constructed.   Seich Sou grows at the northeastern limits of Thessaloniki city and constitutes a significant protective zone (natural protective mantle) against flooding and soil erosion. As a whole, the area of the suburban forest is characterized as hilly and semi-mountainous. The topographical relief is intense with several peaks, descending ridges and small ephemeral streams, which flow towards the city of Thessaloniki. The wildfire of 1997 lasted three days and burned 1664 ha of the total area (3018.84 ha) of the suburban forest. Before the fire, most of the study area was cover by dense Pinus brutia forest. Today the same area is cover mainly by regeneration of Pinus brutia, but with younger trees.
Gneiss is the dominant geological formation of the study area, which presents considerable thickness and is mainly impermeable to water. Most gneiss lithological types are easily weathered and covered by a loose weathering mantle of ranging thickness, resulting in the manifestation of springs of usually low yield in its contact with the intact rock [77]. In May 2001, the Forest Research Institute [78] of Thessaloniki published a detailed study analyzing the soil properties of Seich Sou, which was based on field research and lab analyses. In general, the soils of the area can be characterized as shallow, of medium mechanical composition, mainly sandy loam (SL) and with a distribution of horizons A (B) C. The horizon A rarely exceeds 3 cm, while the horizon B is usually absent and is followed directly by the horizon C, which is quite deep. The pH ranges between 6.66 and 7.35 and the soil depth 0-30 cm.
Based on the available meteorological data, the climate of the region could be characterized as a typical Mediterranean climate, with very dry and warm summers and relatively mild winters ( Figure 2). The mean annual rainfall is 437.4 mm, with a maximum in winter (December) and a secondary peak in spring (May). The mean annual temperature is 16.5 • C, the mean max 19.6 • C and the mean min 13.2 • C.

Methodologies
The methodology of the present work can be divided into three main sections: (1) field work and check-dam mapping and dimensioning, (2) determination and analysis of RUSLE parameters and (3) calibration of WaTEM/SEDEM and modeling of annual erosion rates.

Field Work-Check-Dam Trap Efficiency Estimation
The position of each check-dam was located and recorded, while the coordinates of each dam were obtained using GPS device, in the Hellenic Geodetic Reference System

Methodologies
The methodology of the present work can be divided into three main sections: (1) field work and check-dam mapping and dimensioning, (2) determination and analysis of RUSLE parameters and (3) calibration of WaTEM/SEDEM and modeling of annual erosion rates.

Field Work-Check-Dam Trap Efficiency Estimation
The position of each check-dam was located and recorded, while the coordinates of each dam were obtained using GPS device, in the Hellenic Geodetic Reference System "GGRS_1987". Check-dams were constructed during 2001, almost 3 years after the wildfire. Measurements of check-dams' dimensions (height, length, width, etc.) and detailed records of the sediment wedges behind each check-dam were carried out ( Figure 3A). Sediment samples were taken upstream of each check-dam in order to calculate the dry weight of the transported materials retained by the dams. The samples were dried in a laboratory heating chamber at 105 • C for 24 h to calculate the dry weight of the eroded sediments and, therefore, the soil moisture content ( Figure 3B,C).

Methodologies
The methodology of the present work can be divided into three main sections: (1) field work and check-dam mapping and dimensioning, (2) determination and analysis of RUSLE parameters and (3) calibration of WaTEM/SEDEM and modeling of annual erosion rates.

Field Work-Check-Dam Trap Efficiency Estimation
The position of each check-dam was located and recorded, while the coordinates of each dam were obtained using GPS device, in the Hellenic Geodetic Reference System "GGRS_1987". Check-dams were constructed during 2001, almost 3 years after the wildfire. Measurements of check-dams' dimensions (height, length, width, etc.) and detailed records of the sediment wedges behind each check-dam were carried out ( Figure 3A). Sediment samples were taken upstream of each check-dam in order to calculate the dry weight of the transported materials retained by the dams. The samples were dried in a laboratory heating chamber at 105 °C for 24 h to calculate the dry weight of the eroded sediments and, therefore, the soil moisture content ( Figure 3B,C). The recorded data from the field measurements were used to calculate the trap efficiency (TE). The trap efficiency of check-dams was calculated using Brown Equation (1), which estimates TE based on the original capacity of the structure and the sub-catchment of each dam [10,47]. The method is proposed when there are no inflow data available [79]. The Brown equation [80] describes the trap efficiency as follows: where D is a factor determined by the detention time and size of materials and varies between values 1, 0.1 and 0.046 for large, medium and fine-grained materials, respectively. Brown's equation uses a diagram to determine the relation between TE and the C/W ratio [81]. C is the check-dam effective capacity (m 3 ) and W is the check-dam upstream catchment area (km 2 ). The measurement of the catchment area (sub-catchments) for each dam was performed in GIS environment. The estimation of the effective capacity of the check-dams was implemented according to Kotoulas Equation (2) [82]: where h is the dam height (m), tanϕ and tanα are the natural stream bed angle and the slope of depositions, respectively, L is the deposition length at the crown height (m), p and p1 are the two slopes gradient (left and right, respectively, m/m) and m1 is the width of the stream bed (m). In Figure 4 the main parameters of Kotoulas Equation (2) are depicted.
Land 2022, 11, 1705 where h is the dam height (m), tanφ and tanα are the natural stream bed angle and the slope of depositions, respectively, L is the deposition length at the crown height (m), p and p1 are the two slopes gradient (left and right, respectively, m/m) and m1 is the width of the stream bed (m). In Figure 4 the main parameters of Kotoulas Equation (2) are depicted.

RUSLE Parameters
The rainfall erosivity factor (R factor) was calculated, since the WaTEM/SEDEM uses the RUSLE equation to estimate the soil erosion rate. The rain data of the weather station of the Aristotle University of Thessaloniki were processed in a 30 min time step for the period 1990-2021. The equation of Renard et al. [83] (3) was used: where R is the mean annual rainfall erosivity coefficient, E is the rainfall kinetic energy, I30 is the maximum intensity of rain in a 30 min time step, m is the number of rainfalls for each year and n is the number of observation years used to measure the mean R. The EI30 index for each rainfall event was calculated applying Morgan's methodology [84].

RUSLE Parameters
The rainfall erosivity factor (R factor) was calculated, since the WaTEM/SEDEM uses the RUSLE equation to estimate the soil erosion rate. The rain data of the weather station of the Aristotle University of Thessaloniki were processed in a 30 min time step for the period 1990-2021. The equation of Renard et al. [83] (3) was used: where R is the mean annual rainfall erosivity coefficient, E is the rainfall kinetic energy, I 30 is the maximum intensity of rain in a 30 min time step, m is the number of rainfalls for each year and n is the number of observation years used to measure the mean R. The EI 30 index for each rainfall event was calculated applying Morgan's methodology [84]. The 5 × 5 m resolution DEM was used as an input in the WaTEM/SEDEM and for factors estimation. The choice of high-resolution DEM increases the accuracy of the results, especially when dealing with small catchments [85,86]. Furthermore, spatial analysis was performed using GIS and Idrisi [87] software for the preparation of the input data. All layers were created with a 5 × 5 m resolution as follows: (1) Crop factor map (C factor). The determination of the coefficient was made using the CORINE maps in combination with the proposed values for each EU country [88].
The values of C factor used in the current study are presented in Table 2. (2) Soil erodibility map (K factor). The soil erodibility values were determined using the Soil Conservation Service (SCS) nomogram (Renard et al. 1997). For this purpose, maps of mechanical composition and soil hydrological types of the catchments were constructed based on information obtained from pre-existing field surveys [78,89,90]. The soil characteristics and the respective calculated K factors are presented in Table 3.  (3) River map was extracted from the DEM using a threshold value of 1 km 2 , as previous studies suggested [91,92]. The river network was re-digitalized to calculate the river routing. (4) Pond map was created using the calculated trap efficiencies of each check-dam. Checkdams were introduced to the model as sedimentation areas with their respective trap efficiency [10,47,91,92]. (5) Parcels map were estimated using orthophotos of the study area from the past 30 years' period.

WaTEM/SEDEM Setup and Calibration
The spatially distributed model WATEM/SEDEM was applied to calculate the sediments retained by the check-dams and the annual sediment yield at the outlet of the catchments, under the different land-use scenarios during the last decades. The model predicts the long-term values of mean annual erosion caused by runoff and sediment transport [10]. The mean annual erosion is calculated through the modified form of RUSLE proposed by [10,93,94] (4): where E is the mean annual soil erosion (kg m −2 year −1 ), R is the rainfall erosivity factor (MJ mm m −2 h −1 year −1 ), K is the soil erodibility factor (kg h MJ −1 mm −1 ), LS 2D is the two-dimensional topographic factor, C is the crop and management factor and P is the erosion-control practice factor. At the same time, the model enables the calculation of the annual transport capacity T C for each pixel [95], which corresponds to the maximum amount of transported materials that can be transported from one cell to another [66,67] (5): where T C is transport capacity (kg m −2 year −1 ), ktc is transport capacity coefficient (m), R, K, LS are factors of the RUSLE and s is local slope (m m −1 ). During processing the model uses two pathways: when E > T C , the sediments are deposited, while if E < T C , the sediments are transported to the next downstream pixel [47]. According to the above relations, the model creates the modeling routing.
The WaTEM/SEDEM model was calibrated for the entire study area using two different periods of land-use patterns, as well as the measurements of sediment amounts, which were retained by check-dams as the calibration dataset. The land-use patterns of the periods 2001-2013 and 2013-2021 were used for the calibration of the model as they corresponded to the most representative land-use conditions for the period of the check-dams' construction. The check-dams that were used for calibration were non-silted, in order to obtain more accurate sediment yield calculations, as previous studies also suggest [10,47]. Although non-silted dams indicate lower erodibility rates [10], it could also be attributed to a variety of other factors (i.e., rainfall intensity, time of construction, site selection, size Various ratios of ktc (high and low) values were tested to determine model efficiency. Ktc high values are used for poorly-vegetated surfaces, while ktc low values represent high-vegetated surfaces. The high ktc values are assigned for land uses with C factor above the threshold value of 0.1 (ktc limit); on the contrary, land-use classes with lower C factor use a lower transport capacity value (ktc low).
For each different ktc ratio, the WaTEM/SEDEM exported the sediment depositions for each check-dam. These modeled values were compared with the observed using statistical analysis, in order to examine which ktc ratio showed the best model performance. The Nash-Sutcliffe efficiency (NSE) [96] and root mean square error (RMSE)-observations standard deviation ratio (RSR) goodness-of-fit indexes were used to compare the observed modeled pond deposition values, in order to calibrate the WaTEM/SEDEM. NSE ranges between -∞ and 1, where values between 0 and 1 could be considered as acceptable, while 1 is the optimal value. Additionally, values ≤ 0 suggest that the mean of the observed values is a better predictor than the mean of the simulated values, which indicates unacceptable model performance [97]. RSR optimal value is 0, which indicates that RMSE is zero and, therefore, perfect model simulation. Lower RSR means a lower RMSE, which demonstrates better model performance [97].

Check-Dam Trap Efficiency Field Work Results
The main characteristics of the constructed check-dams are presented in Table 4. As it is previously stated, the construction of check-dams was completed in 2001, three years after the wildfire. It is evident from Table 4 that the amount of retained sediments is very low, approximately 4% of the total original storage capacity. After 20 years of operation, it was observed that the check-dams were non-silted, and most of them retained a small amount of sediments. This fact could be attributed to various factors (catchment geomorphology, soil depth, fire severity), but the main cause seems to be the time elapsed between fire occurrence and the check-dams' construction.

Model Calibration
After the examination of various pairs of ktc (low:high), the common ratio of 1:3.33 was finally adopted. At the early stages of this study, different ratios were tested, such as 1:2 1:3, although it resulted in very poor calibration results. Van Rompaey et al. [66] suggested the average ratio of 1:3.33 between the ktc pairs, which were widely used in previous studies [10,47,74,92,98,99]. In the Mediterranean region (Spain) ratio values ranging between 1:3.33 and 1:3 were used [10,79], while in Italy ratios 1:3.80 and 1:2.20 were used in mountainous and non-mountainous catchments, respectively [74,98,100]. Calibration of the transport capacity parameters is still a very important and complex issue that is required to be extensively studied [74,98,101]. In Table 5, the results of the three proposed ktc ratios (1:3.33, 1:3.80 and 1:2.20) are presented, with the respective NSE and RSR values. The evaluation of the optimal ktc values was achieved comparing the observed (from check-dams) and modeled values of the sediment yield, using the NSE and RSR goodnessof-fit statistic indexes. After the determination of the optimal ktc coefficients, they were kept constant for the rest of the calculations, in order to obtain comparable results. The optimal values were 4 for ktc-high and 1.2 for ktc-low, attributing a model efficiency of 0.72 for the NSE statistic index. Previous WaTEM/SEDEM model applications reported a model efficiency (NSE) range between 0.14 and 0.89 [10,47,[101][102][103][104][105]. According to previous researches, NSE values of 0.55 and higher indicate good model performance, and the simulated values are generally acceptable [106][107][108]. Respectively, the RSR showed a value of 0.51, which could be considered satisfactory for model performance [97].

WaTEM/SEDEM Model Application Results
After calibration, the WaTEM/SEDEM model was run for the period 1990-2000 individually for each year. The decade was divided in two land-use periods: pre-fire 1990-1997 conditions and post-fire 1997-2000. Model results showed a total sediment yield of 170 tons in the first period and 955 tons in the period of 3 years subsequent to the forest fire. Re-spectively, the mean annual erosion rate was 0.0419 t/ha/year for the pre-fire period and 0.816 t/ha/year after the fire. As expected, the modeling results presented a significant increase in sediment production in the burnt areas, which was 20 times higher than that of the pre-fire conditions.
A previous study that was conducted in the same area by our research team displays similar annual erosion rates in the undisturbed areas of the suburban forest of Thessaloniki [15]. Diversifications of the erosion rate values among the studies could be explained by the different climate, soil, geomorphological, land cover and experimental conditions [21,[109][110][111][112][113]. In the current study, the difference could be attributed to the outlier value of rainfall intensity and kinetic energy of precipitation in 1995. The estimation of the mean annual erosion rate, excluding the year 1995, presents a value of 0.024 t/ha/year, equivalent to 0.023 t/ha/year of the previous study in the same area [15]. Moreover, due to the high density of forest cover, the soil losses of the study area can be considered negligible. As it is widely known, highly vegetated Mediterranean sites often present low natural soil losses [114,115]. Moreover, it should be mentioned that in some Mediterranean regions, where the soils are often thin, stony and already degraded, low erosion rates are frequently reported [116][117][118]. Additionally, the relatively low values of R factor in the study area may contribute to the low erosion rates. The utilization of the spatial distributed model (WaTEM/SEDEM) provided more reliable quantitative predictions in comparison to RUSLE, which often resulted in overestimation of mean annual soil erosion in the study area. Previous studies conducted in the Mediterranean region showed the same inaccuracy in quantitative predictions with the RUSLE application [119,120].
Comparing the erosion rate values of the two periods, the post-fire mean annual erosion rate was 20 times higher. The increase in sediment production was exclusively affected by the destruction of the forest cover, as the R factor showed a small increase between the years 1997 and 1998 and then remained relatively stable and below the average R (356) in the post-fire years ( Figure 5). The results are also confirmed by the output maps, in which the undisturbed areas did not show an increase in sediment production ( Figure 6). The unburnt area in the Konstadinidi catchment operated as a natural "fence", decreasing the spatial connectivity of sediment production areas with the river network. Many studies have focused on the variation of the erosion rate 1-3 years after a fire, presenting high differences in the estimations [116,[121][122][123]. Such differences were generated mainly by the fire severity, geomorphology and climate conditions of each study area [110,116,124]. 1990 1991 1992 1993 1994 1995 1996 1997 1997-1998 1998-1999 1999-2000 Sediment export (tons) 14 The results are also confirmed by the output maps, in which the undisturbed areas did not show an increase in sediment production ( Figure 6). The unburnt area in the Konstadinidi catchment operated as a natural "fence", decreasing the spatial connectivity of sediment production areas with the river network. Many studies have focused on the variation of the erosion rate 1-3 years after a fire, presenting high differences in the estimations [116,[121][122][123]. Such differences were generated mainly by the fire severity, geomorphology and climate conditions of each study area [110,116,124]. The results are also confirmed by the output maps, in which the undisturbed areas did not show an increase in sediment production ( Figure 6). The unburnt area in the Konstadinidi catchment operated as a natural "fence", decreasing the spatial connectivity of sediment production areas with the river network. Many studies have focused on the variation of the erosion rate 1-3 years after a fire, presenting high differences in the estimations [116,[121][122][123]. Such differences were generated mainly by the fire severity, geomorphology and climate conditions of each study area [110,116,124]. Excluding the period 1997-1998, in which the R factor was exceptionally very low, the erosion rate in the burned area was estimated to be 0.998 t/ha/year. Post-fire erosion rates in Mediterranean regions often show lower values compared to other climate zones [116,125]. The results of the calibrated model comply with this assumption, since a 0.998 t/ha/year erosion rate is considered to be close to the lower limits of erosion in burned Mediterranean catchments, which usually range between 1-10 t/ha/year, depending on the fire severity and other related factors [110,116,126]. As aforementioned, the soil depth in Seich Sou forest is very thin, stony and already significantly eroded, a fact which explains the low erosion values in all periods. In addition, the mean value of the post-fire erosion rate lies within the acceptable limits of Europe's soil formation (0.3 and 1.4 t/ha/year) [126,127]. However, when calculating the severity of the erosion rate, it should be taken into account that the average rate of soil formation in the Mediterranean region is considered to be relatively slower (0.1 t/ha/year) [110], imposing the need to validate modeling results using field measurements.
Concerning the period 2001-2021, the erosion rates, as expected, showed a decrease. The alteration on the sediments' mobility was caused by the check-dams' construction and the increase of the vegetation cover rate. Previous studies had extensively analyzed the effect of vegetation cover and check-dams in the reduction of soil erosion rates [10,28,32,[128][129][130]. The period 2001-2012 modeling results demonstrated an erosion rate equal to 0.084 t/ha/year. Respectively, the estimation of mean annual erosion rate, for the period 2013-2021, presented a value of 0.08 t/ha/year. Although the vegetation cover the second period has higher density, the erosion rate presents similar values. This condition could be attributed to the higher rainfall kinetic energy of the period 2013-2021.
Concerning the long-term effects of the forest fire on soil erosion, it should be highlighted that 25 years after the fire, the mean erosion rate remains double (0.08 t/ha/year) in comparison with the pre-fire period (0.0419 t/ha/year). Although forest regeneration was very fast and successful in the study area, it seems that the forest ecosystem has not been fully restored to the pre-fire condition. In adverse Mediterranean conditions, the post-fire "window of disturbance" could be significantly prolonged, causing increased soil erosion rates [15,131].

Post-Fire Modeling Application with Check-Dams
The model was also applied for 3 years, under the scenario of check-dam construction immediately after the wildfire, in order to investigate the importance of fast decisionmaking after the fire and the impact of check-dams on soil erosion. Previous studies have highlighted the vulnerability of soil erosion and degradation in the first weeks after the forest fire [132,133], caused by the climatic variables of precipitation (intensity, distribution) under Mediterranean conditions. Model results showed a total amount of 548 tons 3 years after the fire with the check-dams, corresponding to a 0.36 t/ha/year mean annual erosion rate. The check-dams retained 411 tons of sediments, reducing by 42.61% the sediments' export of the study area. Boix-Fayos et al. [10], respectively, presented 34% less sediment delivery caused by the construction of check-dams.
It must be highlighted that in the first year after the wildfire (1997), the rainfall kinetic energy was exceptionally low, a fact which led to less sediment formation and transport. The sediment deposition behind the check-dams for each year individually was 81, 176 and 154 tons for 1997-1998, 1998-1999 and 1999-2000, respectively. The first year R factor was 102.03 MJ mm m −2 h −1 year −1 , followed by two years where the values were the double 228.31 MJ mm m −2 h −1 year −1 and 199.95 MJ mm m −2 h −1 year −1 , respectively. Similar results caused by lower rainfall kinetic energy were presented also by Shakesby et al. [117], in which in the first year after the fire the erosion rate was 1.38 t/ha/year, and on the contrary, the second-year erosion rate was 3.57 t/ha/year. Additionally, Bussi et al. [134] highlighted the significance of rainfall intensity on sediment transport, revealing that 40% of the sediment volume retained behind the check-dam was delivered from one rainfall event of high intensity.
Generally, the construction of check-dams decreases soil erosion rates and is a positive, drastic and efficacious management practice for post-fire sediment control [10,47,135,136]. However, check-dams were characterized as a short-term solution due to their limited life expectancy, referring to their finite storage capacities [47,137].

Conclusions
In the present study, the spatial distributed model WaTEM/SEDEM was calibrated based on field measurements derived from 18 check-dams in the suburban forest of Thessaloniki. The overall aim of this research was to estimate the annual soil erosion rates since 1990 under significant land-use changes (including the catastrophic impact of the wildfire in 1997) and the effect of the constructed check-dams as an erosion-control measure. The results showed very a low erosion rate for the pre-fire period, an expected significant increase in the 3 years subsequent to the wildfire and a constant decrease in the years until 2021. However, it seems that the post-fire regeneration of the forest has not been fully achieved, since the annual soil erosion rate in the post-fire period is double compared to the pre-fire period. Twenty-five years after the forest fire, the "window of disturbance" is still open in the study area.
The kinetic energy of rainfall proved to be an important factor in the accuracy of sediment export. Estimation of annual erosion rates, on a case-by-case basis, should use the analytical R factor method to avoid miscalculations, due to the fluctuation of the kinetic energy of rainfall events.
Check-dams are considered to be an efficient short-term solution for sediment control. In the study area, the scenario of the construction of check-dams immediately after the wildfire presented a 42.61% decrease on the annual erosion rate. Nevertheless, the checkdams were constructed 3 years after the fire, leaving a significant time gap. Our indication is Land 2022, 11, 1705 13 of 18 validated by the field measurements, where the amounts of the sediments retained behind the check-dams were very low. In addition, the calibration and validation of erosion models using field measurements enhance the accuracy of the modeled erosion values. Validated erosion models can be used as a tool for fast decision-making of suitable erosion-control works' determination and establishment, actions which are extremely crucial during the early post-fire period. The calibrated WaTEM/SEDEM model can be applied successfully in other Mediterranean catchments of similar characteristics.

Data Availability Statement:
The data that support the findings of this study are available from the corresponding author upon reasonable request.

Conflicts of Interest:
The authors declare no conflict of interest.