A storm event-based simulation system is developed to reproduce the effects of wildfire on a predefined forest area's rainfall infiltration and runoff. The simulation system reproduces individual surface cell slope, infiltration capacity, and hydrophobicity. A cellular automaton represents space and time. The Horton equation reproduces infiltration capacity evolution. Different equation parameters account for precipitation intensity variations. Infiltration volume and runoff velocity are computed at the cell level. A Digital Elevation Model reproduces the topography. A real rainfall event is reproduced for a section of La Primavera Forest in Jalisco, Mexico. For a subwatershed within the area of interest, wildfire effects are shown on infiltration, overland flow and outlet discharge throughout the rainfall event. It is shown that rainfall infiltration for wildfire affected terrain is lower and that overland flow and outlet discharge is higher. Infiltration and runoff dynamics at a surface cell depend on storm profile, cell slope, cell soil texture, and overall topography. Water dynamics on the area of interest are shown with a movie effect. It is concluded that a cellular automaton-based model can reproduce infiltration and runoff for different soil texture, topographic and rain conditions. This simulation system is geared for an optimisation system to pinpoint the locations of a series of forest land remediation controls that or maximise infiltration or minimise runoff.

Keywords: catchment peak discharge, cellular automaton, infiltrability, infiltration and runoff, wildfire, rainfall simulation

Numerous studies show that storm events high precipitation intensity cause high watersheds peak discharge, storm runoff or total runoff, and sediment yields.^1–5 The authors have carried out research to develop a simulation system to predict infiltration volume and infiltration excess overland flow before and after a wildfire has taken place. The model reproduces a real rainfall event, the forest topography, and individual forest surface cell slope and infiltrability. For a predefined subwatershed, the simulation system presents the wildfire effect on the catchment outlet discharge volume throughout the rainfall event. Rainfall infiltration and Hortonian overland flow are shown at the forest surface cell level throughout the rainfall event. Water runoff velocity is computed at the cell level with Manning’s equation adapted for surface runoff.⁶ The model is designed to reproduce rainfall evolution at a middle scale, not at a specific spot, nor at a large scale –which is normally treated by Geographical Information Systems. Cellular Automata (CA) applications to environmental modelling generate large-scale pattern from small-scale local processes and the resulting representation could provide simpler and more efficient approximations of spatially complex processes compared to continuous approximations.⁷ CA provide better results in terms of accuracy at a comparable computational cost with a continuous approach. CA have been widely used in environmental modelling to study spatially extended dynamics.^8–16 TreeMig, also utilises CA.¹⁷

This paper presents a study carried out to reproduce high severity wildfire effects on rainfall infiltration and runoff. The topography and soil infiltration profile of a forest section which is part of Bosque la Primavera in Jalisco, Mexico, were provided as input. A real rainfall event was reproduced. Simulation results obtained for a subwatershed located within the area of study are comparable with catchment peak flows of other studies’ field measurements with similar soils to this study. The primary aim of this work is model design, model development and tuning. This research is geared for a future heuristics-based optimisation system. This system will pin point the location of land remediation controls in a way that infiltration is maximised or Hortonian overland flow is minimised. As a separate objective, the optimisation system is to support reforestation projects: the arrangement of forest land remediation efforts is to be designed in a way that rainfall infiltration increases the moisture level of otherwise comparatively dry spots. This optimisation system –which is in process– will help prevent newly planted trees from reaching a permanent wilting point.

Assumptions

The present simulation system makes the following assumptions:

1. A normalised scale in meters is used for cell altitude. The altitude of the lowest terrain cell displayed by the 3D mesh graph is assumed to be zero.
2. A surface terrain devoid of vegetation cover is considered. Even though forest soil has no vegetation cover after a high severity wildfire, sediment loss is not considered for this initial version of the simulation.
3. The simulation system can be used on a forest topography different from Bosque La Primavera only if it is parametrised for the forest specific conditions.
4. The system does not take into account the geological soil structure other than soil texture and slope. Also, it is not considered that a particular surface cell slope direction –North, South, East or West– could influence rainfall infiltration or runoff.
5. Horton equation parameters¹⁸ used for unburned and high severity wildfire terrain, respectively, are calculated with data obtained at previous studies in forest soil with the same texture as Bosque La Primavera.
6. Even though the model has the flexibility of reproducing different infiltration capacity for distinct surface areas –up to the surface cell level- the present study considers high and uniform hydrophobic conditions after a high severity wildfire in the area of interest. It has been shown that high severity wildfire affected areas exhibit a higher presence of hydrophobic soils.^19,20
7. Infiltration capacity is assumed to be homogeneous within a forest surface cell, and hydrophobicity is assumed to be homogeneous across a cell.

Simulation settings

A real forest area which is part of Bosque La Primavera in Jalisco, Mexico, was used for this simulation study. The topography, rainfall precipitation data, and soil infiltration profile were provided as input for the simulation. This forest is selected because it is a protected natural area which is prone to fire and requires rapid intervention. Wildfire events occur on an annual basis. The vertices of the selected area are located at the following coordinates:  for the south-west vertex,   for the north-west vertex,  for the south-east vertex,  for the north-east vertex.

CA–based representation

The simulation utilises the CA approach which deals with the discrete elements that make up the phenomenon of interest. The CA paradigm may model the complexity of the behaviour of a phenomenon through the simple rules that govern its individual elements.⁷ For the simulation study herein presented the forest area surface cells are the constituent elements, which are obtained by dividing the area of interest into10,000 identical square-shaped units, 10 m large. A 3D mesh graph based on a DEM with 10 m horizontal resolution, 1 m vertical resolution, is created to represent the topography. The 3D mesh graph cells have a one to one relationship with the terrain surface cells. The cell under consideration – (i,j) – and the four cells of the Neumann neighbourhood – (i+1,j), (i,j+1), (i,j-1), (i-1,j) – are constantly shifting. The cellular automaton 10,000 cells are analysed every iteration. Figure 1 depicts a section of the 3D mesh graph. Thompson et al.²¹ indicate that a DEM with a resolution higher than 10 m would represent a more computationally intensive program execution without necessarily contributing to the simulation performance. This model uses a CA absorbing boundary: runoff that flows beyond the study area outermost surface cells is removed from the simulation. The CA can thus be defined as:

$CA=[R, L, Q, f, I]$

where:

$R=\left\{\left(i,j\right)|1\le i\le N,1\le j\le M\right\}$  a two dimensional  lattice.

L is either the Von Neumann or Moore neighbourhood.

Q is the finite set of cell state values.

$f:Q X Q^{\left|L\right|}\to Q$  is the state transition function.

$I:R\to Q$  is the initialisation function.

where:

$R=\left\{\left(i,j\right)|1\le i\le N,1\le j\le M\right\}$  is a two dimensional lattice.

L = Either the Von Neumann or Moore neighbourhood.

Q = The finite set of cell state values.

$f:Q X Q^{\left|L\right|}\to Q$  = The state transition function.

$I:R\to Q$  = The initialisation function.

The cells

Each automaton cell represents a portion of the terrain. Viewed from above, a surface cell is a square-shaped area. Surface cells differ in slope, soil texture, burned-unburned condition, and their specific location within an irregular topography. The finite set of cell states Q can be defined as follows:

$Q=\left\{W, T, S, A, I_c, F_c, R_o, L_f\right\}$

Where

$W$  is the cell cumulative water infiltration.

$T$  is the cell superficial water volume.

$S=S_{ns}{S}_{ew}$  is a vector defining the surface cell north-south slope and east-west slope.

$A=a_1, a_2, a_3, a_4$  is a vector defining the normalized altitude for all four cell vertices.

$I_c$  is the initial infiltration capacity for the cell soil texture.

$F_c$  is the final infiltration capacity for the cell soil texture.

$R_o$  is the runoff received by the cell within the current iteration.

$L_f$  is the elapsed time since rainfall or runoff was first received.

Where

$W$  The cell cumulative water infiltration.

$T$  = The cell superficial water volume.

$S=S_{ns}{S}_{ew}$  ; a vector defining the surface cell north-south slope and east-west slope.

$A=a_1, a_2, a_3, a_4$  ; a vector defining the normalized altitude for all four cell vertices.

$I_c$  = The initial infiltration capacity for the cell soil texture.

$F_c$  = The final infiltration capacity for the cell soil texture.

$R_o$  = The runoff received by the cell within the current iteration.

$L_f$  = The elapsed time since rainfall or runoff was first received.

The simulation system incorporates a discrete approach in the space and time scales. On the time scale, the simulation span is constituted by a succession of time steps or iterations. One iteration is equivalent to 1 minute of rainfall elapsed time. Update rules are defined within the CA approach to determine the changes which occur every iteration in the state of surface cells.⁷ Update rules consider the state of single cells, which are arranged on a grid, and specify the conditions that must be met for a cell state change to occur. Update rules may contemplate the physical laws that govern the relation between cells, or approximate relations between them. For example, Douvinet et al.²² used three hydrological update rules for flood assessment.

The new state for a particular cell is determined by its current state and neighbouring cells' state.⁷ This simulation system utilises the Neumann neighbourhood which defines four contiguous cells as the neighbours of the cell under consideration (Figure 1). Rainfall and runoff are the main factors that influence a cell state change. A cell state change in the form of cumulative infiltration is caused by water absorption. A cell state change in the form of superficial water volume is caused either by rainfall, infiltration or runoff. Every iteration all 10,000 surface cells are analysed to determine infiltration and runoff which take place, and the 3D mesh graph gets displayed. Although both topography and infiltration are three-dimensional phenomena, surface water flow is a quasi-horizontal process. Therefore, update rules operate in the two-dimensional domain.

Figure 1 Roving window on the 3D mesh surface.

Manning equation adapted for surface runoff

Previous research indicates that soil infiltrability is diminished for steep terrain,^23–25 and that forest severe wildfire increase runoff significantly at the steep slope.^26–29,7 This simulation system computes runoff velocity at the cell level with Manning equation adapted for surface runoff.⁶ The slope, rainfall intensity and post wildfire condition are thus taken into account. Manning equation adapted for surface runoff is described next.

$\upsilon =\left(1/n\right)d^{2/3}{S}^{1/2}$   (1)

where

$\upsilon$  : runoff average velocity.

$n$  : coefficient of soil roughness.

$d$  : flow depth.

$S$  : surface slope. Original and destination cell superficial water altitude is used to compute surface slope.

Where

$\upsilon$  = The runoff average velocity

$n$  = The coefficient of soil roughness

$d$  = The flow depth

$S$  = The surface slope. Original and destination cell superficial water altitude is used to compute surface slope.

Horton equation

Soil infiltration capacity is computed every iteration for each individual surface cell throughout the storm event using Horton's infiltration equation.¹⁸ This equation (equation 2) establishes that rainfall infiltration rate diminishes as the rainfall event lengthens until a minimum threshold is reached, which corresponds to the saturated hydraulic conductivity.³⁰

$f_p=f_c+(f_0-f_c)e^{(-kt)}$   (2)

Where

$f_p$  : forest surface cell infiltration capacity at time t.

$f_c$  :final infiltration capacity: that is, saturated hydraulic conductivity.

$f_0$  : initial infiltration capacity

$k$  : a constant for the rate of decrease in infiltration capacity.

$t$  : elapsed time since rainfall or runoff is first received by the forest surface cell.

Where

$f_p$  = The forest surface cell infiltration capacity at time t.

$f_c$  = The final infiltration capacity; that is, saturated hydraulic conductivity.

$f_0$  = The initial infiltration capacity.

$k$  = A constant for the rate of decrease in infiltration capacity.

$t$  = The elapsed time since rainfall or runoff is first received by the forest surface cell.

Graphic display

The Gnuplot graphic program facilitates 3D mesh display, which is used to represent the topography of the area under study. The gnuplot-iostream library was used within the C++ simulation system because it incorporates command-line driven 3D mesh graph plotting facilities. Forest surface cells have a one to one relationship with the 3D mesh cells. The simulation system displays the terrain surface once for every iteration. A darker blue shade on a surface cell reveals a higher water content. 3D mesh graph cell shading changes gradually producing a video effect.

Simulation main cycle

Algorithm 1 describes the main cycle of the simulation execution.

Infiltration process

Soil infiltration capacity is computed every iteration at the surface cell level throughout the storm event using Horton infiltration equation.¹⁸ This equation (equation 2) was selected because it provides a good fit for the infiltration time series obtained for this study.^30,31 Northern Colorado Front Range soil texture is similar to Bosque La Primavera soil texture and Benavides-Solorio and MacDonald³² runoff measurements were utilised to derive infiltration capacity time series for this study. A minute-by-minute infiltration capacity time series was calculated using equation 3 and Benavides-Solorio and MacDonald³² runoff time series for plot Bobcat No. 5, plot Bobcat No. 16 and plot Bobcat No. 12. Initial and final infiltration capacity values are extracted from the derived infiltration capacity time series. To find the diminishing constant that provides the best fit, a C++ computer program was developed to make an exhaustive search.

Benavides-Solorio and MacDonald³² Bobcat unburned plot No. 12 runoff time series was used to determine Horton equation parameters for unburned terrain and Bobcat plots No. 5 and No. 16 runoff time series were used to determine equation parameters for high severity wildfire terrain. For unburned terrain, the values obtained for initial and final infiltration capacity are1.3 and  and the diminishing constant is . For high severity wildfire terrain, the values obtained are as follows: (a) for plot Bobcat 5, initial and final infiltration capacity are and and the diminishing constant is ; (b) for plot Bobcat No. 16, initial and final infiltration capacity are 1.56 and  the diminishing constant is . Horton exponential decay curves for the same soil texture as Bosque La Primavera for high wildfire affected terrain (a) and unburned terrain (b). Bars show minute-by-minute precipitation intensity for the rainfall event which occurred on June 27, 2016.

Figure 2 shows the 92-min rainfall event precipitation profile reproduced by the simulation (bars) and Horton exponential decay curves for unburned terrain (Figure 2b) and wildfire affected terrain (Figure 2a). The data series obtained from the measurements at plot Bobcat No. 5 is used to derive the exponential decay curve in Figure 2.a.1. The data series obtained from the measurements at plot Bobcat No. 16 is used to derive the exponential decay curve in Figure 2.a.2. The simulation study considers a particularly severe wildfire throughout the study area and post fire hydrophobicity is assumed to be homogeneous within and across all surface cells. Huffman et al.¹⁹ and Benavides-Solorio and MacDonald³² performed analysis in sites with the same soil texture as Bosque La Primavera. They explain that burn severity is among the most important factors to control the strength of soil hydrophobicity. Bobcat plot No.12 is used for unburned terrain because it is more representative given the dry soil and high intensity rainfall conditions. For this study soil infiltrability is assumed to be homogeneous for all surface cells given the extensive ground cover.

$Infiltration=Rainfall-Runoff$   (3)

Figure 2 Infiltration capacity decay curves.

It should be indicated that the model herein presented includes versatility to reproduce the specific hydrophobicity of distinct zones within the area of study, up to the surface cell level. Also, the model may reproduce any infiltration pattern other than the Horton equation used in this study.

Data validity

Van de Giesen et al.³³ found that for shorter plots there is relatively higher runoff –and lower infiltration– than for longer plots. However, they indicate that the importance of this phenomenon depends of the plot length with respect to rain duration. They argue that for long or intense rainfall events both short and long plots reach equilibrium and the importance of spatial variability almost disappears. Stomph et al.³⁴ analysed the scale effects of Hortonian overland flow. They concluded that the scale effect is related to the differences in time needed to reach the equilibrium phase in the hydrograph. They performed analysis for  and hillslopes lengths and found that rainfall exceeding  causes the scale effect to gradually decline and eventually disappear. We consider that the  runoff analysis on  plots with rainfall intensity greater than  which was carried out by Benavides-Solorio and MacDonald (2001)– is valid for the  precipitation and  rainfall span reproduced in the  surface cells of this study because in both cases the duration and intensity is significant enough to reach equilibrium and several of the runoff and infiltration processes are similar for both scales.^33,34

Infiltrability adjustment for precipitation variations

Infiltrability is adjusted for precipitation intensity variations in this study. As proposed by Dunkerley³⁵ it is assumed that a more complete inundation of the soil surface takes place in periods of high precipitation intensity and that infiltrability diminishes. In periods following significant hiatuses in rainfall it is assumed that infiltrability recovers. A 5-min precipitation intensity moving average was calculated for the 92-min precipitation time series of the rainfall event reproduced in this study. A precipitation intensity burst is considered to take place when the moving average gets around twice the previous 5-min average and stays at a similar level for no less than five minutes. A precipitation intensity contraction is considered to take place when the moving average is 90 percent below the previous 5-min average, or lower, and stays at a similar level for no less than five minutes. The rainfall event reproduced in this study exhibits precipitation bursts at minute 25 and 78 and a precipitation contraction at minute 74.

Sample cells infiltration progression

The simulation system determines the infiltration and runoff progression for any surface cell of interest. Figure 3 shows the infiltration progression for two sample cells: (a) cell located at the hillslope (column 55, row 62, in reference to the north–west vertex of the area of study), (b) cell located at col 74, row 85, below the subwatershed discharge area). Infiltration volume determined for (a) surface cell at the hillslope (column 55, row 62); (b) surface cell below the subwatershed outlet discharge (column 74, row 85). Infiltration at the hillslope is constrained by infiltration capacity and surface water availability. As Figure 3a depicts infiltration diminishes when precipitation intensity declines. Except for rainfall start, infiltration below the outlet discharge area is only constrained by the infiltration capacity because it is located at a water accumulation spot.

Figure 3 Infiltration rate samples.

Rainfall runoff process

Superficial water amount which gets transferred from one surface cell (original cell) to a destination cell in one iteration span depends on runoff velocity, which is calculated with Manning Equation adapted for surface runoff.⁶ High runoff velocity takes place at the steep slope and on high severity wildfire affected terrain. The higher the runoff velocity, the subsequent lower superficial water volume available for infiltration. This is consistent with the findings of Campbell et al.,²⁶ Inbar et al.,²⁷ Benavides-Solorio and MacDonald³² Robichaud et al.,²⁸ Mondal et al.,²³ Arnous²⁴ Arnous²⁵ Robichaud et al.²⁹ Manning Equation adapted for surface runoff (equation 1) takes into account surface water depth, surface slope and the terrain roughness coefficient. This study uses a roughness coefficient value of 0.04 for high severity wildfire affected terrain and a value of 0.10 for unburned terrain. Algorithm 3 describes the runoff process:

Eligible water for runoff is calculated as a fraction of the total volume of water that would flow to a contiguous cell if it were not limited by the span of one iteration. This fraction, which is called the runoff percentage, gets computed for each destination cell. Runoff destination cells are determined by the ground level slope of the original cell. For cells with a slope in one cardinal point direction, a water column is transferred in a way that the destination cell superficial water altitude gets to be equal to the original cell water altitude, unless water eligible to flow runs out. For cells with a slope in two cardinal points direction, epsilon thick sheets of water are repeatedly transferred to the destination cell whose water altitude is lowest until the superficial water for the three cells is levelled, or superficial water eligible to flow runs out. This is similar to the "lowest cell first fill" approach used by Ghimire et al.³⁶ For cells with no slope, all four neighbour cells are destination cells. Epsilon thick sheets of water are repeatedly transferred to the neighbour whose water altitude is lowest until the superficial water for all five cells is levelled, or superficial water eligible to flow runs out. Guidolin et al.³⁷ use a faster water transfer algorithm, Parsons and Fonstad³⁸ developed a cellular automata model for surface water flow; however, their models are not focused on soil infiltration. Cirbus and Podhoranyi³⁹ developed CA to predict the spread of river water mass; however, their model is not focused on soil infiltration either.

Simulation animations

An internet link is provided to watch the 3D graph overland flow simulation, throughout June 27, 2016 storm event in Bosque La Primavera.⁴⁰ An infiltration simulation is also included. The simulation displays the 3D mesh graph repeatedly, once for every minute of rainfall elapsed time, producing a film effect. The rainfall intensity varies throughout the storm event (Figure 2) leading to minute-by-minute variations in cell surface water volume, cell infiltration and runoff. Cell-specific cumulative water amount is reflected by cell colouring: An increasingly dark shade depicts higher cell water amount. For the runoff animation, surface colouring depicts superficial cell-specific water volume. Water courses are shaped by cell colouring, turning darker as Hortonian overland flow finds its way to ravines and sub watersheds, accumulates within them, and flows downhill, replicating the way runoff is conducted in natural conditions. Figure 4 shows overland flow for high wildfire affected terrain at rainfall simulation minute 23, a moment of high precipitation intensity. Cell-specific water depth at storm minute 23 is shown for high severity wildfire terrain. June 27, 2016 storm event is reproduced. The infiltration animation shows infiltration within the shallow surface. It shows that rainfall permeates through the hydrophobic layer of high severity wildfire terrain at around minute 50 after rainfall start. The hydrophobic condition of high wildfire affected terrain and the high intensity of the rainfall event reproduced cause the high level of runoff and relatively low infiltration observed.^19,20

Figure 4 Cell-specific surface water.

Validation

Cumulative infiltration is shown for unburned and wildfire affected terrain at the sub watershed (a). Overland flow is shown for unburned and wildfire affected terrain at the sub watershed (b). June 27, 2016 rainfall event is reproduced. Figure 5 shows the subwatershed cumulative infiltration volume and overland flow volume for unburned terrain and high severity wildfire affected terrain throughout the storm event. It can be appreciated that while cumulative infiltration increase throughout the simulation span (Figure 5a) the watershed overland flow diminishes continuously after rainfall minute 46 (Figure 5b) when rainfall precipitation intensity goes down. Overland flow volume is significantly higher, and infiltration volume is significantly lower, for wildfire affected terrain due to its lower infiltration capacity. This is consistent with Doerr et at.⁴¹ Minute-by-minute discharge volume is shown for "El Coyote" subwatershed. June 27, 2016 rainfall event is reproduced for unburned terrain and high wildfire affected terrain. Subwatershed overland flow volume increases almost continuously for the rainfall event minute  when precipitation intensity is most significant (Figure 5b). The catchment overland flow volume peaks at simulation minute 36 and it amounts to  for unburned terrain, and 581,209 L for wildfire affected terrain, which is 37.1% higher. The subwatershed catchment outlet discharge (Figure 6) peaks at simulation minute 59 for unburned terrain, at 10, 270 L min^-1 For high severity wildfire terrain it peaks at minute 58, and it amounts to 14,467 L min^-1, which is 40.87% higher than unburned terrain. Storm-end cumulative catchment outlet discharge is 537, 292 L for unburned terrain, and 884, 796 L for wildfire affected terrain, which is 64.7% higher (Figure 7). This is consistent with Moody et al.^41–43

Figure 5 Cumulative infiltration and overland flow.

Figure 6 Catchment outlet discharge.

Figure 7 Cumulative outlet discharge.

Catchment cumulative outlet discharge for unburned terrain and high severity wildfire affected terrain. For wildfire affected terrain Figure 8 shows minute–by–minute evolution of rainfall, infiltration, and catchment outlet discharge volumes throughout the storm event. Instant catchment overland flow volume change (increase or decrease) is also shown. As long as precipitation volume is higher than infiltration volume, the catchment overland flow and/or the outlet discharge increases; indeed, overland flow volume increases dramatically at moments of very intense precipitation (around minute 5 and around minute 20 of this storm event). The catchment overland flow volume variation curve turns negative (the catchment overland flow volume decreases) when precipitation volume is lower than infiltration volume (mostly after minute 46). The infiltration volume curve approximates an exponential decay curve because infiltration is restricted by Horton infiltration capacity calculation. Infiltration volume is also restricted by superficial water availability. Catchment minute-by-minute total rainfall volume (bars), infiltration volume, and outlet discharge. Instant catchment overland flow variation –volume increase or decrease– is also shown. For unburned terrain the infiltration curve is higher than it is for wildfire terrain, which leads to lower subwatershed outlet discharge values. The lower infiltration volume exhibited by the simulation for high severity wildfire terrain is consistent with Huffman et al.¹⁹ and Benavides-Solorio and MacDonald.²⁰

Figure 8 Wildfire affected terrain water dynamics.

Infiltration comparison unburned vs high wildfire terrain

For the 4.82 ha subwatershed of this study the simulation shows that the cumulative infiltration amounts to 1,525,200 L for burned terrain at rain storm end. For unburned terrain, the cumulative infiltration is 1,872,710 L, which is  higher. This is consistent with diverse studies which show that in natural conditions burned terrain exhibits a diminished infiltration capacity.^44,1,5,42,43 Forest burning causes an expansion of runoff and catchment water discharge volume.^41–43 This study shows that at the rainfall event-end, the cumulative catchment discharge volume for unburned terrain is 537, 292 L, whereas it is 884,796 L for wildfire affected terrain, which is 64.7% higher. Figure 7 shows minute-by-minute watershed cumulative outlet discharge volume.

This simulation is considered satisfactory even for a 4.82 ha subwatershed. The infiltration volume calculated by the simulation is consistent with the subwatershed size and topography leading to catchment discharge volumes comparable with other studies field measurements (Table 1 and Figure 9). The simulation is relevant because it shows the extent to which infiltration, runoff and outlet discharge may be affected for a specific watershed after a wildfire disturbance, given the catchment topography and soil infiltration profile.

Figure 9 Peak discharge volume.

Catchment peak discharge

Several studies show that there is a relation between rainfall intensity and the watershed peak discharge, storm runoff or total runoff, and sediment yields.^1–5 Devastating wildfires have produced large amounts of runoff and enormous peak discharges which have led to water erosion with a proportion that ranges from a relatively small amount to several orders of magnitude higher; and even up to a 200–fold increase in the erosion rate.² Some of the known factors that influence peakflows are rainfall intensity, slope steepness and soil water repellency.^45,3 However, Kunze and Stednick⁵ found that maximum 30min rainfall intensity I₃₀ is a good predictor of the variability observed in peak discharges, storm runoff and sediment yields up to two years after the forest fire takes place; it has been used widely to explain the expansion in runoff levels.

Table 1 includes a comparison of catchment peak flows of different studies with similar soils to this study.^{26,46,47,27,2–4} Watershed areas differ significantly and unit area peak discharge can range from several litres⁴⁷ to a few cubic meters,^26,27,5 and even up to values as high as 24m³ S^–1 km^–2, as shown by Moody and Martin² Moody and Martin³ found a good relationship between rainfall intensity and unit-area peak discharges using data from different watersheds. They also found a threshold of peakflow expansion when 30min rainfall intensity reaches 10mmh^–1. This graph shows  and the corresponding peak discharge volume for this study (centre) and for t­he field observations carried out by Moody and Martin², and Kunze and Stednick⁴ which are indicated in Table 1. The storm event selected for this research has a 59 mm h^–1 I₃₀ and the model obtained a peakflow value of 5.00 m³ S^–2 km^–2 for the selected watershed, which is between 3.9 and 24m³ S^–1 km^–2, the values measured for and 90mm h^–1 I₃₀ storms, both of which took place after devastating Colorado wildfires.^2–4 Figure 9 shows a good fit for a distribution that includes  and the corresponding peak discharge volume for this simulation and for the field observations carried out by Moody and Martin et al.^2–4 The slope steepness data measured at the different watersheds was quite similar and it ranged from 20% to 35% considering the total watershed and the main channel. The slope for the subwatershed of this study is 28.9%. This research determines a 40.87% watershed peak flow expansion after a high severity wildfire, a value considered inside the range of observed field data. The model herein presented includes appropriate versatility to reproduce the topography, soil infiltration profile, overland flow, and specific hydrophobicity of distinct zones within the study area. Even though the model presents consistent results, as shown in the previous sections, model parametrisation needs to be carried out in order to use it for a forest topography different from Bosque La Primavera.

The simulation system presented provides a means to reproduce the spatio–temporal rainfall infiltration and runoff for a terrain surface with irregular topography. The CA approach provides a good means for rainfall runoff modelling as it flows on flat land and along natural waterways and canals. Rainfall infiltration evolution, represented by Horton equation (equation 2), keeps a close parallel with actual infiltration measurements which translates into realistic infiltration simulation dynamics over the area. Specific surface cell location, slope, soil texture and hydrophobicity determine infiltration rate and span (Figure 3). This leads to reliable infiltration and overland flow values at any location within the study area, both for unburned and high severity wildfire affected terrain. Minute-by-minute and total infiltration, overland flow, and catchment discharge values can be reproduced throughout a storm event and after it has ended (Figure 5,6,7,8).

The model developed is adequate to simulate spatio–temporal forest surface water dynamics. This initial version of the simulation program acknowledges no soil moisture vestige from previous rainfall events. Future versions of the simulation system will include inter-event surface cell moisture evolution, sediment loss and deposition, and subsurface runoff. The primary aim of this work is model design, model development and tuning. The model is suitable to provide a fitness function for a metaheuristic in charge of optimising land remediation efforts spatial distribution to maximise infiltration and minimise runoff for burned terrain. The development of this optimisation is in progress. The simulation model may also provide a fitness function for a metaheuristic to optimise the spatial distribution of land remediation efforts to support reforestation projects. Threshold soil moisture which helps maintain forest reforestation alive is to be attained. That is, the spatial distribution of forest land remediation efforts can be designed in a way that rainfall infiltration increases the moisture level of otherwise comparatively dry spots.

This work was carried out thanks to the support received from the Consejo Nacional de Ciencia y Tecnología (CONACYT), México. Bosque La Primavera topographic map was provided by Organismo Público Descentralizado Bosque La Primavera.

Javier Eugenio Vergara–Blanco. Carried out the research itself including state of the art research, model development from scratch, data collection, software development from the ground up, software debug, tune-up and execution, manuscript writing.

Jérôme Leboeuf–Pasquier. Provided advice on the following items: research method, model development, software coding and development, software tune up, data collection, and manuscript writing. Juan de Dios Benavides–Solorio. Provided advice on the following items: watershed research method and criteria, watershed modelling, wildfire effects on infiltration and runoff, manuscript writing.

The authors declare no conflict of interest.

1. Martin DA, JA Moody. Comparison of soil infiltration rates in burned and unburned mountainous watersheds. Hydrological Processes. 2001;15(15):2893–2903.
2. Moody JA, DA Martin. Initial Hydrologic and Geomorphic Response Following a Wildfire in the Colorado Front Range. Earth Surface Processes and Landforms. 2001a;(26):1049–1070.
3. Moody JA, DA Martin. Post-fire, rainfall intensity–peak discharge relations for three mountainous watersheds in the western USA. Hydrological Processes. 2001b;15(15):2981–2993.
4. Neary DG, GJ Gottfried, PF Ffolliott. Post-Wildfire Watershed Flood Responses. 2003. p. 1–8.
5. Kunze MD, JD Stednick. Streamflow and suspended sediment yield following the 2000 Bobcat fire, Colorado. Hydrological Processes. 2006;(20):1661–1681.
6. Hergarten S, HJ Neugebauer. Homogenization of Mannings formula for modeling surface runoff. Geodynamics- Physics of the Lithosphere. 1997;24(8):877–880.
7. Breckling B, G Pe’er, YG Matsinos. Cellular Automata in Ecological Modelling. Chapter 8. Fred Jopp, Hauke Reuter, & Broder Breckling, Berlin Heidelberg, Springer. 2011. p. 105–117.
8. Mathey A, J Nelson. Decentralized Forest Planning Models – a Cellular Automata Framework. Vol. 15, cap. 1. Klaus von Gadow, & Timo Pukkala, Dordrecht, Heidelberg, London, New York. 2008. p. 169–185.
9. Mathey A, E Krcmar, S Dragicevic, et al. An object-oriented cellular automata model for forest planning problems. Ecological Modelling, Elsevier. 2008;212(3-4):359–371.
10. Altartouri A, L Nurminen, A Jolma. Spatial neighborhood effect and scale issues in the calibration and validation of a dynamic model of Phragmites australis distribution - A cellular automata and machine learning approach. Environmental Modelling & Software. 2015;(71):15–29.
11. Ghisu T, B Arca, G Pellizzaro, et al. An optimal Cellular Automata algorithm for simulating wildfire spread. Environmental Modelling & Software. 2015;(71):1–14.
12. Ortigoza GM, A Lorandi, I Neri. ACFUEGOS: An unstructured triangular cellular automata for modelling forest fire propagation. In: 6th International Conference on High Performance Computer Applications, ISUM 2015. Mexico City, Mexico. 2016. p. 132–143.
13. Sun T, L Zhang, W Chen, et al. Mountains Forest Fire Spread Simulator Based on Geo-Cellular Automaton Combined With Wang Zhengfei Velocity Model. IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing. 2013;6(4):1971–1987.
14. Vahidnia MH, AA Alesheikh, S Behzadi, et al. Modeling the spread of spatio-temporal phenomena through the incorporation of ANFIS and genetically controlled cellular automata: A case study on forest fire. Scopus. 2013;6(1):51–75.
15. Pak SI, T Hayakawa. Forest fire modeling using cellular automata and percolation threshold analysis. San Francisco, CA; United States. 2011. p. 293–298.
16. Jellouli O, A Bernoussi, M Amharref, et al. Modeling of wind flow and its impact on forest fire spread: Cellular automata approach. In: International Conference on Cellular Automata. Fez, Morocco. 2006. p. 269–279.
17. Lischke H, NE Zimmermann, J Bolliger, et al. TreeMig: A forest-landscape model for simulating spatio-temporal patterns from stand to landscape scale. Ecological Modelling, ELSEVIER. 2006;199(4):409–420.
18. Horton RE. An approach toward a physical interpretation of infiltrationi capacity. Soil Science Society Proceedings. 1940;(5):399–417.
19. Huffman EL, LH MacDonald, JD Stednick. Strength and persistence of fire induced soil hydrophobicity under ponderosa and lodgepole pine, Colorado Front Range. Hydrological Processes. 2001;(15):2877–2892.
20. Benavides-Solorio, J de D, LH MacDonald. Post-fire runoff and erosion from simulated rainfall on small plots, Colorado Front Range. Hydrological Processes. 2001;15(15):2931–2952.
21. Thompson JA, JC Bell, CA Butler. Digital elevation model resolution: effects on terrain attribute calculation and quantitative soil-landscape modeling. Geoderma. 2001;100(1-2):67–89.
22. Douvinet J, MJ Van De Wiel, D Delahaye, et al. A flash flood hazard assessment in dry valleys (northern France) by cellular automata modelling. Natural Hazards. 2015;75(3):2905–2929.
23. Mondal S, AC Pandey, RD Garg. Groundwater Prospects Evaluation Based on Hydrogeomorphological Mapping using High Resolution Satellite Images: A Case study in Uttarakhand. Indian Society of Remote Sensing. 2008;(36):69–76.
24. Arnous MO. Groundwater potentiality mapping of hard-rock terrain in arid regions using geospatial modelling: example from Wadi Feiran basin, South Sinai, Egypt. Hydrogeology Journal. 2013;24(6):1375–1392.
25. Arnous MO, DR Green. GIS and remote sensing as tools for conducting geo-hazards risk assessment along Gulf of Aqaba coastal zone, Egypt. Journal of Coastal Conservation. 2011;15(4):457–475.
26. Campbell RE, MB Baker Jr, PF Ffolliott, et al. Wilfire effects on a ponderosa pine ecosystem: An Arizona case study. U.S. Department of Agriculture, Forest Service, Rocky Mountain Forest and Range Experiment Station., Research Paper RM-191, Fort Collins, CO, USA. 1977. p. 12.
27. Inbar M, M Tamir, L Wittenberg. Runoff and erosion processes after a forest fire in Mount Carmel, a Mediterranean area. Geomorphology. 1988;24(1):17–33.
28. Robichaud PR, FB Pierson, RE Brown, et al. Measuring effectiveness of three postfire hillslope erosion barrier treatments, western Montana, USA. Hydrological Processes. 2007;22(2):159–170.
29. Robichaud PR, SA Lewis, JW Wagenbrenner, et al. Post-fire mulching for runoff and erosion mitigation Part I: Effectiveness at reducing hillslope erosion rates. Catena. 2013;(105):75–92.
30. Hillel D. Environmental Soil Physics (1st ed.). Academic Press. San Diego CA. USA. 1998. p. 771.
31. Brooks KN, PF Ffolliott, JA Magner. Hydrology and the Management of Watersheds. Wiley-Blackwel. Chichester, West Sussex, UK. Oxford, UK. Iowa, USA. 2012. p. 552.
32. Benavides-Solorio, J de D, LH MacDonald. Measurement and prediction of post-fire erosion at the hiilslope scale. Colorado Front Range. International Journal of Wildland Fire. 2005;(14):457–474.
33. Van de Giesen NC, TJ Stomph, N de Ridder. Scale effects of Hortonian overland flow and rainfall - runoff dynamics in a West African catena landscape. Hiydrological Processes. 2000;(14):165–175.
34. Stomph TJ, N de Ridder, TS Steenhuis, et al. Scale Effects of Hortonian Overland Flow and Rainfall-Runoff Dynamics: Laboratory Validation of a Process-Based Model. Earth Surface Processes and Landforms. 2002;(27):847–855.
35. Dunkerley D. An approach to analysing plot scale infiltration and runoff responses to rainfall of fluctuating intensity. Hydrological Processes. 2017;31(1):191–206.
36. Ghimire B, AS Chen, M Guidolin, et al. Formulation of a fast 2D urban pluvial flood model using a cellular automata approach. Journal of Hydroinformatics. 2013;15(3):676–686.
37. Guidolin M, AS Chen, B Ghimire, et al. A weighted cellular automata 2D inundation model for rapid flood analysis. Environmental Modelling & Software. 2016;(84):378–394.
38. Parsons JA, MA Fonstad. A cellular automata model of surface water flow. Hydrological Processes. 2007;(21):2189–2195.
39. Cirbus J, M Podhoranyi. Cellular Automata for the Flow Simulations on the Earth Surface, Optimization Computation Process. Applied Mathematics & Information Sciences. 2013;7(6):2149–2158.
40. Vergara-Blanco JE. Internet link for simulation animations. 2020.
41. Doerr SH, RA Shakesby, LH MacDonald. Soil Water Repellency: A Key Factor in Post-fire Erosion? Vol. Fire Effects on Soils and Restoration Strategies., cap. 7, Artemi Cedra, & Peter M. Robichaud, 2009. p. 197–223.
42. Moody JA, DA Martin, SL Haire, et al. Linking runoff response to burn severity after a wildfire. Hydrological Processes. 2008;22(13):2063– 2074.
43. Coombs JS, JM Melack. Initial impacts of a wildfire on hydrology and suspended sediment and nutrient export in California chaparral watersheds. Hydrological Processes. 2013;27(26):3842–3851.
44. Pierson FB, PR Robichaud, KE Spaeth. Spatial and temporal effects of wildfire on the hydrology of a steep rangeland watershed. Hydrological Processes. 2001;15(15):2905–2916.
45. DeBano LF, DG Neary, PF Ffolliott. Fire's effects on ecosystems. John Wiley & Sons, Inc. New York, NY, USA. 1998. p. 352.
46. Scott DF, RE Schulze. The Hydrological Effects of a Wildfire in a Eucalypt Afforested Catchment. South African Forestry Journal. 2010;160(1):67–74.
47. Troendle CA, GS Bevenger. Effect of fire on streamflow and sediment transport in Shoshone National Forest, Wyoming. In: Second Biennial Conference on the Greater Yellowstone Ecosystem, Yellowstone National Park, Wyoming, USA COLOCAR AQUÍ DATOS EDITORIALES Y PÁGINAS. 1996.