Hydrological Response of Natural Mediterranean Watersheds to Forest Fires

Abstract

Fires are common in forested Mediterranean-climate watersheds. Forest fires cause abrupt land use/cover (LULC) changes affecting soil properties and hydrological processes within and across watersheds. A major forest fire in Attica, Greece, that affected the Lykorrema stream experimental watershed provided the opportunity for an in-depth study of the impact of forest fires on the hydrological balance of natural Mediterranean watersheds. To this end, detailed hydrometeorological data recorded for five years before and for five years after the fire incidence were utilized. SWAT model was also used to consider the potential influence of meteorological conditions temporal variability on the results of the analysis. Specifically, SWAT model was parameterized calibrated and validated for the pre-fire and the post-fire conditions using the corresponding detailed hydrometeorological data for the respective periods. Then the two versions of the model were applied for the entire period providing comprehensive time series for all the flows and storages in the studied watershed. In this way, the post-fire LULC and soil properties changes were the only influencing factors driving the alterations in the hydrological balance allowing an impartial comparison. The obtained results highlighted the considerable impact of forest fires on the watersheds’ hydrological functioning. Specifically, the maximum direct runoff depths and the maximum flow rates were substantially higher in the post-fire conditions. In contrast, actual evapotranspiration was reduced, when the effect of fire was considered. The obtained results indicate that the altered post-fire LULC and soil properties are major drivers of the watershed’s hydrological balance changes. SWAT model performed sufficiently well for both the pre- and post-fire conditions and provided a deeper insight into the impact of forest fires on the hydrological functioning of natural Mediterranean watersheds.

1. Introduction

One of the most common hazards of the forested watersheds in the Mediterranean basin is forest fires. Anthropogenic activities are considered as a driving force of wildfires, however, the lack of forest monitoring and climate change play also a significant role [1,2,3]. In 2018 the Intergovernmental Panel on Climate Change (IPCC) reported that human-induced warming has already increased to 1 °C above pre-industrial levels. Studies have shown that with increasing warming pattern the burnt area is also expected to increase [4].

Forest fires have numerous severe outcomes. Alterations in vegetation cover, soil properties disturbance, changes in the hydrological behavior of watersheds as well as probable life and property loss are some of the results of a forest fire [5]. Loss in vegetation cover has a direct impact in the hydrological response of a watershed with increased overland flow caused by the reduction of infiltration as well as through decreased evapotranspiration [6,7,8,9,10,11,12,13,14,15]. The hydrological response of an area can be also indirectly affected by the disturbances of the soil properties of the area. The destruction of the top soil organic matter that leads to the destabilization of the soil structure and the increased ground coverage from ash may exacerbate the impacts of water repellency [1,16,17]. This phenomenon leads to the reduction of infiltration and subsequently to the increase of overland flow [1,18]. The decrease of evapotranspiration may also result in increased groundwater recharge [19,20,21].

The effect of forest fires on the hydrological behavior of a certain area is obvious for around one to two years after the fire event [1,18,22]. However, several studies have revealed that the impact of the forest fires in the hydrological behavior can be apparent even five to seven years after the fire incident [1,6,18,23,24]. The required time period for the hydrological processes of a watershed to return to pre-fire conditions is strongly dependent on the rate by which the vegetation cover is restored. In more detail, the amount of time needed for full runoff restoration depends on the different recovery rates of each species existing before the fire in the area [25,26,27]. However, influencing factors of the hydrological process restoration rate are the meteorological and hydrological conditions of the area as well as the erodibility of the soils and the steepness of the slopes [23,28,29]. The results of Wine et al. [30] indicated that Mediterranean ecoregion divisions may present relatively fast post-fire recovery rates.

In a post-fire regime, the difference in the water balance can be the result of land use changes or/and of climatic interferences that may occur across a range of time scales [21,30,31,32]. Several studies have been conducted in order to determine how forest fires affect the hydrology of natural river basins. As an example, Zhou et al. [18] in their study concluded that the hydrological behavior of the studied watershed in southeast Australia was severely altered because of wildfires. More specifically, for the first 15 years after the fire event the streamflow of the watershed increased. This resulted from the decrease in evapotranspiration and soil infiltration. Another study that draws the same conclusions regarding the effect of forest fires is the one conducted by Pereira et al. [13]. In the period after the fire incident, the peak runoff resulting from heavy rainfall is higher than in the pre-fire period as well as the runoff during the dry periods is lower than in the pre-fire period. However, they concluded that the average increase in daily runoff was around 5% on average and 20% at most.

Generally, the analysis of the effect of forest fires on the hydrological response of natural watersheds is very challenging. The main difficulties are: a) the lack of detailed hydrometeorological data for adequate periods for both the pre-fire and the post-fire periods and b) the effect of meteorological conditions variability that makes the two periods not directly comparable as the postfire recovery period is much shorter than the typical 30 years reference period [18,33]. Additionally, the effects of climate variability could be confused for wildfire effects. Lagged watershed responses may also complicate any further analysis [21,30]. For this reason, hydrological modeling is routinely used to evaluate the effect of forest fires on the hydrological response of natural watersheds and to assess the associated hydrological risks (e.g. [12,13,14,15,18,34,35,36,37,38,39,40,41,42,43,44,45,46]. However, in most relevant studies, model parameterization for the post-fire conditions is made arbitrary depending on the fire severity or the expected land cover changes, instead of calibrating and validating the model based on actual hydrological observations for the pre-fire and post-fire periods [18].

Earth Observation (EO) satellites provide multitemporal, high-resolution data, with systematic revisit and synoptic view of the earth surface in local, regional, and national scale. Therefore, they consist of an accurate mean to ideally map the LULC temporal changes before and after fire events [47,48].

A major forest fire in Attica, Greece, which occurred in August 2009 and affected the Lykorrema stream experimental watershed [1,17,33,48,49,50,51,52,53] provided the opportunity to study the effect of forest fires on the hydrological balance of natural watersheds utilizing detailed hydrometeorological data obtained both before and after the fire incidence. A preliminary assessment carried out by Soulis et al. [1] in the same study area soon after the forest fire, revealed an apparent increment of direct runoff depths and peak flow discharges. In this study, hydrological modeling and detailed hydrological and meteorological data that are available for the study area for the time period between 2004 and 2014 covering almost equally the pre-fire and post-fire conditions are utilized to make an in-depth investigation of the impact of forest fires on the hydrological processes and the hydrological balance of natural Mediterranean watersheds. In order to investigate the effect of the post-fire land use changes as the only influencing factor of the alterations in the water balance, the soil water assessment tool (SWAT) model was used. Specifically, SWAT model was calibrated and validated for the pre-fire and the post-fire conditions using the corresponding detailed hydrometeorological data for both periods. Then the two versions of the model were executed for the entire period providing comprehensive time series for all the flows and storages in the studied watershed. In this way, the post fire land use and soil properties changes were the only influencing factors of the alterations in the water balance allowing an impartial comparison. In this way it was made possible to consider the influence of meteorological conditions variability in the pre-fire and the post-fire period, to provide a deeper insight on the forest fire impact on the hydrological processes and to quantify its effect on the hydrological balance components.

2. Materials and Methods

2.1. Study Area and Input Data

The study was executed based on the small-scale natural Mediterranean experimental watershed of Lykorema stream. The watershed lays in the east side of Pedeli Mountain in Attica, Greece (coordinates: UL 23°53′33″ E–38°04′13″ N, LR 23°56′00″ E–38°02′28″ N) and it covers an area of 7.84 km² (Figure 1). The elevation of the watershed ranges between 950 m.a.s.l and 280 m.a.s.l with an average elevation of 560 m.a.s.l. The average slope of the basin is 36%, while overall the watershed is characterized by a relatively sharp relief.

The climate of Attica, where the basin is located, is characterized as Mediterranean with mild, wet winters and hot, dry summers. The biggest part of precipitation falls in the autumn–spring period in liquid state, with average rainfall for the study period (2004–2014) being around 668 mm. The average maximum temperature of 31.3 °C occurs during July while the average minimum temperature of 3.4 °C occurs in January. The average annual temperature for the study time period was 16.3 °C. Because of the small size of the studied headwater watershed, its steep relief and the limited contribution of groundwater flow from deep aquifers, the watershed responds promptly to rainfall inputs [53].

The dominating geological formation of the watershed is schists covering 96% of the area, while the rest is covered by marbles. The geological formations of the basin are strongly tectonically fractured with eroded upper layers [54]. Soulis et al. [17,53] executed a detailed study of the soil characteristics of the watershed including an analysis of the spatial and temporal variability of hydraulic conductivity during the study period [17]. In general, it is described that the basin is mainly covered by coarse soils while smaller parts with medium texture soils, both with high hydraulic conductivities (Figure 1). Soil hydraulic conductivity in the studied watershed was profoundly decreased after the effect of the forest fire of 2009, but the measurements seven years after the incidence were similar to the pre-fire period [17].

The Lykorema stream experimental watershed is operated by the Agricultural University of Athens and it is instrumented and systemically studied since September 2004. The hydro-meteorological network that covers the area is relatively dense. It includes one hydrometric station at the outlet of the basin and it encloses five rain gauges, four temperature-relative humidity recorders, and one full meteorological station with a soil moisture sensor array (Figure 1). The hydrometric station in the outlet of the basin consisted of a sharp crested-weir accompanied by a digital water-level recorder. Soon after the fire of 2009, the station was expanded with a larger broad-crested weir, in order to cover the higher flow rates and the increased sediment movements. Furthermore, a digital water-level recorder was installed at the middle of a long rectangular culvert just before the weir to secure the measurement of any unexpectedly high peak flows. More details about the hydrometeorological network and the calibration of the hydrometric stations are provided by Soulis and Dercas [49] and Soulis [33].

The meteorological input for the study was provided by four rain gauges and one full meteorological station, which measured the data for rainfall, temperature, wind speed, relative humidity, and solar radiation. All stations are distributed within the watershed, perform measurements with 10 minutes temporal step and were provided by the Agricultural University of Athens (AUA) and the National Technical University of Athens (NTUA). A quality check of the meteorological data was performed including the analysis of the correlation between nearby stations and the identification of gaps and obviously erroneous values, such as unrealistic values or values beyond specific thresholds, considering also the metadata about malfunctions in the monitoring network. Missing data and erroneous values within the meteorological time series were reproduced via spatial interpolation from the available data of the nearest station. The temporal step of the dataset used in this study is daily and it covers the time period between 2001 and 2014.

Lykorema stream watershed was selected for this study, as it experienced a major wildfire incident in 2009 during the operation of the experimental watershed. This wildfire affected the larger part of the watershed (Figure 1) and provided the opportunity for an in-depth study of the impact of forest fires on the hydrological balance of natural Mediterranean watersheds utilizing the detailed hydrometeorological data obtained for adequate time before and after the fire incidence.

The land cover classification for the conditions before and after the wildfire of 2009 was performed with the assistance of remote sensing and extensive field surveys. More specifically, images from Landsat-5 Thematic Mapper satellite system were acquired through the Earth Explorer system of the United States Geological Survey (https://earthexplorer.usgs.gov/) [48]. The Maximum Likelihood classification technique was applied using ground control points collected from the field survey. The boundaries of the obtained classes were refined using very high-resolution orthophoto maps delivering a highly accurate LULC map.

Consequently, the watershed was covered mainly by a mixed vegetation cover consisting of young pine trees, pasture, and a few scattered tufts of old trees. A small part of the basin was covered by bare rock, as well as from a network of forest roads [50,53]. After the fire incident of 2009, almost 87% of the basin area was affected, because of the flammable land cover and the weather conditions at the time. However, because of the topography of the basin, consisting of steep slopes, there are parts of the basin that remained unburnt or semi-burnt [33].

2.2. Model

Hydrological modeling is routinely used to evaluate the effect of forest fires on the hydrological response of natural watersheds [18,34,35,36,37,38,39,40,41,42,43]. In this study, the soil and water assessment tool (SWAT) was used to quantify the effect of forest fires on the various hydrological balance components taking into consideration the potential influence of the different meteorological conditions prevailing the period before and after the fire incident on the results of the current analysis. Specifically, the most recent version of ArcSWAT 2012 interface (https://swat.tamu.edu/software/arcswat/) for ArcGIS 10.5.l (ESRI) was used and all geospatial analysis was also made in the later GIS software platform.

SWAT is a process-based semi-distributed, conceptual and continuous hydrological model which is developed to predict the impacts of land management practices on the environment variables in large complex watersheds with varying soils, land use, and management conditions over long periods of time [55]. It was selected because it is a robust, well-established, public domain model that can be used for this purpose and that is very well suited in the case of natural mountainous watersheds. SWAT model has been applied successfully in many relevant studies [56,57,58,59,60,61,62,63].

SWAT divides the basin into several subbasins to achieve better simulation of the physical processes in the watershed. The use of subbasins is of great benefit particularly in watershed areas that are dominated by different land uses soil types and hydrological behavior. Furthermore, SWAT divides the watershed into hydrological response units (HRUs), which are lumped areas within the watershed with the same land use, soil type, and slope. The model simulates all hydrological processes in HRU level. In order to get accurate predictions on the studied process, SWAT first investigates the water balance of the studied area. Water balance is the driving force of all the processes in the watershed area [55].

2.2.1. Model Setup

The SWAT model was set up with two different configurations during this study. In the first case, the model configured based on the information of the watershed before the fire incident (Oct 2004–Aug 2009) and in the second case regarding the basin’s information after the fire (Aug 2009–Dec 2014). For both set-ups, the same digital elevation model (DEM) was occupied. The DEM was generated using the corresponding 1/5000 scale topographic diagrams of the Hellenic Military Geographical Service and has a spatial resolution of 5 m × 5 m. Moreover, the flow direction and flow accumulation were computed for both set-ups with the same threshold, which resulted in the division of the whole area into 69 subbasins. As a next step land use, soil and slope maps were imported to the model, which were required for the formation of the HRUs. Each model had different land use and soil maps corresponding to the pre-fire and the post-fire conditions. This whole procedure resulted in the formation of 244 HRUs. Furthermore, the meteorological information of each period was imported to the two models. In both models Penman-Monteith equation [64,65] was selected as the reference evapotranspiration calculation method, while surface runoff and infiltration were calculated through the SCS-CN methodology [33,48,51,52,66]. Lastly in both models the first 3 years were used for model warm up, while the model output was provided in daily time step.

One important aspect of the whole model setup procedure was the assignment of the land use properties of the Lykorema watershed to the corresponding land uses of the USA region. SWAT provides a table of properties for different land use types [55]. Knowing the land use properties of our watershed we assigned each land use type to the corresponding type of the SWAT table. This way it was made certain that the land use input was prepared exactly as the model needed it. Regarding the soil properties, it was needed to enrich the SWAT soil table with three new soil classes that were representing the soil classes of the Lykorema watershed. Based on the detailed information over the Lykorema watershed’s soils obtained in previous studies [17,53], the sufficient enrichment of the SWAT table was performed.

2.2.2. Sensitivity Analysis, Calibration, and Validation

Sensitivity analysis, calibration, and validation were carried out with the assistance of the SWAT Calibration and Uncertainty Procedures (SWAT-CUP) software. SWAT-CUP (https://swat.tamu.edu/software/swat-cup/) is a calibration tool fully compatible with the output file format of SWAT, thus is easy to implement when using SWAT as the hydrological model. For this study SUFI-2 [67] was selected as the calibration algorithm. SUFI-2 serves as a combination of the calibration and the uncertainty analysis and thus parameter uncertainty is considered for all sources of uncertainty, like conceptual model, driving variable or measured data.

The first step was to perform a sensitivity analysis for both models. Sensitivity analysis (SA) aims to identify the key parameters that affect model performance and it plays important roles in model parameterization, calibration, optimization, and uncertainty quantification [68]. In our study, the one-at-a-time sensitivity analysis method was used. Through sensitivity analysis model parameters affecting the surface runoff, the groundwater, the baseflow, the evapotranspiration, the soil properties, and the response of the basin were tested. As a result, the group of sensitive parameters were identified for the calibration procedure.

An important fact is that the curve number (CN) was not included in the sensitivity analysis and the calibration procedure, even if it is a key influencing parameter for the model. Specifically, for the CN parameter the spatially distributed CN values that were directly estimated using measured rainfall and runoff data by Soulis [33] for the studied watershed for the pre-fire and the postfire conditions using the Soulis and Valiantzas method the [51,52] were utilized to reduce the number of calibration parameters. For each soil and land cover complex of the watershed, the corresponding CN number was assigned as described in the study of Soulis [33].

After defining the sensitive parameters for the model, the next step was the model calibration. For both models, three years were used as the calibration period. For this time period, the parameters that resulted from the sensitivity analysis were occupied and were assigned with physical ranges. For the execution of the calibration, several iterations of 2000 simulation runs each were performed until adequate model performance to be achieved. It should be noted that three different calibration approaches were attempted and compared in this study to facilitate a more representative comparison between the pre-fire and post-fire conditions. In the first attempt, the two model configurations were calibrated independently. In the second attempt, the pre-fire configuration of the model was initially calibrated. Then the post-fire model was calibrated using the parameter values that were specified for the pre-fire configuration and allowing as free parameters only those that could be possibly affected by the forest fire (e.g., soil hydraulic conductivity and parameters related with the vegetation cover) and only for the affected parts of the watershed. As it was explained above the CN values of the affected areas were also changed accordingly. In the third attempt, the postfire configuration of the model was initially calibrated. Then the pre-fire configuration of the model was tested with the parameter values that were specified for the post-fire configuration. The most suitable calibration approach was selected by comparing the performance achieved by each one of them and considering as well their ability to provide representative results for the two examined scenarios. As it is explained in detail in Section 3, the most representative model configurations were obtained with the second approach.

The next step was the model validation to verify the ability of the models to adequately describe the pre-fire and post-fire hydrological functioning [69]. In order to execute the model validation, the parameter values that were specified during the calibration procedure were occupied. Each model was run for two years and the assessment of the model performance was performed.

Model performance, during the calibration and validation procedures, was evaluated statistically based on the Nash–Sutcliffe efficiency and the modified Nash–Sutcliffe efficiency (NSE and MNSE, −∞ to +1, values close to +1 indicate better model performance) and the percent bias (PBIAS, indicator of under- or over-estimation, values close to 0 indicate better model performance). The performance was considered satisfactory if NSE > 0.5 and PBIAS < ±25% [70].

2.2.3. Model Application and Data Analysis

The main goal of this study is to investigate the impact of forest fires on the hydrological balance of natural watersheds. The study utilizes 10 years of detailed hydrometeorological data, 5 before and 5 after the fire incident, as well as detailed field and physiographic data. However, impartial comparison of the hydrological balance for the two periods is not feasible as the weather conditions for these short time periods before and after the fire vary making it very difficult to differentiate the effect of the forest fire from the effect of weather conditions variability.

Accordingly, in order to study the effect of the post-fire changes as the only influencing factor of the alterations in the water balance, the two configurations of SWAT model that were calibrated and validated for the pre-fire and the post-fire periods were applied from 2001 to 2014 covering the entire study period (2004–2014) and three additional years as the warmup period (2001–2003). In this way, it was made possible to make an impartial comparison of the hydrological balance for the two periods considering the influence of the weather conditions variability in the pre-fire and the post-fire period and additionally to provide a deeper insight on the forest fire impact on the hydrological processes and to quantify its effect on the hydrological balance components.

The obtained 10-year data series for the main hydrological balance components (discharge, surface, lateral and groundwater flow, actual evapotranspiration, and percolation) were statistically analyzed and compared in daily, monthly, and annual time steps. The flow duration curves (FDC) for the two simulated data series and the observed data were created and compared as well. FDC is a plot of discharge versus the percentage of time that each particular discharge is equaled or exceeded (exceedance probability).

3. Results

3.1. Model Calibration and Validation

To obtain a general overview of the behavior of the two model configurations, the initial model performance without calibration was examined. Both the pre-fire model and post-fire model configurations performed poorly without calibration. Specifically, they show significant overestimation of the high flows and underestimation of baseflow. This information provided a solid ground to further parametrize the models in a sensible way during the calibration procedure.

Next step was the sensitivity analysis. Numerous parameters, which could potentially affect the direct runoff response and the base flow (e.g. surface, lateral and groundwater flow, deep percolation, and evapotranspiration processes), were tested. As stated before, a one-at-a-time sensitivity analysis was performed in order to detect the effect of each parameter alone, without the influence of any other. The sensitivity analysis produced identical sensitive parameters for both models. The parameters that were specified to be used in the calibration procedure are presented in

Table 1. The more sensitive parameters of the model and their ranges as they resulted by the sensitivity analysis.

Parameter Name	Parameter Description	Value Change Method	Range
ESCO	Soil evaporation compensation factor.	replace	(0, 1)
EPCO	Plant uptake compensation factor.	replace	(0, 1)
SOL_AWC	Available water capacity of the soil layer.	relative	(−0.5, 0.5)
SOL_K	Saturated hydraulic conductivity.	relative	(−0.5, 0.8)
SOL_Z	Depth from soil surface to bottom of layer	relative	(−0.5, 1.5)
ALPHA_BF	Baseflow alpha factor.	absolute	(−0.048, 0.95)
GDRAIN	Drain tile lag time.	absolute	(0, 100)
GW_DELAY	Groundwater delay.	absolute	(−10, 300)
GW_REVAP	Groundwater "revap" coefficient.	absolute	(−0.02, 0.18)
GWQMN	Threshold depth of water in the shallow aquifer required for return flow to occur.	absolute	(−500, 2000)
RCHRG_DP	Deep aquifer percolation fraction.	absolute	(−0.05, 0.95)
REVAPMN	Threshold depth of water in the shallow aquifer for “revap” to occur.	relative	(−0.6, 0.4)
LATTIME	Lateral flow travel time.	absolute	(0, 180)
CNCOEF	Plant ET curve number coefficient.	absolute	(−0.5, 1)

.

However, during the model calibration procedure, it was observed that the addition of some additional parameters could improve the model performance, such as depth to the impermeable layer (DEP_IMP), and delay of the lateral flow (LAT_TIME). It was also observed that the calibration of the soil hydraulic properties parameters independently for each soil type achieved superior performance than their calibration in lumped form.

The calibration of the CN parameter was also attempted to examine if the adjustment of the CN values directly estimated by rainfall and runoff measurements [33] could provide a noticeable improvement in the models’ performance. The obtained results indicated that the optimum calibrated values of the CN parameter for each soil-cover complex slightly changed (less than ±1). Any improvements in model performance were also negligible.

For the pre-fire configuration of the model, the available period for calibration and validation was from Oct 2004 until the fire event in Aug 2009. This time period included almost 5 hydrological years, of which the first 2 years were wet and the last 2 years were very dry. Therefore, the hydrological years 2004–2005, 2007–2008, and 2008–2009 were used for the calibration and the hydrological years 2005–2006 and 2006–2007 for validation in order to include at least one wet year for both calibration and validation and to increase the homogeneity between the two periods. For the post-fire model configuration, the available time period for calibration and validation was from Oct 2009 until Dec 2014, covering almost 5.5 hydrological years. In this case, the distribution of wet and dry years was better compared to the pre-fire period. Therefore, for simplicity, the first 3.5 hydrological years were randomly selected as the calibration period, while the remaining 2 as validation period.

The model performance, for calibration and validation, and its suitability for the purposes of this study were evaluated both quantitatively based on the performance metrics and qualitatively based on a graphical representation of the simulated versus the observed discharge. The model performance metrics obtained with the three approaches described in the methodology section are presented in

Table 2. Model performance metrics for the three calibration approaches, (a) independent calibration of both configurations, (b) independent calibration of the pre-fire model and calibration of only the parameters affected by the forest fire for the affected parts of the watershed (using pre-fire model calibrated parameters), (c) independent calibration of the post-fire model and use of the common parameters for the pre-fire model (using post-fire model calibrated parameters). NSE is the Nash–Sutcliffe efficiency MNSE is the modified Nash–Sutcliffe efficiency and PBIAS is the percent bias.

	Pre-Fire Model		Post-Fire Model
	Independent Calibration (a,b)	Using Post-Fire Model Calibrated Parameters (c)	Independent Calibration (a)	Using Independent Pre-fire Model Calibrated Parameters (b)
Calibration
NSE	0.63	0.47	0.84	0.82
MNSE	0.54	0.33	0.61	0.57
PBIAS (%)	−2.8	17.6	4.8	−5.6
Validation
NSE	0.56	0.42	0.73	0.7
MNSE	0.51	0.28	0.57	0.53
PBIAS (%)	−3.9	13.9	6.3	−7.9

.

As can be seen in this table the obtained NSE values for the independent calibration and the corresponding validation of the pre-fire model configuration can be considered as satisfactory. The obtained PBIAS values were also sufficient while the MNSE values were only acceptable indicating that estimation of the low flows was the most challenging. It can be also seen that using the calibrated parameters from the post-fire configuration the performance of the model is not acceptable; however, it is interesting that it is not very far from that. As regards the post-fire configuration, the performance achieved for the independent calibration is very good while performance achieved with the calibration of only the parameters that could possibly be affected by the forest fire and only for the affected parts of the watershed was only slightly inferior.

Based on the above results, the model configurations obtained with the second approach, i.e., the independent calibration of the pre-fire model and calibration of only the parameters affected by the forest fire for the affected parts of the watershed for the post-fire model, was selected for the investigation of the impact of forest fires on the hydrological balance of the studied watershed even if it obtained slightly inferior performance for the post-fire period.

The main parameters that varied in the post-fire version of the model were the saturated soil hydraulic conductivity (SOL_K, slightly reduced), the available water capacity of the soil layer (SOL_AWC, slightly reduced), the lateral flow travel time (LAT_TTIME, reduced), the soil evaporation compensation factor (ESCO, increased), the plant uptake compensation factor (EPCO, reduced), and the curve number (CN, increased) for the two soil types of the affected area. The new parameters values are related with reduced infiltration, faster lateral flow, and reduced evapotranspiration from deeper layers of the soil, which is following the expected behavior of the model for the post-fire conditions.

Using similar parameterization for both periods at least for the characteristics and the parts of the watershed that were not affected by the fire, may provide a better comparison between the two periods and a better insight on the changes of the hydrological processes due to the fire. Considering also that precipitation was higher during the post-fire period, the selected calibration approach verifies additionally the ability of the model that was calibrated for the pre-fire period to provide adequate predictions for the higher precipitation depths observed during the post-fire period. A graphical representation of the simulated daily discharge obtained by the selected calibration approach versus the observed discharge is provided in Figure 2.

As can be seen in the graphical representation of the simulated versus the observed discharge for the pre-fire period (Figure 2a), the simulated discharge follows closely enough the observed one; however, in some cases, the model was unable to reproduce the peak flows. In contrast, the model predicts more accurately the peak flows for the post-fire period (Figure 2b) and follows closely the observed hydrograph as well, resulting in superior performance for the post-fire period. Based on the quantitative and the qualitative evaluation it can be deduced that the two configurations of the model reproduce sufficiently the general patterns of discharge variation for the pre-fire and the post-fire conditions.

3.2. Hydrological Response

From the analysis of the observed rainfall-runoff data for the pre-fire and the post-fire periods, it can be deduced that the wildfire of August 2009 had a considerable effect on the hydrological response of the case study watershed. As can be seen in Figure 2, the observed daily peak discharges are much higher in the post-fire period. Specifically, the maximum observed daily discharges during the pre-fire and the post-fire periods were 0.54 and 4.94 m³/s correspondingly. At the same time, the maximum daily rainfall depths during the pre-fire and the post-fire periods were 111.4 and 149 mm correspondingly. In the same figure, it can be also seen that similar rainfall depths resulted in much higher discharges in the post-fire period. Increased discharges are observed during the entire studied post-fire period.

A similar picture is obtained when the maximum hourly discharges are considered (Figure 3). Specifically, the maximum peak discharge observed after the forest fire was 81.25 m³/s, which was almost 11.8 times higher than the maximum peak discharge observed before the fire (6.9 m³/s). This behavior is more obvious in Figure 3 where the peak hourly discharges for the events that took place in the pre-fire and the post-fire period are plotted in comparison with the total rainfall depth and the maximum hourly rainfall intensity for each event. As can be seen, even if there is substantial scatter, it is obvious that higher peak discharges are observed for similar rainfall depths and similar rainfall intensities. Soulis et al. [53] concluded that higher runoff values are not correlated with increased antecedent soil moisture; accordingly, higher peak discharges cannot be attributed solely to more wet initial conditions during the post-fire period.

However, the observed deviations of the weather conditions before and after the fire (Figure 2) hamper the comparison and make the deduction of safe conclusions very difficult. To this end, the model configurations corresponding to the pre-fire and the post-fire conditions were applied for the entire 10-years study period (Oct 2004–Dec 2014). Three and half additional years (2001–2004) were used as warmup period. The obtained results in daily time step are summarized in

Table 3. Statistics of the observed and the simulated daily discharges. Q1 to Q95 are the percentile flows, which are the flows equaled or exceeded at least at the corresponding percentage of the time.

	Observed Entire Period	Observed Pre-Fire	Observed Post-Fire	Simulated Entire Period Pre-Fire Conditions	Simulated Entire Period Post-Fire Conditions
Statistic	Discharge (m3/s)
Average	0.027	0.019	0.033	0.020	0.034
Maximum	4.942	0.616	4.942	1.343	4.801
Minimum	0.003	0.003	0.004	0.001	0.002
Q95	0.004	0.004	0.005	0.002	0.003
Q90	0.005	0.004	0.006	0.003	0.004
Q75	0.007	0.006	0.009	0.005	0.006
Median (Q50)	0.013	0.010	0.017	0.010	0.010
Q25	0.025	0.015	0.031	0.022	0.020
Q10	0.047	0.045	0.049	0.047	0.048
Q5	0.069	0.071	0.069	0.064	0.117
Q1	0.154	0.130	0.219	0.132	0.388

. As can be seen the maximum and the average daily discharges are profoundly higher in the post-fire conditions. In general, obvious differences are observed in higher flows; specifically, flows over the 5th percentile flow (flow equaled or exceeded at least 5% of the time). The statistics of the observed discharges for the entire period lay somewhere between the simulated flows for the two investigated scenarios. However, both model configurations seem to underestimate the lower flows.

A clearer picture is presented in the flow duration curves depicted in Figure 4. All depicted FDCs are generally steep which is reasonable for a small watershed in a Mediterranean climate. The FDC for the post-fire conditions (Figure 4b), though, is steeper indicating significant high flows for short periods of time and higher flow variability. The middle and lower parts of the FDCs corresponding to the simulated discharges for the two scenarios are very similar. The differences observed for very low flows are difficult to be explained as the magnitude of these flows is similar to possible flow measurements errors and is also lower than the model accuracy. The FDCs obtained from the analysis of the observed discharges before and after the fire incident indicate that discharges are generally higher in the post-fire period. However, this fact cannot be attributed to the effects of fire alone as the post-fire period was wetter than the pre-fire period (

Table 4. Summary of the water balance components at an annual time step based on the model results for the two examined scenarios. PREC: precipitation, REF ET: reference evapotranspiration, PERC: percolation, ACT ET: actual evapotranspiration, SED YIELD: sediment yield.

			Conditions without the Effect of Forest Fire				Conditions with the Effect of Forest Fire
Year	PREC	REF ET	PERC	ACT ET	WATER YIELD	SED YIELD	PERC	ACT ET	WATER YIELD	SED YIELD
	mm	mm	mm	mm	mm	t/ha	mm	mm	mm	t/ha
2005	778.6	767.9	323.4	351.0	137.8	0.02	360.2	283.5	201.0	0.03
2006	692.5	824.1	263.7	372.3	111.0	0.02	319.4	259.0	149.6	0.03
2007	541.5	833.9	103.3	399.2	57.8	0.02	180.3	267.6	95.5	0.02
2008	484.0	883.1	76.4	354.4	36.6	0.01	143.0	231.2	75.5	0.02
2009	530.1	792.3	167.7	366.7	54.7	0.02	240.2	262.5	90.5	0.08
2010	669.0	872.6	155.6	402.5	62.9	0.02	249.7	290.5	111.0	0.02
2011	758.5	838.9	262.7	429.0	120.1	0.03	311.7	301.1	163.2	0.45
2012	796.0	878.2	124.8	374.4	58.1	0.02	266.2	249.2	154.8	0.05
2013	687.1	881.6	264.7	384.2	129.0	0.02	374.2	267.1	168.2	0.17
2014	744.3	828.6	192.9	414.2	62.2	0.02	298.6	298.8	150.3	0.03
Per Year	668.2	840.1	193.5	384.8	83.0	0.02	274.4	271.0	136.0	0.09

).

In Figure 5, the monthly runoff depth temporal distribution against the monthly precipitation depth is presented. As can be seen the simulated runoff depth with the pre-fire model configuration follows the observed one during the pre-fire period (until Aug 2009), but it is lower than that during the post-fire period (from Sep 2009). The opposite behavior is observed for the simulated runoff depth with the post-fire model configuration. This observation indicates that the two versions of the model describe sufficiently well the general pattern of the runoff response for the pre-fire and the post-fire conditions correspondingly. Compering the simulated monthly runoff depths for the two scenarios it can be observed that during the wet season runoff depths under the pre-fire conditions are obviously higher than the corresponding runoff depths under the post-fire conditions. During the dry season, simulated runoff depths are similar in both scenarios.

The use of SWAT in this study had the additional aim to provide information about wildfires effect on the hydrological processes and to quantify its effect on the various hydrological balance components. A summary of the water balance at an annual time step based on the model results for the pre-fire and the post-fire periods is presented in Table 4 and a graphical representation of the average annual values of the simulated water balance components is provided in Figure 6.

As can be seen in Figure 6, while the lateral flow is similar in both scenarios, surface flow is much higher in the post-fire conditions (almost three times higher). The groundwater flow is negligible in both cases. Accordingly, the higher water yield after the effect of fire is the result of the higher surface runoff during storm events. This fact explains the higher peak discharges during the post-fire period and the small differences in the low flows.

An additional interesting finding is the noticeably increased sediment yield (4.5 times higher) under the post fire conditions, especially during the wet years (Table 4), which is related to the increased surface flow. It should be noted though that the parameters related to soil erosion (soil erodibility factor) were not calibrated because of the lack of measured sediment data.

In Figure 6, the actual evapotranspiration is profoundly lower by about 29.5% in the post-fire conditions than in the pre-fire conditions. Under natural conditions, the watershed is mostly covered with bushes and trees typical to Mediterranean environments (Figure 7a) that have very deep roots and are able to extract water from the deeper soil layers or even the sallow aquifer formed in the tectonically fractured and eroded upper layers of the parent material [71]. The forest fire destroys the vegetation cover (Figure 7b) having as a result the reduction of actual evapotranspiration. The first years after the fire the watershed is covered by herbaceous vegetation during the wet period (Figure 7c); however, this type of vegetation cannot extract water from the deeper soil layers and they die during the dry period limiting the actual evapotranspiration. The lower actual evapotranspiration results in higher water yield by about 38.9% and higher percolation by about 29.4% (Figure 6 and Table 4).

Similar observations can be made in the temporal variation of the main hydrological balance components for the two examined scenarios that are illustrated in Figure 8 and Figure 9. A remarkable observation is that after the effect of the forest fire, runoff response is dominated by surface runoff mostly during wet years (Figure 8b). In dry years lateral flow is generally higher than surface flow as is also the case during the entire study period for the pre-fire conditions (Figure 8a).

Groundwater flow is negligible in all cases as the watershed is small with very sharp relief. The only exception is the last year of the study period for the post-fire conditions (Figure 8b). In this instance, noticeable groundwater flow is observed after a series of wet years and due to the decreased actual evapotranspiration. However, it should be considered that the parameters related with groundwater flow were not sensitive in the pre-fire period as the groundwater flow was negligible at that period, limiting thus the confidence in the resulted groundwater flow estimations for the post-fire period.

The temporal variation of the main hydrological balance components for the pre-fire and the post-fire conditions is illustrated in Figure 9. Interestingly, even if actual evapotranspiration is generally higher in the pre-fire scenario, the difference seems to be higher in dry years than in wet years (Figure 9a). This observation supports the assumption mentioned above that decreased evapotranspiration in the post-fire conditions could be the result of the replacement of xerophytic bushes and trees with herbaceous vegetation that is unable to extract water from deeper layers in dry years.

As it was expected, percolation is higher in wet years and especially in successive wet years. Moreover, the observed differences in percolation depths between the two scenarios also seem to be higher in dry years than in wet ones (Figure 9b) following the corresponding differences of evapotranspiration depths. A final observation is that the runoff hydrograph is sharper after the effect of the forest fire especially in wet years (Figure 9c) when surface runoff dominates the runoff response (Figure 8b). As a result, the simulated sediment yield is much higher in years characterized by surface runoff dominance (Table 4 and Figure 8b, hydrological years 2010–2011 and 2012–2013).

4. Discussion

The initial model performance indicated that, at least in this application, SWAT model is unable to describe the hydrological response of the watershed for both scenarios without calibration even if detailed soil, land covet, topographic and field data were used for its set-up. The lack of clear physical basis in models’ parameterization (i.e., need for calibration) is a major problem in ungauged basins [17,72,73,74,75], considering also that hydrological data are in many cases scarce. As also reported by Merheb et al. [75] watersheds show very heterogeneous responses over time and space, resulting in limitations in hydrological modelling and large uncertainties in predictions. Hydrological modeling is routinely used for the investigation of land use land cover changes and forest fires impact on hydrological response [12,13,14,15,18,34,35,36,37,38,39,40,41,42,43,44,45]. However, in most cases, because of the lack of observed hydrological data, especially for the post-fire period, the models are used without calibration and validation or with limited calibration. Accordingly, detailed hydro-meteorological data for an adequate period before and after a wildfire incident in a well-studied watershed are valuable for the investigation of the effect of forest fires. Though, similar datasets are rare as, in addition to the huge effort needed to secure the continuation of the measurements after the wildfire, the accidental occurrence of a wildfire is also required as it was the case in this study. For this reason, some researchers also exploit prescribed burning in similar studies [6,25,46,76].

SWAT model performed adequately well, especially in the post-fire conditions. Interestingly the calibration of the soil hydraulic properties parameters independently for each soil type achieved superior performance than their calibration in lumped form. The studied watershed is characterized by very high spatial variability [17] and for this reason, a fine discretization was used. As an example, Soulis and Valiantzas [51] suggested that the use of the SCS-CN method in lumped form may reduce the performance of the method in watersheds characterized by increased spatial variability. In this study, the spatially distributed CN parameter values that were directly estimated using measured rainfall and runoff data by Soulis [33] for the studied watershed were successfully used. In this manner, the number of calibration parameters was reduced. Further calibration of CN parameter did not provide any improvements in the model performance. Accordingly, using a method for the direct estimation of CN parameters, such as the method proposed by Soulis and Valiantzas [33,51,52], may provide a good alternative facilitating the calibration of hydrological models.

Temporal variability of weather conditions seems to create additional challenges in the model calibration. In particular, superior performance was achieved in wet years than in dry years in both model configurations. This observation may indicate that hydrological functioning may vary between wet and dry years as a result of corresponding changes in vegetation cover and soil properties. However, it should be also considered that high flows have a greater influence in the calibration objective functions than low flows. Tufekcioglu et al. [45] in their study, attributed the unsuccessful calibration of SWAT model, to the spatial and temporal variability, to the small size of the watershed, and the very high percolation. These conditions are very similar to the characteristics of the studied watershed in this study. Similar challenges are also reported by other researchers [12,77,78,79].

An additional challenge is that even after successful calibration of SWAT model, there is still a reservation on the ability of the model to accurately represent hydrological processes and quantify hydrological balance components as is typically the case for complex models involving numerous parameters [74]. The obtained results revealed that different sets of calibration parameters resulted in similar model performances revealing the increased uncertainty involved in the estimation of the parameters of the models. For this reason, in this study, it was decided to use the same basic model configuration for both scenarios and to separately calibrate only the parameters affected by the forest fire such as the soil hydraulic properties and evapotranspiration-related coefficients for the burned areas. Future research could focus on the utilization of additional data for the model calibration such as field or remotely sensed soil moisture variability observations.

Despite the above limitations, SWAT model is proven to be a valuable tool for distinguishing the effect of forest fires from the effect of weather conditions variability and for investigating the wildfire effect on the hydrological balance components. SWAT model has been also applied successfully in many relevant studies [12,14,56,57,58,59,60,61,62,63,80,81,82].

The results of the current analysis highlight the profound effect of the 2009 fire incident on the hydrological processes of the studied watershed. In most studies, the loss of vegetation cover is considered as a key factor for forest fires impact on hydrological response [6,7,8,9,10,11,12,13,14,15,18]. The post-fire alteration of soil hydraulic properties, and especially of soil hydraulic conductivity, is an additional important factor leading to reduced infiltration and increased overland flow [1,16,17,18,44]. However, it should be noted that reduced hydraulic conductivities are crucial at the short term but they are observed only at the upper few centimeters of the soil profile and only for a narrow period of time [1,17,44]. Soulis et al. [17] reported that the observed soil hydraulic conductivity values in the studied watershed were significantly reduced after the effect of the forest fire but seven years after the wildfire they were similar to the pre-fire period. Accordingly, the lasting effects of forest fires that are observed many years after the fire cannot be attributed to soil hydraulic properties alone. For example, Zhou et al. [18] found that for the first 15 years after a bushfire in southeast Australia there was a substantial increase in streamflow; this was attributed to initial decreases in evapotranspiration and soil infiltration rates resulting from the fires, followed by logging activity. Pereira et al. [13] also observed increased peak runoffs in the post-fire period and lower runoff during the dry periods in comparison to the pre-fire period. Other researchers also reported profound reductions in evapotranspiration and increases in water yield [6,7,8,9,10,11,18].

The obtained results agree with the previous studies showing a profound reduction of the actual evapotranspiration (29.5%) and increase in water yield (38.9%) and percolation (29.4%) in the post-fire conditions in comparison to the pre-fire conditions. Increase of groundwater recharge is also reported in previous studies [19,20,21]. The observed reduction of the actual evapotranspiration seems to have a dominant role in the alteration of the water balance. As can be seen in Figure 6, surface runoff is significantly increased due to the reduction of infiltration as a result of the wildfire; however, this increase is much smaller than the reduction of actual evapotranspiration. Nevertheless, even moderate increases of surface runoff can result in considerable increases in peak flows as they consider short periods of time during the storm events. As it was explained in Section 3, the decreased evapotranspiration in the post-fire conditions could be related to the destruction of the typical vegetation consisting of xerophytic bushes and trees with deep roots and adaptations allowing them to extract water from deep soil layers and sallow aquifers and its replacement by grasses and herbaceous vegetation that have not this ability. Poon and Kinoshita [15] draw similar conclusions; however, more research at smaller scales is required to investigate this assumption.

The post-fire vegetation recovery in the studied watershed was initiated soon after the fire. During the first autumn after the fire, the bigger part of the watershed was already covered with grasses, and some oak shrubs were started to regenerate. The following years, xerophytic bushes such as thymus, organum and cistus started to glow and progressively replacing grasses. However, the pine forest started to recover very slowly and only in small parts of the watershed. The gradual recovery of the vegetation denotes that forest fire impacts are expected to be time varying as suggested by Bart [83].

The runoff hydrograph is also much sharper with considerably higher peak discharges after the effect of the forest fire especially in wet years, when runoff response is dominated by surface runoff. This is in accordance with the findings of most relevant studies which attribute this behavior to lowered evapotranspiration, lowered infiltration, and increased and more rapid overland flow [1,16,18,44]. The comparison of the hydrological balance simulated for the two scenarios leads to a similar conclusion, i.e., that at least in larger time steps and in the long term the main driver is the decrease of evapotranspiration. The vast increase of peak discharges, however, cannot be attributed solely to the decrease of evapotranspiration, but it may be the combined result of many factors such as reduced soil hydraulic conductivities, fast surface and channel flow due to vegetation destruction, and increased initial soil moisture due to evapotranspiration reduction. Other disturbances such as grazing, soil disturbance during reforestation campaigns, forest logging, and other could also be involved [17,33]. The identification of the mechanisms behind the increased runoff and peak discharges is really challenging and further detailed research is required as it may provide a better insight on the hydrological response recovery process and the length of the post-fire recovery period as well as on the effectiveness of the various post-fire mitigation measures. The estimation of the length of the recovery period is also important as it defines the window of opportunity for the occurrence of a big storm in the burnt watershed with possibly catastrophic consequences. For example, in the studied watershed such a storm event occurred on 30 December 2012 resulting in the highest recorded daily discharge (4.94 m³/s, 7.8 times higher than the pre-fire period) and the highest recorded peak discharge (81.25 m³/s, an order higher than the pre-fire period). This event caused many damages in roads, bridges, and other infrastructures in the study region. The same storm could be less damaging in case that it has occurred before the fire incident or after the recovery period. In contrast, even higher peak discharges and more catastrophic damages could appear if a similar storm has occurred in the first winter after the fire. Candela et al. [34], in their study investigating the effects of forest fires on flood frequency curves in a Mediterranean watershed, concluded that flood frequency is greatly increased during the post-fire period. These results are in line with the results obtained in the present study (Figure 4). Nalbantis and Lymperopoulos [36] also investigated methods to assess flood frequency after forest fires in ungauged basins.

Additional research is also required to investigate the combined effects of ecological disturbances such as wildfires and climate change. Previous studies by Wine and Cadol [21] and Wine et al. [30] demonstrated that “at small spatial scales wildfire contribution to stream-flow is intermittent, whereas if a threshold of aerial extent continues to be exceeded at larger spatial scales this patchwork of wildfire histories integrates to sustained increases in streamflow due to wildfire effects that may rival or exceed projected climate change impacts.”

A final finding is the noticeable increase of sediment yield (4.5 times higher) under the post fire conditions. Increased erodibility and sediment transport can be the combined result of vegetation destruction and soil structure destabilization and increased surface flow. Several researchers reported considerable increase of sediment yield as well [7,8,11,12,25,26,27,29,43,47].

5. Conclusions

In this study, the impacts of forest fires on the hydrological response of natural Mediterranean watersheds was investigated using as an example a small experimental watershed. Detailed hydrometeorological data recorded for the adequate period before and after the fire incidence were used in conjunction with SWAT model to separate the effects of the wildfire from the effects of weather variability. The main findings of this study are the following:

• SWAT model was unable to describe the hydrological response of the watershed without extensive calibration for both the pre-fire and the post-fire conditions even if detailed soil, land covet, topographic, and field data were used for its set-up.
• Detailed hydro-meteorological data for an adequate period before and after a wildfire incident in a well-studied watershed are valuable for the investigation of the effect of forest fires.
• Despite the serious limitations, SWAT model was proven to be a valuable tool for the investigation of the forest fire effects on the hydrological response given that sufficient hydrometeorological data for the model calibration are available.
• A profound reduction of the actual evapotranspiration (29.5%) and increase in water yield (38.9%) and percolation (29.4%) in the post-fire conditions were observed.
• Runoff hydrograph is much sharper with considerably higher peak discharges after the effect of the forest fire especially in wet years when runoff response is dominated by surface runoff.
• The main drivers behind the alterations in hydrological response seem to be that of the vegetation cover change and the reduction of soil hydraulic conductivity; however, the involved mechanism are not clear and other disturbances should be also considered.
• A noticeable increase of sediment yield under the post fire conditions was observed.

Finally, even if the above findings provide additional information that may help to improve the current understanding on the effect of forest fires on hydrological response, the obtained results indicated that more research is still required in many aspects of this study. For instance, future research may focus on improving SWAT model parametrization and calibration considering the spatial and temporal variability in the watershed characteristics and weather conditions making use of additional data such as the constantly improving remotely sensed soil moisture spatial distribution. Additional research is also required to investigate the combined effects of ecological disturbances such as wildfires and climate change. The mechanisms behind the increased surface flow, peak discharges and sediment yield, and the post-fire recovery processes should be investigated in greater detail as well.