Event-based quickflow simulation with OpenLISEM in a burned Mediterranean forest catchment

Keywords: autocalibration, catchment scale, eucalypt, event-based modelling, maritime pine, post-fire, rainfall-runoff modelling, soil moisture content, soil water repellency, surface runoff.

During the last decades, wildfire occurrence has increased in the Iberian Peninsula, and the burnt area has, on average, exceeded 100 000 ha every year (San-Miguel-Ayanz et al. 2013). Wildfire frequency is expected to further increase in the future, owing to changes in climate and socio-economic drivers (Turco et al. 2014, 2016; Viedma et al. 2015; Calheiros et al. 2021).

The direct impacts of wildfires typically include the consumption of the vegetation and litter layer as well as heating-induced changes in topsoil properties (Shakesby and Doerr 2006; Moody et al. 2013; Ebel 2020). Thereby, wildfires decrease interception and increase both effective rainfall and rainfall erosivity (Cerdá and Doerr 2005; Stoof et al. 2012). Among the topsoil properties affected by fire, soil water repellency (SWR) has received much research attention because of its expected role in post-fire runoff generation (Doerr et al. 1996; DeBano 2000; Ferreira et al. 2000). SWR expresses the resistance of soils against wetting for periods ranging from seconds to weeks (DeBano 2000; Malvar et al. 2016; Mao et al. 2019). SWR is a highly dynamic property, varying between dry and wet seasons in both unburnt and burnt soils (Doerr et al. 2000; Santos et al. 2016; Martins et al. 2020). SWR severity has often been related to soil moisture content (SMC) but the fairly few high-resolution time series that exist to date have revealed complex relationships in burnt areas (Keizer et al. 2008; Malvar et al. 2016; Robichaud et al. 2016). This complexity has been explained by the existence of a so-called transition zone in SMC where soils can be both repellent and wettable (Crockford et al. 1991; Dekker et al. 2001). Because of SWR’s dynamic nature, and the lack of non-destructive field measurement methods, its contribution to post-fire runoff generation has been hard to quantify (Shakesby and Doerr 2006; Shakesby 2011).

The role of SWR in post-fire runoff generation has also been poorly explored using rainfall runoff models (Lopes et al. 2021). Post-fire hydrological modelling has focused on coarse temporal resolutions, with the bulk of the studies using a daily time-step (Dun et al. 2009; Seibert et al. 2010; Mahat et al. 2015; Moussoulis et al. 2015), and very few studies using a (sub-)hourly resolution (Canfield et al. 2005; Seibert et al. 2010; Van Eck et al. 2016; Rengers et al. 2016; Wu et al. 2021a, 2021b). Overall, efforts to assess model suitability for post-fire conditions have suffered from the lack of field data, as only a limited number of experimental catchments have been instrumented and monitored during the full window of disturbance period (Shakesby 2011; Moody et al. 2013). To date, only a handful of post-fire hydrological modelling studies have explicitly included the role of SWR. Liu et al. (2021) correlated temporal increases in saturated hydraulic conductivity with decreases in SWR due to the wet pre-event conditions. Vieira et al. (2014) simulated the role of SWR in the seasonal runoff response by reducing the effective infiltration depth in the Morgan–Morgan–Finney model with increasing SWR severity. Nunes et al. (2016) used the Thornthwaite–Mather model to simulate the role of SWR in post-fire infiltration through changes in daily soil moisture content, decreasing field capacity with increasing SWR severity. McGuire et al. (2018) assessed the use of spatial saturated hydraulic conductivity at catchment scale to simulate post‐wildfire hydrological response and debris flow initiation under repellent conditions, finding the influence of sorptivity on infiltration rates to be negligible. Van Eck et al. (2016) is the only study at catchment scale that employed a process-based model to explore how variation in SWR affected the runoff response.

The present study is a follow-up of Van Eck et al. (2016) to improve the understanding of the hydrological response of recently burnt catchments and, in particular, the role therein of SWR, through the application of a process-based model (OpenLISEM). The role of SWR is explored by calibrating infiltration capacity for events with dry and repellent as well as wet and non-repellent antecedent conditions. To this end, an experimental catchment was selected that was small and entirely burnt, where SMC and SWR had been measured (Vieira et al. 2016, 2018), and for which the model input data requirements of OpenLISEM were met (see Van Eck et al. 2016).

The specific objectives were to: (i) quantify how well OpenLISEM could simulate the hydrological response at catchment scale for events with contrasting initial conditions of SMC and SWR; (ii) optimise OpenLISEM performance through a generalised calibration approach, adjusting maximum infiltration capacity as a function of initial SMC; (iii) assess how well within-catchment spatial runoff predictions agreed with 1- to 2-weekly field observations, in particular of overland flow at the microplot scale.

On 27 August 2008, a wildfire consumed almost 68 ha of forest near the village of Colmeal, in the municipality of Góis, north-central Portugal (40°08′42″N, 7°59′16″W; 490 m above sea level (asl)). This burnt area included a small catchment of 10 ha, dominated prior to the fire by maritime pine (Pinus pinaster Ait.) stands and eucalypt (Eucalyptus globulus Labill.). Land management in this catchment has been extensively described in Vieira et al. (2016) and can be summarised as a combination of three forest types (pine and shrubs, eucalypt and shrubs, and eucalypt plantations), and four types of soil mobilisation operations (none, contour ploughing, downslope ploughing and terracing; Fig. 1). As described in Vieira et al. (2016, 2018), vegetation and hydrological recovery in the catchment after the 2008 fire were limited as a consequence of past disturbances, especially previous wildfires and the soil mobilisation operations mentioned above.

Fig. 1.  Colmeal catchment: (a) land use; (b) land operations; and (c) elevation, soil moisture sensor locations (sm 1, sm 2), outlet, roads and waterlines.

The climate of the study area can be characterised as humid mesothermal (Köppen Csb, Peel et al. 2007), with prolonged dry and warm summers. The mean annual temperature and precipitation at the nearest meteorological station (Góis (station 13I/01G); 10 km) are 12°C and 1133 mm, respectively (SNIRH 2011). The occurrence of precipitation in this area is most frequent between November and February, whereas rainfall events between March and May are less frequent but can have high intensities (Vieira et al. 2018).

The parent material in the study catchment consists of pre-Ordovician schists and greywackes (De Brum Ferreira 1978; Pimentel 1994), which have resulted in shallow soils that are typically mapped as Humic Cambisols (Cardoso et al. 1971). The A horizons in the study catchment have a coarse sandy loam texture (sand > 70%) and a high stone content (40–46%).

According to field indicators (i.e. canopy and woody and litter consumption, ash colour and mineral soil), the burn severity was low to moderate for vegetation and soil, as tree canopies and logs were only partially consumed, the litter layer was fully consumed, the ash was black and the mineral soil was unaffected (Hungerford 1996; DeBano et al. 1998).

After the fire, the catchment was instrumented with a hydraulic channel (H flume) that was equipped with an ultrasonic water level sensor and a rain gauge (0.2 mm tips) connected to a Campbell data logger (CR1000) to measure rainfall and streamflow with 10-min resolution during 4 post-fire years. Additional field measurements of SMC, surface runoff, ground cover and SWR were taken at several points in the burned area. SMC was registered continuously at 3–5 cm soil depth (using Decagon EC-5 soil moisture sensors at three locations, one close to a eucalypt plantation (sm1), another near the pines (sm2), and one located at the outlet (Fig. 1c). Surface runoff was assessed in 12 microplots (0.25–0.5 m²), and monitored with a weekly frequency depending on rainfall (Vieira et al. 2016, 2018). At those same microplots, ground cover was also assessed monthly, whereas soil texture, soil depth and soil roughness were characterised once (Vieira et al. 2016, 2018). SWR was measured in the field at monthly intervals in a eucalypt plantation and a pine stand. The SWR monitoring sites were selected immediately downstream from the experimental catchment because of the extensive soil disturbance caused by the measurements, but both sites were comparable with the eucalypt and pine stands within the catchment, including vegetation and soil burnt severity. The monthly SWR data were used to identify the periods when soils were likely to be repellent (Supplementary Material S1). Finally, soil texture, soil depth and soil surface roughness were measured once in all microplots.

The hydrological processes in the catchment were simulated using the Limburg Soil Erosion Model version 2.01 (De Roo et al. 1996). OpenLISEM is a process-based and spatially distributed model that simulates rainfall runoff events using a geographic information system (De Roo et al. 1996). The processes incorporated in OpenLISEM are precipitation, interception, surface storage, infiltration and water flux into the soil, surface runoff, and channel flow (De Roo et al. 1996). Optionally, erosion processes can also be simulated, which was not performed in the present study because it focused on the hydrological processes in the catchment. In this study, infiltration was simulated with the Green–Ampt one-layer method, while overland flow was routed using the steepest descent flow method.

The rainfall input data for OpenLISEM was obtained by summing the registered tips (of 0.2 mm) of a rain gauge located at the outlet over 10 min intervals and this was used as rainfall intensity (mm h⁻¹). In addition, OpenLISEM input maps with a 5 m spatial resolution were created with respect to: location of the catchment outlet; terrain (digital elevation model (DEM)); land cover characteristics (ground cover (–), leaf area index (m² m^–2), plant height (m), width of impermeable roads (m)); land surface characteristics (random roughness (cm), Manning’s n (–), stone fraction (–)); soil physical characteristics (cohesion of bare soil (kPa), cohesion by roots (kPa), aggregate stability (–), median soil texture (µm), soil depth (mm)); soil hydrological characteristics (saturated hydraulic conductivity (mm h⁻¹), soil water tension at the wetting front (cm), saturated soil moisture content (cm³ cm⁻³), initial SMC (cm³ cm⁻³)). The required input data maps are spatially variable, based on differences in land use and management (Table 1). Cohesion by roots, aggregate stability and plant height were estimated to be zero as a consequence of the wildfire. Cohesion of bare soil was set at 20 kPa, while median soil texture was estimated to be 80 µm. The width of impermeable roads was set to 3 m for all grid cells that included the forest tracks in the catchment (Fig. 1). Initial SMC (θ_i, cm³ cm⁻³) varied according to what was recorded by the field soil moisture sensors (sm1, sm2 and outlet). Input maps were the same for each event, except for the calibrated parameters, which varied over time and between events (see Event selection criteria and modelling aproach). The spatial resolution of all maps was 5 m and a fixed model time step of 5 s was used.

Table 1.  Model inputs according to land use and land management characteristics.

The outputs used from OpenLISEM simulations for each event were the time series of the simulated quickflow at the outlet, and the spatially explicit map of the surface runoff generated in the catchment. Following Reitz and Sanford (2019), quickflow was defined as the rapidly varying runoff portion of the discharge hydrograph, with the total discharge being the sum of quickflow and baseflow. This is illustrated in Fig. 2. Quick flow was derived from the observed streamflow data using the baseflow separation technique proposed by Arnold et al. (1995).

Fig. 2.  Schematic representation of hydrological event with a single peak. Discharge flow as light blue line, overlapped by baseflow in dark blue area, and quickflow (discharge minus the baseflow) in grey area. Three panels identify antecedent conditions (left), event (middle), and post-event periods (right), as well as their criteria for event separation and indicators used. The left panel illustrates the time frame in which baseflow and initial soil moisture content (θ_i) were determined.

A total of 16 rainfall runoff events were selected for this modelling exercise (Table 2). All occurred during the first 2 post-fire years, as overland flow generation was most pronounced during this initial phase of the window of disturbance (Vieira et al. 2018), and were separated by at least 1 h without rainfall. The 16 events were selected because they had corresponding major quickflow peaks of more than 8 L s⁻¹, and contrasting initial SMC conditions (Fig. 2). The events were classified into two classes of initial SMC conditions, with ‘dry’ or ‘wet’ events having mean SMCs smaller or larger than 17% (Table 2). This threshold value corresponded to the mean SMC conditions over the entire study area. This resulted in seven dry and nine wet events (Table 2). The dry and wet events involved similar median rainfall amounts but with a greater variability in the case of the wet events (Supplementary Fig. S2d). The wet events involved consistently higher median and maximum amounts of total quickflow, higher peak quickflows and higher quickflow coefficients than the dry events (Supplementary Fig. S2). Further details on the differences between the wet and dry events are given in Supplementary Material S3.

Table 2.  List and characteristics of all the simulated events.

The modelling procedure included calibration, validation and uncertainty analysis (Fig. 3; see section 2.4). Model calibration, focused on the outlet quickflow, was performed by adjusting three OpenLISEM parameters for each individual event, i.e. Manning’s n (n, –), saturated soil moisture content (θ_s, cm³ cm⁻³), and saturated hydraulic conductivity (Ksat, mm h⁻¹) (Table 1). The calibration procedure involved multiplying the respective input maps by a multiplication factor. A similar calibration was applied by Vieira et al. (2014) and Nunes et al. (2016) for burned areas; however, these prior studies involved seasonal and daily time steps as opposed to the sub-hourly time step used here.

Fig. 3.  Modelling approach scheme.

The events were randomly separated into a calibration and a validation dataset, by alternately selecting events from wet and dry datasets (Table 2). An autocalibration was then performed and an optimal parameter set was obtained for the wet and dry calibration events respectively (Fig. 3, Table 2). To assess the feasibility of a calibration validation split test for distinct moisture classes, the optimal parameters were then applied for the validation events (Fig. 3, Table 2), and the resulting model performance was compared with an independent autocalibration for the validation events alone. The split test is frequently used in hydrological modelling studies to provide a test of the statistical significance of the model performance metrics (Liu et al. 2018).

Besides the simulation of repellent conditions, other improvements were implemented in the current study following the simulations of Van Eck et al. (2016), such as the inclusion of forest roads, the resolution increase (10 vs 5 m grid size), which led to the decrease of the simulation time step (5 s), the selection of more events to be simulated (5 vs 16), the inclusion of continuous soil moisture sensor data (Fig. 1) and the uncertainty analysis. Nevertheless, other methodological options taken by Van Eck et al. (2016) were kept, such as the land use and the land operations delineation areas as well as most of the physical soil parameter values (Table 1).

In addition to the model calibration and validation focused on quickflow at the outlet, we also validated the spatial maps of runoff generated by OpenLISEM, by comparing them with runoff measurements obtained from microplots within the catchment (0.25–0.5 m²).

Model performance for quickflow simulations was evaluated for: (i) total quickflow (m³); (ii) peak flow (L s⁻¹) until a maximum of four peaks per event; and (iii) timing of the peaks (min) until a maximum of four peaks per event. Three commonly used statistical indicators (Moriasi et al. 2015) were calculated for each of the three model outputs (i–iii) to assess model performance:

1. Nash–Sutcliffe efficiency (NSE) – determines the relative magnitude of the residual variance compared with the measured variance. NSE values larger than 0.5 indicate satisfactory model performance and values > 0.8 indicate very good model performance (Moriasi et al. 2015)

2. Coefficient of determination (R²) – describes the proportion of data variance explained by the model. R² ranges from 0 to 1, with values higher than 0.5 indicating reasonable model performance (Santhi et al. 2001; Van Liew et al. 2003).

3. Percentage bias (PBIAS) – indicates the magnitude of model errors compared with measurements. Positive PBIAS values indicate model underestimation and negative values model overestimation (Vijai et al. 1999). PBIAS values smaller than ±15% for runoff are considered satisfactory (Moriasi et al. 2015).

It should be noted that these metrics and thresholds referenced in Moriasi et al. (2015) are indicative for daily time-steps, whereas in the current study, model performance was evaluated for hydrographs at the catchment outlet using measured data with 10 min temporal resolution, and by comparing them with the modelled hydrograph.

Autocalibration for the split-test (calibration and validation, Fig. 3, Table 2) was performed executing a constrained fitting of the model to quickflow with R. This fit was executed with a model fitting function (modFit), under the Pseudo-random Search Optimisation Algorithm (Price 1977). Parametric uncertainty for each event was obtained by applying a Markov Chain Monte Carlo (MCMC) function in R (modMCMC). Both model fitting and MCMC procedure consisted of 1000 iterations of the chosen parameters under the lower and upper provided bounds, executed under the Flexible Modelling Environment (FME) package in R (Soetaert and Petzoldt 2010), which has been previously applied for hydrological modelling with Soil and Water Assessment Tool (SWAT) (Wu and Liu 2012).

Finally, three statistical models were built to predict Manning’s n, Ksat and θ_s: one model for all the events (overall), and two others considering wet and dry conditions separately. This was intended to relate the autocalibration results to several additional explanatory variables in an attempt to understand which variables better explain the obtained model calibration. This was done in Rstudio with the linear model (Chambers and Hastie 1992) function lm {stats}, considering as response variables the calibrated Manning’s n, Ksat and θ_s, and as possible explanatory variables the Antecedent Precipitation Index from the 10 days prior to the simulated event (API, –), baseflow at the beginning of the event (L s⁻¹), event rainfall duration (min), maximum rainfall intensity (mm h⁻¹), initial soil moisture content (θ_i, cm³ cm⁻³), time since fire (days) and total event rainfall (mm). A compilation of these explanatory variables can be found in Supplementary Table S2.

Autocalibration with OpenLISEM allowed the prediction of the hydrological response of individual events at the catchment outlet with a good to very good model performance (NSE ranging from 0.74 to 0.97; Table 3, Fig. 4). When evaluating OpenLISEM performance for all calibrated events combined (overall), results showed very good model performance in predicting total quickflow, peak flows and time of the peak (Table 4). However, when splitting these events by moisture class (see Supplementary Material S3), OpenLISEM performed better for wet (NSE ranging from 0.83 to 1.00) than for dry conditions (NSE between 0.52 and 1.00), with the latter showing a satisfactory (NSE 0.52) model performance for peak flow, and very good for total quickflow (NSE 0.95) and time of the peak (NSE 1.00).

Table 3.  OpenLISEM calibration range inputs and model performance for quickflow calibration of individual events.

Fig. 4.  Measured and OpenLISEM simulated (a, b) total quickflow (n = 16); (c, d) peak flows (n = 28); and (e, f) time to peak (n = 28) for calibration (a, c, e) and validation (b, d, f) respectively. Details on modelling approach and data used for calibration and validation can be found in Fig. 3 and Table 2.

Table 4.  OpenLISEM model efficiency for calibration and validation datasets.

The validation results were not entirely satisfactory as can be seen by the decrease in model performance, especially for total quickflow and peak flow (Table 4, Fig. 4). Autocalibration results with this specific dataset also showed lower model performances in comparison with the calibration ones (Table 4), especially for total quickflow by obtaining a model performance below the satisfactory level (NSE 0.46) for the overall dataset. Despite that, the resulting decrease in model performance for validation events was still far from the autocalibrated one (Table 4), suggesting that either OpenLISEM is very sensitive for small changes between event characteristics or that the calibrated dataset was not robust enough for a split test.

The linear regression for all optimal parameter sets using the observed auxiliary variables (Supplementary Table S1) showed that wet and dry events were affected by distinct combinations of variables besides initial SMC alone. Still, splitting the data by SMC facilitates the calibration procedure as shown by the increased adjusted R² for wet and dry conditions when compared with all events combined (Supplementary Table S1).

The uncertainty analysis performed for each event through the MCMC procedure also showed that a small input variation can result in substantial changes in the quickflow predictions (Supplementary Figs S5, S7, S9, S11). In the illustrated event c (Fig. 5), quickflow prediction varied between 0 and 10 L s⁻¹ for the first quickflow peak, and 3 and 12 L s⁻¹ for the second quickflow peak, with small variations of Manning’s n (0.39 – 0.59) and Ksat (5 – 9 mm h⁻¹). Nevertheless, despite these uncertainties in predicting total quickflow and peak flows, the timing of the peaks was constantly predicted in line with the observed timing of the flow.

Fig. 5.  OpenLISEM quickflow predictions variation for event c simulation following 1000 MCMC iterations with NSE > 0.5. Note minimum and maximum quickflow predictions in light grey range, quickflow predictions standard deviation in dark grey range, best quickflow prediction as black line, and observed quickflow as blue line.

The most notable result was that very low or sometimes no surface runoff was predicted by OpenLISEM on two slopes where field microplots were installed and runoff was observed (Supplementary Figs S6, S8, S10, S12). Those plots were located at Eucalypt contour ploughed and Eucalypt downslope ploughed sites. On the remaining plots (Pine unploughed site), model simulations presented a better agreement with the field visits during wet events than when compared with the dry events, but in general OpenLISEM simulated lower runoff amounts than those measured by all field plots (Fig. 6). Additionally, the connectivity between the upper part of the catchment and the bottom part near the outlet was only verified in 63% of the simulated events (e.g. events a, b Fig. 7), which is less than was verified with field observations (75%) concerning those specific events, whereas the remaining events presented a contributing area limited to the Pine site and the road adjacent to the outlet (event i, Supplementary Fig. S10).

Fig. 6.  Rainfall for weekly field monitoring (dark blue, upper axis) vs measured runoff at microplots (orange for dry and blue for wet, bottom axis), and event rainfall (light blue, upper axis) vs simulated event runoff (grey, bottom axis) with the OpenLISEM model at the Pine unploughed site. Note field measurements nos 31, 56, 62 and 63 had more than one event in that week.

Fig. 7.  OpenLISEM spatial runoff predictions for events (a) (left), and (b) (right) for optimal parameterisation set.

The validation of the model simulations with the field data at microplot scale revealed limitations due to the different time steps used (event vs weekly). The primary reason for those limitations was the fact that the rainfall event contribution to the total hydrological response of that week was highly variable (17–97%).

In this study, we calibrated OpenLISEM through varying Manning’s n, Ksat and θ_s according to an SMC threshold for catchment dryness or wetness at the beginning of each event. This study hypothesised that such classification could potentially be used as a proxy for SWR. However, in the same way as the preceding study (Van Eck et al. 2016), the main difficulties found in the simulation of quickflow with OpenLISEM were associated with the occurrence of SWR and its spatial and temporal expression, as shown by the lower model performance for drier conditions (Table 3).

Our results showed that the current classification for wet or dry conditions may require some improvement, by detailing such classification for each land use or operation, or by considering an SMC gradient along the catchment. There are some approaches for SMC or soil hydrological properties based on topographic analysis that can provide a reasonable estimation to improve hydrological modelling (Nunes et al. 2009; Tavares Wahren et al. 2016).

To overcome the low resolution of SWR measurements, exploring continuous SMC data as an indicator of SWR was expected to be a promising alternative, primarily because of the well-known relationship between SWR and SMC (Doerr and Thomas 2000; Malvar et al. 2016; Robichaud et al. 2016), and secondly, because it would be available through the use of sensors, independently of the event timing. However, it seems that other variables besides SMC interfere in the hydrological response at this scale, as demonstrated by the multiple significant relationships between the calibrated inputs and other auxiliary variables such as API, time since fire, baseflow or rainfall characteristics (Supplementary Table S1). For example, the rainfall conditions in the wet events were similar between the calibration and validation events, while dry events showed greater dispersion in the auxiliary variables. This may have been one of the reasons why model performance on validation was not entirely successful. Therefore, we hypothesise that future modelling exercises at these spatial and temporal scales may require a greater number of events for calibration, thus allowing integration of larger hydrological variability.

Despite the limited results, we still consider that OpenLISEM could be improved by including an SWR module, which would limit the water infiltration capacity in the presence of SWR conditions or the presence of low SMC. This could be done through the reduction of θ_s throughout the entire event, as performed in this study, or only during a short initial stage of the rainfall event, as verified in several field observations, whereas the SWR barrier breaks at a certain rainfall amount or SMC threshold (Crockford et al. 1991; Doerr and Thomas 2000; Malvar et al. 2016; Mao et al. 2019); this approach was successfully tested by Nunes et al. (2016) for SWR breaking after a sequence of wet days.

The calibration results showed good agreement between observed and simulated total flow and peak flow at the outlet. However, PBIAS analysis showed some level of underestimation of the total flow and peak flow but still with good model performance. OpenLISEM had greater difficulties in simulating dry than wet events, whereas peak flows were underestimated for dry conditions, as shown by the satisfactory PBIAS metrics. Similar results were also found by Van Eck et al. (2016), where OpenLISEM underestimated flows under prevalent SWR conditions, but overestimated flows when SWR was less prevalent. The validation results, however, showed that the number of events used was not sufficient to provide a robust calibration with OpenLISEM regarding total flow and peak flow. Nonetheless, it seemed sufficient to accurately predict the peak timing, which is an excellent result for predicting extreme events and helps to improve preparedness in facing possible hydrological and erosive risks downstream.

In terms of spatial simulation of runoff, OpenLISEM revealed poor agreement with field measurements, by simulating limited or no runoff in two out of three monitored slopes (Figs 6, 7). This also happened in other studies with this same model (Jetten et al. 2003), with this problem being attributed to an inadequate infiltration estimation or problems in the flow network, which was not identified in the current study. This disagreement occurred for two slopes (eucalypt) that had been subjected to land operations (contour and downslope ploughing), and this difference between simulations and reality came from the elevated infiltration simulated in these locations, as a consequence of the higher surface roughness and the greater soil depth. Additionally, OpenLISEM also presented difficulties in addressing catchment connectivity, especially in wet conditions, when baseflow was larger. The fact that OpenLISEM cannot simulate baseflow, even when an accurate SMC at the beginning of the event is provided, may lead to a limited connection between slopes, especially near the streams where soils are likely nearly saturated in wet conditions (Supplementary Fig. S8).

The methodological approach used in this study to address post-fire conditions and SWR through the calibration of events according to their moisture class allowed good overall model performance to be obtained but also revealed several difficulties in adjusting infiltration through θ_s. A similar approach was taken by Moussoulis et al. (2015) who assessed the effect of SWR in their runoff predictions by calibrating the threshold depth for infiltration with Système Hydrologique Européen (MIKE SHE) (Abbott 1986) in seven partially burned catchments in Greece. Moussoulis et al. (2015) also calibrated Manning’s n for the peak flow size and timing, obtaining a satisfactory to good model performance at a monthly scale, whereas at a daily time-step, MIKE SHE was not satisfactory. Without addressing SWR specifically but only addressing post-fire conditions, Seibert et al. (2010) suggested the best calibration to address fire-induced changes would be by reducing water storage capacity within the Hydrologiska Byråns Vattenbalansavdelning (HBV) model (Bergstrom 1992), which is in line with the proposed changes in the present study. Alternatively, Canfield et al. (2005), Goodrich et al. (2005), Rengers et al. (2016, 2019) and Wu et al. (2021a, 2021b) suggested that post-fire hydrological calibrations should be focused on Ksat and Manning’s n, which would lead to higher runoff peaks without a dramatic increase in runoff volume. The calibration of Canfield et al. (2005) resulted in good to very good NSE performances for individual events; however, events resulting in an NSE below 0.7 were removed from their analysis. Wu et al. (2021a, 2021b) obtained satisfactory to very good model performance with automatic calibration for quickflow, whereas Rengers et al. (2016, 2019) assessed model performance focusing on peak flow timing for both calibration and validation.

The fairly high post-fire Manning’s n estimations obtained after the autocalibration procedure (Table 3) agree with the estimations obtained by Wu et al. (2021a, 2021b) for burned and unburned conditions for the south of Portugal. They can be explained by the deep soil interventions found in these locations such as contour ploughing and terracing (Supplementary Fig. S1), which are frequently found in the Mediterranean landscape (Martins et al. 2013), and also supported by the roughness measurements obtained in situ (Vieira et al. 2016).

For Ksat parameterisation, Ebel and Martin (2017) and Ebel (2020) estimated that Ksat for burned areas should decrease substantially after fire, as a consequence of changes in soil organic matter and soil structure, followed by an increase to background levels during the 4–5 post-fire years. Although the authors could not attribute this recovery to a specific cause, they suggest similarities with mechanisms documented in other studies (e.g. Robichaud et al. 2016), where Ksat recovery is driven by vegetation and ground cover recovery. However, during the 2 post-fire years of this study, we only found a significant inverse relationship between Ksat and time since fire (estimate −0.19) for dry events (P-value < 0.001, Adjusted R² 0.95, Supplementary Table S1). Wu et al. (2021a, 2021b), using the OpenLISEM model in a larger catchment in southern Portugal, did not find a decrease in calibrated Ksat values for post-fire versus pre-fire events.

The application of OpenLISEM to the Colmeal catchment allowed us to accurately predict post-fire hydrological response in a small catchment at the event scale. The approach adopted in this study agrees with the methodologies applied in the limited amount of studies that address post-fire hydrological modelling at catchment scale (Canfield et al. 2005; Goodrich et al. 2005; Van Eck et al. 2016; Rengers et al. 2016; Wu et al. 2021a, 2021b), but may still require additional improvement in the classification of SMC and its spatial variability. Nevertheless, our results showed that OpenLISEM can be a promising tool in determining the time of the peak, while less efficiently determining total quickflow and peak flows. Moreover, these predictions were achieved in great detail, whereas most studies addressing this matter use daily to monthly time-steps (e.g. Dun et al. 2009; Mahat et al. 2015; Seibert et al. 2010), and often involve the study of partially burned catchments (e.g. Moussoulis et al. 2015; Nunes et al. 2018; Basso et al. 2020). In this sense, we believe this study can provide new insights about post-fire hydrological modelling by studying the hydrological response from several individual rainfall events in an entirely burned catchment. Although the Colmeal catchment is small (10 ha), the knowledge created in this study could be scaled up to bigger catchments, which have often been highlighted to potentially cause serious off-site effects after a fire (Moody et al. 2013; Ebel 2020; Santi and Rengers 2020).

This study aimed to improve the understanding of the hydrological response of recently burnt catchments and, in particular, the role of SWR through the application of the process-based rainfall runoff OpenLISEM model.

After this modelling exercise, we conclude that calibration of OpenLISEM for post-fire hydrological response at catchment scale resulted in accurate simulations at the outlet, especially for quickflow, peak flows and timing of the peaks. Still, the poor validation results for quickflow and peak flows reveal OpenLISEM’s high sensitivity to specific event conditions and indicate that the number of calibrated events was not enough to perform a robust calibration.

Model calibration for two distinct initial SMCs showed lower model performance in predicting post-fire hydrological response for dry conditions, which presumably are under higher influence of SWR. Splitting the dataset for wet and dry conditions allowed determination of which other variables influenced model parameterisation besides initial SMC, and also allowed verification that future model calibration can benefit from such splitting.

The validation of the OpenLISEM spatially distributed runoff simulations revealed comparability problems and poor agreement with field surface runoff measurements within the monitored slopes. Additionally, field evidence shows that OpenLISEM failed to properly describe catchment connectivity between ploughed slopes and near streams, likely owing to an overestimation of the infiltration simulated in these slopes, and also because the model does not simulate baseflow.

This modelling exercise revealed that OpenLISEM has strong potential to simulate post-fire hydrological response, especially if an SWR component can be considered in the simulations, and if more detailed information is provided to the model, namely high resolution of SWR measurements and spatially distributed SMC.

The data that support the findings of this study are available from the corresponding author, Diana Vieira, on reasonable request.

The authors declare no conflicts of interest.

M. Basso is funded by FCT – Fundação para a Ciência e a Tecnologia, I.P., through a PhD fellowship (SFRH/BD/146356/2019). J. P. Nunes is funded by FCT through IF-FCT research grant IF/00586/2015, and J. J. Keizer is funded by FCT through IF-FCT research grant IF/01465/2015. This study was carried out in the framework of the COST Action ES1306: Connecting European Connectivity research, by the FCT national projects EROSFIRE II (PTDC/AGR-CFL/70968/2006) and FEMME (PCIF/MPG/0019/2017). Thanks are due for the financial support to CESAM (UIDP/50017/2020+UIDB/50017/2020), to FCT/MCTES through national funds (PIDDAC) and co-funding by the FEDER, within the PT2020 Partnership Agreement and Compete 2020.

Supplementary material is available online.