Global sensitivity analysis and scale effects of a fire propagation model used over Mediterranean shrublands

Ecological Modelling 136 (2001) 175– 189 www.elsevier.com/locate/ecolmodel

Global sensitivity analysis and scale effects of a fire propagation model used over Mediterranean shrublands R. Salvador a,*, J. Pin˜ol a, S. Tarantola b,1, E. Pla a a

Center for Ecological Research and Forestry Applications (CREAF), Uni6ersitat Auto`noma de Barcelona, 08193 Bellaterra, Barcelona, Spain b Joint Research Centre of the European Commission, ISIS, TP 361, 21020 Ispra (Va), Italy Accepted 12 October 2000

Abstract A Global Sensitivity Analysis (GSA) and an analysis of scale effects have been carried out over the equations given by Rothermel (1972), with some additional modifications. Data mainly derived from Mediterranean shrublands and a spatially close meteorological station have been used to derive probability distribution functions for the variables involved. In spite of the abundant non-linearities contained in the equations studied, scale effects found have been relatively unimportant, thus supporting the current usage of spatial averages (i.e. fuel models) in spatially distributed models based on the Rothermel equations. On the other hand, the results of the GSA clearly showed the negligible effect of the variability of three of the model variables (the low heat content, the particle density and the mineral content) on the output values. However, all other input variables had some noticeable effect over the variability of the output, and they cannot be ignored if an optimal use of the model is desired. Finally, caution is recommended if results have to be extrapolated to other types of vegetation or climate. © 2001 Elsevier Science B.V. All rights reserved. Keywords: Rothermel equations; Global sensitivity analysis; Scale effects; Fire propagation; Mediterranean shrublands

1. Introduction Models are being increasingly used to make predictions in many fields of the environmental sciences. Models, though, have a series of limitations from which the lack of accuracy and/or a * Corresponding author. Tel./fax: +34-93-5811312. E-mail addresses: [email protected] (R. Salvador), [email protected] (S. Tarantola). 1 Tel.: + 39-0332-785639; fax: + 39-0332-785733.

proper validation have the most dramatic consequences. However, other restrictions do generally occur when applying models. Frequently, information on a high number of input variables is required to run a model, this need being specially restrictive for models with a spatially explicit expression (since these require data for variables in each site). Fire propagation models are a good example of data demanding (data restricted) models. Although very different approaches have been used

0304-3800/01/$ - see front matter © 2001 Elsevier Science B.V. All rights reserved. PII: S 0 3 0 4 - 3 8 0 0 ( 0 0 ) 0 0 4 1 9 - 1

176

R. Sal6ador et al. / Ecological Modelling 136 (2001) 175–189

to model the spread of forest fires, these may be grouped in three main categories (Perry, 1998): the physical, the empirical and the semi-physical models. Physical models heavily rely on theoretical aspects of fire propagation such as the heat transfer process (see Frandsen, 1971; Albini, 1985; Weber, 1989 as instances). Nevertheless, while the theoretical aspects may bring some insight in the mechanisms of fire spread, physical models are not currently operational as they require information of many parameters, which are almost impossible to be measured in the field (Beer, 1991). In addition, the heat transfer processes modelled are neither spatially nor temporally constant. On the other hand, some authors have developed empirical models to overcome the problems derived from the complexity of the theoretical approaches. These models have been largely developed in Australia (McArthur, 1966; Green et al., 1983; Cheney et al., 1993) and Canada (Stocks et al., 1989; McAlpine, 1995) where extensive databases of wild and prescribed fires are available. Such empirical models are, in essence, statistically fitted models, and as a main drawback they should be used with great caution under other conditions than those of the original fires. Finally, the semi-physical approach combine physical and empirical techniques, looking for operational models which can be used under a wide range of conditions. Specifically, while the main structure of model equations is based on physical relationships these equations are adjusted and modified considering the results of experimental studies in the laboratory. Hence, due to their intermediate features, semi-physical models are currently applied with predictive purposes in real fire events. The model developed by Rothermel (1972) is, by far, the most used semi-physical model. In fact, although there are several spatially explicit models being currently applied to predict the behaviour of wild and prescribed fires, such as the FARSITE (Finney, 1998), the FIRESTATION (Lopes et al., 1998) and the CARDIN simulators (Martı´nezMilla´n and Saura, 1998), the majority have in their core the set of equations given by Rothermel.

However, such spatially explicit models have been applied rather uncritically, and there is still a lack of evidences of their reliability, being probably limited by several aspects. Among these aspects, two are directly related to the model development and structure: (a) the arising of spatial scale effects; and (b) a significant number of input variables in the model. (a) Rothermel equations were derived empirically in a laboratory at a spatial scale of few square meters but are now currently used to predict the fire behaviour over wide areas. Since exhaustive spatial information of variables involved is not available, a simplified approach based on maps has been usually used. These maps include patches (e.g. polygons in a geographical information system) to which an average value is assigned to the input variables. Indeed, estimates of such averages are readily available in the literature for the main vegetation types (these sets of estimates are known as fuel models). However, Rothermel equations contain many non-linear functions of the input variables and such non-linearities may lead, by the use of mean values through wide areas with some internal variability, to biased outputs (i.e. to noticeable spatial scale effects). (b) Data for many input variables is needed to apply the Rothermel equations, and some of these variables may be very costly and difficult to obtain. Hence, a reduction in the number of variables is highly advisable in order to achieve an operational use of the model. This requires a previous knowledge of the relative influence of all input variables, which may be attained through an analysis of the sensitivity of the model output. Although some attempts considering the combined variability of one or few variables have been made on the Rothermel equations (Salazar, 1985; Burgan, 1987), they have shown not to be adequate to the task since equations involved are complex and variables interact. Consequently, other methods allowing the joint variability of all input variables during the assessment of their relative weights, such as the Global Sensitivity Analysis (Saltelli et al., 2000), will be more suitable.

R. Sal6ador et al. / Ecological Modelling 136 (2001) 175–189

In this study (i) we analyse the possible emergence of scale effects due to the supposition of a homogeneous spatial value for the input variables of the Rothermel equations; and (ii) we carry out a Global Sensitivity Analysis of these equations to look for the most influential variables. Data and information used mainly comes from Mediterranean shrublands.

2. Methodology

2.1. A summary of the fire propagation model The Rothermel model (Rothermel, 1972) has two main outputs: the rate of spread of a point in

177

the fire front (R, given here in cm s − 1) and the reaction intensity (IR), which is the energy released per unit area for a specific period of time (given here in kW m − 2). Since 1972, several modifications have been proposed for these equations (Wu et al., 1996). In this study we have included those of Albini (1976), affecting the net fuel loading and the optimum reaction velocity, and those of Catchpole and Catchpole (1991), based on Wilson (1990) which eliminate the use of the moisture content of extinction and give a new formula for the heat of preignition. The list of the 10 input variables considered in the model is given in the first column of Table 1. The main equations of the model analysed are given in Eqs. (1) and (2) (see Rothermel, 1972 for a detailed description of

Table 1 Some information concerning the 10 input variables finally included in the analysesa. The variable Fuel loading (w0 ) is not included in the table as it was given as a function of the Fuel depth (l) (see Fig. 1) Variable

Symbol and units

p.d.f.

Fuel depth Fuel particle area-to-volume ratio Fuel particle low heat content

l (cm) | (cm−1)

Log N(3.06, 0.517) Log N(3.31, 0.928)

h (Kcal kg−1)

N(4842, 308)

4842

6.4

Literature

zp (D.W. g cm−3)

Log N(8.48, 0.063) N(0.567, 0.125)

4842 0.567

6.4 0.125

Literature

Log N(−0.592, 0.219)

0.567

0.125

N(1.18, 0.377)

1.18

31.9

508 field measurements

N(0.8, 0.112) md (H2O g D.W. g−1) Semi-empirical

0.8 0.19

15 24.9

ST (min. gD.W. g−1)

N(0.049, 0.011)

0.049 22.4

Simulated scenario 16383 measurements on relative humidity applied to an empirically fitted model Measured in 36 individuals

U (km h−1)

Empirical

3.23

tan ƒ P

U% = 6.9U N(0.38, 0.186) Log N(−2.19, 0.64)

6.19 60.3 0.38 48.4 0.138 71.7

16383 hourly data (correction factor: 0.4) Simulated scenario 1517 points derived from a MDE The same as for the fuel depth

N(0.5, 0.15)

0.5

Simulated scenario

Ovendry particle density

Moisture content of the live fuel Moisture content of the dead fuel Fuel particle total mineral content Wind speed at midflame height Slope Dead fuel loading vs. total loading

ml (H2O g D.W. g−1)

Mean

C.V. (%)

Sources of information

24.3 28.6

55.3 30

Data on 23 field plots (3×3 m) Literature and theory

60.3

30

a D.W., dry weight; p.d.f., probability density function; N, normal p.d.f.; Log N, lognormal p.d.f.; C.V., coefficient of variation; min., mineral weight

R. Sal6ador et al. / Ecological Modelling 136 (2001) 175–189

178

the original equations and the way they were obtained). IR = G%wn hpM ps R=

IRx(1+ƒw + ƒs) zbmQig

(1) (2)

All terms included are, in turn, functions of the 10 input variables listed in Table 1: Y%(|, l, w0, zp): optimum reaction velocity (s−1). wn(w0, ST): net fuel loading (kg m−2). h: fuel particle low heat content (kcal kg−1). pM(ml, md, P): moisture damping coefficient. ps(ST): mineral damping coefficient. x(|, l, w0, zp): propagating flux ratio. ƒw(|, l, w0, zp, U): wind coefficient. ƒs(l, w0, zp, tan ƒ): slope coefficient. zb(l, w0, zp): ovendry bulk density (kg m−3). m(w0): effective heating number. Qig(md): heat of pre-ignition (kJ kg−1).

Since the fraction of fuels with smaller diameters have the highest influence in the fire behaviour (Rothermel, 1972) only fuels with a diameter smaller than 6 mm (which are known as 1 h fuels) were considered.

2.2. The geographical scope of the study Data required to carry out this study had different origins, mainly field sampling, meteorological measurements, and information from literature. The first type of data was raw data gathered in previous studies and unpublished works, carried out in several shrublands located in the north east coastal region of the Iberian Peninsula. These shrublands can be included in the fuel models 5 (brush) and 6 (chaparral) of the fuel model classification given by Rothermel (Anderson, 1982). Meteorological data were retrieved from a meteorological station located in the same area. In consequence, the conclusions derived from this study will be mainly concerned with fires on Mediterranean shrublands.

2.3. The estimation of the scale effects Rothermel equations suppose stationarity of the fire front (i.e. they are not intended to explain the fire behaviour while input variables vary). Therefore, besides some restrictions usually applied to ensure a continuous fire front when using spatially distributed models (e.g. Huygens’ ellipses), fire evolution in 2-D fire spread models is basically treated as an additive process of stationary states. Hence, the order of occurrence of such states will not be crucial for the analysis of scale effects. As pointed out in the introduction, the use of mean values through wide areas with some internal variability (i.e. the fuel model approach) may lead to scale effects. The occurrence of such effects can be checked by means of simulations. Specifically, the output derived from the use of the mean of the input variables (which will be called uncorrected mean) should be compared to the average of the output values obtained by the simulations (the corrected mean). If there are no scale effects the following equality will be held: E[S(X1, X2, …, Xn )]= S(E[X1], E[X2], …, E[Xn ]). (3) That is to say, the expected value (E) of the output of the model S (i.e. the average of all simulations on R or IR) will equal the output obtained when applying the mean value of each input variable (E[Xi ]) to the equation S. The probability density functions (p.d.f.s) of all input variables (X1, …, Xn ) are needed to carry out the simulations. The way these functions were obtained is explained in Section 2.5 below, 106 simulations were carried out for each analysis of scale effects.

2.4. The global sensiti6ity analysis The objective of sensitivity analysis is to assess the importance of each input variable on the values of the output variable of the model. Among the many ways to perform such analysis (see Helton, 1993; Hamby, 1994 for a revision) the global sensitivity analysis (GSA) allows quantifying the significance of interactions between the

R. Sal6ador et al. / Ecological Modelling 136 (2001) 175–189

input variables which, due to the high complexity of many models, will frequently be of paramount importance. The most effective way to implement GSA is based on the decomposition of the variance of the output in a unique additive series of variances of independent functions of increasing dimensionality Var[S(X1, …, Xn )] =% Var[Ai (Xi )]+ % Var[Bij (Xi, Xj )] i

+ %

iBjBk

iBj

Var[Cijk (Xi, Xj, Xk )] +···

+ Var[H(Xi, Xj, …, Xn )].

(4)

The terms Var[Ai (Xi )] are the variances due to the main effects for each one of the input variables, the terms Var[Bij (Xi, Xj )] are the variances due to the pairwise interactions between all pairs of input variables, the Var[Cijk (Xi, Xj, Xk )] give the effect of the three component interactions, and so on until the term considering the interaction effects of all n variables together is reached (Var[H(Xi, Xj, …, Xn )]). This decomposition permits assessing the effect on the model output uncertainty of each input variable Xi through the variances of the main effects, and of any group of input variables through higher order variances (interaction effects). For instance, the second order variance Var[Bij (Xi, Xj )] represents the contribution to the model output uncertainty due to (Xi, Xj ) that cannot be explained by the sum of the individual effects Var[Ai (Xi )] and Var[Aj (Xj )]. This variance decomposition has its origin on an expansion of the Haje´k projection (Haje´k, 1968) proposed by Efron and Stein (1981) and which has been described more extensively in sections A1 and A2 of the Appendix A. Various estimators of the terms in Eq. (4) are available: some have been derived from the analysis of frequencies (Cukier et al., 1975; Saltelli et al., 1999), others are based on Monte Carlo integration (Sobol’, 1990). Due to the high number of interactions arising even for a moderate number of independent variables, only the value of few of the variances is

179

usually derived. For each input variable, we have calculated the variance of the main effect and, following Homma and Saltelli (1996), the sum of variances in Eq. (4) where the variable was present (i.e. its main effect plus its interaction terms, which is named total effect). Specifically, a small value for the total effect will be enough to state that the variable analysed is not important, while main effects will be useful to quantify the direct effect of variables over the output (small values of the main effects are not sufficient guarantee to reject a variable as it may have a strong influence through the interactions). Finally the variances calculated for the main and total effects have been divided by the total model output variance, obtaining the so-called sensitivity indices, which have a much easier understanding, and allow the variables to be ranked quantitatively. However, functions involved in the computations of the main and total effects may be very complex to be analytically derived, and numerical methods are frequently used. In section 3 of the Appendix A there is a description of the procedure to estimate these quantities. The same simulations carried out for the analyses of the scale effects (see Section 2.3) were used in the GSA (the p.d.f.s required are described in Section 2.5). The main drawback of the variance-based GSA is its requirement for stochastic independence of the input variables. To minimise this restriction we have either supposed total independence or total dependence (by means of a mathematical function fitted empirically), depending on the degree of tightness between pairs of variables.

2.5. Estimation of the probability density functions As mentioned in Section 2.2, estimates of the p.d.f.s of input variables were derived from data coming from weather stations and extensive field measurements previously obtained in several shrublands. Information from literature was also used when necessary. All relevant features of the input variables and their p.d.f.s are summarised in Table 1.

180

R. Sal6ador et al. / Ecological Modelling 136 (2001) 175–189

Fig. 1. Bivariate scatterplot of the fuel loading and the height of the fuel bed for 23 plots of a Mediterranean shrubland (row data from Pla and Roda`, 1999). Since there is a clear dependence between both variables, a sigmoidal curve has been fitted to the data.

Fig. 2. Histograms for 16 383 values of: (a) the relative humidity of the air; and (b) the wind speed (in km h − 1) measured in a meteorological station located in the study area.

The fuel depth and the fuel loading were the only pair of variables considered to be too close to be treated as completely independent. Data for both variables had been obtained from the same field plots (Pla and Roda`, 1999), and a sigmoidal curve was fitted to them (Fig. 1). By means of this curve, the fuel loading could be given as a func-

tion of the fuel depth values obtained in the simulations. Since dead and live fuels were measured separately the proportion of dead vs. the total fuel could be derived for each field plot. The live fuel moisture content had been measured in the field (Viegas et al., in press) and the dead fuel moisture content was indirectly obtained from measurements of the relative humidity of the air (Fig. 2a) through an empirical model (Pook and Gill, 1993). Hourly averages of wind speed were obtained from the same weather station, and due to the difficulty to fit a theoretical p.d.f. to its histogram (Fig. 2b) values for the simulations were retrieved directly from the sample. As suggested by Andrews (1986) for the fuel types 5 and 6, a correction factor of 0.4 was applied to obtain the wind speed at midflame height. Data of slopes were extracted by sampling from a file of slopes derived from a 45 m resolution Digital Elevation Model (DEM) of a wide area of the north-east of the Iberian Peninsula. Finally, the fuel particle total mineral content had been previously obtained from 36 individuals of different Mediterranean shrub species (Paradell, unpublished data). Data for the fuel particle low heat content, the ovendry particle density and the fuel particle surface-area-to-volume ratio were not measured. Since the low heat content and the particle density are fairly constant within species, we used the lists of values given by Elvira and Hernando (1989); Panov (1994), Papio´ (1990) on Mediterranean species to derive a mean and a standard deviation for them. Because p.d.f.s for both variables could not be guessed from the literature we applied both the normal and the lognormal p.d.f.s. A similar procedure was carried out for the surface-area-to-volume ratio. The data from the revision on Mediterranean species made by Panov (1994) was used and a lognormal p.d.f. was chosen after some theoretical considerations. Since the number of simulations carried out was high (106 simulations for each analysis) limits were set to all simulated values of the input variables to avoid feeding the model with extremely unrealistic values.

R. Sal6ador et al. / Ecological Modelling 136 (2001) 175–189

2.6. Additional scenarios Data for all input variables of the model were gathered all over the year and, in consequence, results obtained in the analyses have to be considered in a yearly basis. However, wild fires (especially those affecting large extensions) frequently occur under specific vegetation and/or meteorological conditions. In order to assess the behaviour of the model when these conditions are met we have carried out some additional analyses. Specifically, we have simulated three new scenarios modifying the original p.d.f.s of some variables: 1. Strong wind: Uncorrected wind speeds have been multiplied by 6.9 to have 10% of the observations with values higher than 100 km h − 1. As a consequence, the mean uncorrected wind speed has raised from 8.1 to 55.2 km h − 1. The rate of spread has been the only output variable analysed because the wind

181

speed is not included in the formula of the reaction intensity (see Eq. (1)). 2. High water deficit: High water deficit conditions have been simulated by: (a) reducing the mean value of the fuel particle live moisture content from 1.18 to 0.8 (and the coefficient of variation from 31.9 to 15%); and by (b) increasing the proportion of the dead fuel loading from 0.138 to 0.5 (with a new coefficient of variation of 30%). 3. Strong wind and high water deficit: Modifications of the two previous scenarios were combined in a new scenario. Again, as IR results were not affected by the wind, only results for R were analysed. All three scenarios were considered for both the analyses of scale effects and the GSA.

3. Results Histograms of the outputs of R and IR carried out under the original conditions are clearly asymmetrical (Fig. 3), with highest frequencies at low values and a long tail spanning to the right side (i.e. most of the fires simulated have low rates of spread and intensities, and only few of them become really active and intense). These results are in concordance to what could be expected from ignitions occurring all over the year, since just a few of them will become active, dangerous fires.

3.1. Results of the analyses of scale effects In all the scenarios there are divergences between the corrected and the uncorrected means (Table 2). Specifically, while the R is systematically underestimated by the uncorrected means, the IR shows an opposite behaviour. However, in all cases the absolute degree of the bias observed is similar and rather low (with all differences lower than 25%). Fig. 3. Histograms of the outputs of 5000 simulations considering the original conditions for: (a) the rate of spread (R, given in cm s − 1); and (b) the reaction intensity (IR, given in kW m − 2).

3.2. Results of the GSA In all analyses carried out (Fig. 4 and Fig. 5) there are three variables (h, zp and ST) in which

R. Sal6ador et al. / Ecological Modelling 136 (2001) 175–189

182

Table 2 Values of the corrected and uncorrected means obtained under the original conditions and the three additional scenarios Scenario

Output variable

Original

R IR

Strong wind

R

Uncorrected mean

Corrected mean

Ratio uncorrected–corrected

0.371 cm s−1 38.6 kW m−2

0.464 cm s−1 32.9 kW m−2

0.800 1.173

1.61 cm s−1

2.07 cm s−1

0.778

−1

−1

High water deficit

R IR

1.08 cm s 109.1 kW m−2

1.27 cm s 91.2 kW m−2

0.850 1.196

Wind+water deficit

R

4.71 cm s−1

5.63 cm s−1

0.837

their variability has no noticeable weight over the variability of the model output. In fact, this result is rather stable since it is found in all scenarios, it occurs for R and IR, it happens in the main effects and in the total effects, and it is not affected when the p.d.f.s are modified. However, the rest of the input variables usually contribute to some extent to the variability of the outputs. Another consistent result concerns the similarities between the main and the total effects. Thus, for each scenario and output variable, the rank and relative degree of importance among variables is always kept when comparing main and total effects (obviously, total effects are always equal or bigger than main effects). The GSA carried out over R, under the original conditions, points to l, md and ml (given in this order) as the most influential variables (see Fig. 4a and Fig. 5a), while the results for IR are quite similar, but in this case ml is set as the dominant input variable (Fig. 4d and Fig. 5d). In both analyses the main effects account for a big portion of the variance of the output (almost 70% in R and 85% in IR) and little variability remains to be assigned to the interaction effects. The increase in wind speeds simulated under the strong wind scenario leads to an enlargement of the importance of the wind speed (U) (see Fig. 4a and Fig. 5a). This variable becomes as important as l in the main effects, and is the most influential variable when the total effects are analysed. However, the effects of the increase in the wind speeds are not easily explained. Indeed, they are not directly proportional to the multiplying factor (U%=6.9*U) and the interactions have

now a stronger presence than they had under the original conditions (the main effects explain here only 50% of the output variance). In addition, the effect of | is much noticeable than before, becoming more important than both measures of moisture content. High water deficit levels lead to changes in the relative importance of md and ml (Fig. 4b, Fig. 4d, Fig. 5b and Fig. 5d). On one hand the effect of ml has almost disappeared, probably as a consequence of two features of this scenario. Specifically, the mean value of P (the proportion of dead fuels) has been raised to 0.5 and, although the mean value of ml has been lowered to 0.8 (leading to faster and more intense fires), its variance has been lowered too. On the other hand, the importance of md has increased (probably due to the same reasons). These results suggest that the behaviour of fires occurring under high water deficit levels will be mainly controlled by the relative humidity of the air and by the values of l. The amounts of variance explained by the main effects are similar to these found under the original conditions (72% for R and 87% for IR). Finally, the results for the scenario with high wind and high water deficit are a mixture of those obtained in the previous two paragraphs (see Fig. 4c and Fig. 5c). Thus, while the effect of U and | has been probably risen as a consequence of the increase in the wind speeds, the weight of ml drops down to very low values as it happened in the scenario with a high water deficit. Furthermore, the amount of variance explained by the main effects has also a value between those of the high wind and the high water deficit scenarios (62% for R).

R. Sal6ador et al. / Ecological Modelling 136 (2001) 175–189

183

4. Discussion In spite of the abundant non-linearities contained in the model analysed, the scale effects found have been relatively low. Moreover, the amount of variability used in the simulations should be considered as an upper bound for the variability to be

Fig. 5. Results of the GSA concerning the total effects. Sensitivity indices for the original conditions and the additional scenarios are compared in the different graphics: wind on R (a), water deficit on R (b), wind and water deficit on R (c) and water deficit on IR (d).

Fig. 4. Results of the GSA concerning the main effects. Sensitivity indices for the original conditions and the additional scenarios are compared in the different graphics: wind on R (a), water deficit on R (b), wind and water deficit on R (c) and water deficit on IR (d).

expected inside a polygon of a Mediterranean shrubland during a fire event. In fact, the simulations carried out included the whole range of variation of the variables describing the fuel structure and the temporal variability of the weather conditions through the year (which is adequate for

184

R. Sal6ador et al. / Ecological Modelling 136 (2001) 175–189

the GSA but probably overestimates the variability in the analysis of the scale effects). In addition, the magnitude of biases caused by the scale effects is still rather small compared to the uncertainity and the humble prediction abilities of the current fire behaviour models. Hence, the results of the analysis of the scale effects do not go against the use of fuel models, as coarse extrapolations of the spatial information of the fuel beds in spatially distributed fire behaviour models containing the Rothermel equations in their core. However, the rather stable and systematic biases found in this study suggest the possibility to apply correction factors to the parameters in order to minimise the scale effects. This would probably require the addition of estimates of variability together with the mean values of the input parameters in the fuel models. On the other hand, the results of the GSA clearly show the negligible effect of the variability of the low heat content, the particle density and the mineral content on the output values. Nevertheless, this is not an unexpected result since models based on the Rothermel equations (e.g. the BEHAVE simulator Rothermel, 1983) frequently give a constant value for these parameters (although this may partially be caused by a lack of available data for them). Besides these three variables, all other input variables have some noticeable effect over the variability of the output and they cannot be ignored if an optimal use of the model is desired. Nonetheless, although the weight of each variable may change when the scenario is changed, there are some variables that frequently account for the biggest portion of the output variance. Particularly, the height of the fuel bed and the moisture content of the dead fuel are always prominent while the slope and the proportion of dead fuels are never found among the most explicative variables. This information could be used to develop a simplified model containing only five variables (the height of the fuel bed, the surface area-to-volume ratio, the moisture content of the live and dead fuels, and the wind speed). However, although the slope never explains a big proportion of the variability of R and IR, it may have a prevalent effect on the

direction of maximum spread under low wind conditions. This will be especially important for the spatially distributed models based on the Rothermel equations. Consequently, it is not advisable to reject the slope from the list of input variables. An alternative way to reduce the complexity of the model is given by the prevalence of the main effects in some of the scenarios analysed in this study. Specifically, such prevalence suggests the possibility to create equations made up by the addition of univariate functions of the input variables. However, this approach has some serious limitations. The analytical derivation of the univariate functions may be very complex (or impossible), and although this could be overcome by fitting a function empirically to values extracted from simulations, it is frequently hard to guess the type of curve to be fitted from the scatterplots derived. Furthermore, as it has been observed in the scenario of strong wind, interaction effects may be quite important under specific conditions (where a model based on univariate functions would perform rather poorly). Next, there are some considerations that should be made in order to analyse and understand properly the results obtained in this study. First of all it has to be reminded that the GSA requires full independence of variables, which is something rarely fully accomplished. In fact, other techniques are available to estimate main effects in the case of correlated inputs (McKay, 1996), but, these may have a tremendous computational cost. The usage of fuzzy techniques has been also suggested, as their simpler rules do not include stochastic independence (Ferson and Kuhn, 1992). However, the assumptions of independence made in this study are probably quite close to the reality (it should be remembered that the effect of the fuel loading has been included as a function of the height of the vegetation). On the other hand, although the analyses of the scale effects carried out are rather robust and conclusive, more complex studies considering the usual constraints applied to ensure a continuous fire front may bring further evidences.

R. Sal6ador et al. / Ecological Modelling 136 (2001) 175–189

The analysis of wind effects also has some limitations. Since Rothermel equations did not explicitly consider the wind direction, winds were always given in the direction of the maximum slope, and this may have led to an overestimation of wind effects. Finally, if the outcomes of this study are to be considered for practical uses it should be always born in mind its geographical scope. Results obtained have been derived from Mediterranean shrublands and any extrapolation to other types of vegetation or climates has obvious risks (because the p.d.f.s of some of the input variables can be rather different). It is also crucial to acknowledge that no evidence for the predictive ability of the equations analysed on the behaviour of fire may be derived from the present study. This is an objective that can only be accomplished by comparison of predictions with data retrieved from real fires.

5. Conclusions In summary, it may be concluded that in spite of the abundant non-linearities contained in the equations studied, scale effects found have been relatively low, thus supporting the current usage of spatial averages (i.e. fuel models) in spatially distributed models based on the Rothermel equations. On the other hand, the results of the GSA clearly show the negligible effect of the variability of the low heat content, the particle density and the mineral content on the output values. However, all other input variables had some noticeable effect over the variability of the output, and they cannot be ignored if an optimal use of the model is desired. Indeed, besides fixing the values of the three variables with no noticeable effect, no other strategy to simplify the model seems to bring clear improvements to the original equations. Finally, caution is recommended if results have to be extrapolated to other types of vegetation or climate.

185

Raymond Salvador. The work of Eduard Pla was funded by a doctoral grant of the spanish Ministerio de Educacio´n y Ciencia, and by the CICYT projects AMB95-247 and AGF97-533. Appendix A. The basis of the GSA In this appendix it is shown how the output of a model can be described as a sum of functions of increasing dimensionality, with null expectation and additive variance (sections A1 and A2). The estimation of the main and total effects of the input variables using simulations is also described (section A3). Sections A1 and A2 develop the concise description of the variance decomposition given in Efron and Stein (1981). Section A3 uses the information of Parzen (1962) and Saltelli et al. (1993). A.1. The output may be described as a sum of functions of increasing dimensionality, with null expectation Let the model output be described as a function S : R n “ R. The n input variables (X1, …, Xn ) are given as random variables (r.v.) and, in consequence, S(X1, …, Xn ) is also a r.v. For a realisation of the input variables we have a specific output value S(x1, …, xn ) which may be decomposed in a sum of functions of increasing dimensionality: S(x1, …, xn ) = v+ % Ai (xi )+ % Bij (xi, xj ) i

+ %

iBjBk

iBj

Cijk (xi, xj, xk )+ ···+ H(xi, xj, …, xn ).

These functions expectations:

(a1) are,

in

turn,

sums

of

v= E(S), a fixed value (where S stands for S(X1, …, Xn )) Ai (xi )= E[S Xi = xi ]− v,

Acknowledgements This work was supported by a postdoctoral grant of the Fundacio´ Catalana per la Recerca given to

Bij (xi, xj ) = E[S Xi = xi, Xj = xj ]− E[S Xi = xi ] −E[S Xj = xj ]+ v,

(a2)

R. Sal6ador et al. / Ecological Modelling 136 (2001) 175–189

186

and so on with functions of 3, 4, …, n variables (in a way analogous to the ANOVA decomposition). By substituting Eq. (a2) in Eq. (a1) most of the terms cancel out and the following equality is obtained (Efron and Stein, 1981), S(x1, …, xn )=E[S Xi =xi, Xj =xj, …, Xn =xn ] (a3) which is completely true (thus proving the decomposition for a single realisation). Given that this equality holds for any realisation, a broader statement can be made:

E[Q1,2, … , p (X1, X2, …, Xi, …, Xp ) X1, X2, …, Xi − 1, Xi + 1, …, Xp ]= 0,

where the expectation is computed by integrating over Xi. As an example, for a second-order term: E[Bij (Xi, Xj ) Xi ] = E[E[S Xi, Xj ] Xi ]− E[E[S Xi ] Xi ]

being E[E[S Xi, Xj ] Xi ]

= v + % Ai (Xi )+ % Bij (Xi, Xj ) iBj

i

=E

iBjBk

&

&

(k"i "j) S(x1, …, xn )

Cijk (Xi, Xj, Xk ) +···

+ H(Xi, Xj, …, Xn ),

(a4) =

where:

& &

&

Xi

n

5 f(xk ) dxk



k"i"j

5 f(xk ) dxk

(k"i" j) S(x1, …, xn )



k"i"j

j

f(xj ) dxj = E[S Xi ],

Ai (Xi )=E[S Xi ]− v, Bij (Xi, Xj )= E[S Xi, Xj ]− E[S Xi ]− E[S Xj ] + v,

E[E[S Xi ] Xi ]= E[S Xi ] (E[S Xi ]

and so on (in other words conditional expectations become r.v.). Next, since expectations of conditional expectations are unconditional expectations (e.g. E[E[S Xi ]]= v, E[E[S Xi, Xj ]] = v, E[E[S Xi, Xj, Xk ]]= v, etc.), all functions in Eq. (a4) have null expectations: E[Ai (Xi )]=E[E[S Xi ]− v]= E[E[S Xi ]]− v =0

is already an univariate functio n of Xi ), E[E[S Xj ] Xi ]= E[E[S Xj ]]= v (Xi is not a variable of the function E[S Xj ]) , and E[v Xi]= v

(v is a constant).

A.2. The 6ariance of the output may be gi6en as the sum of the 6ariances of functions included in the decomposition

E[Bij (Xi, Xj )] = E[E[S Xi, Xj ]]−E[E[S Xi ]]− E[E[S Xj ]] + E[v] =v− v− v +v = 0 etc.

(a6)

− E[E[S Xj ] Xi ]+ E[v Xi ]= 0

S(X1, …, Xn )

+ %

evaluated over any of its own variables, is null. Specifically, it can be proved that:

(a5)

As a consequence of Eq. (a4) and Eq. (a5), the expectation of each term [say Q1,2,…,p (X1, X2, …, Xp ))] in the decomposition,

When considering any pair of functions in the decomposition (Eq. (a4)) there is always (at least) one input variable that belongs to one of the two functions only. Thus, taking into account Eq. (a5), the covariance between two functions P and Q, which share the variables Xi, …, Xp, is given by

R. Sal6ador et al. / Ecological Modelling 136 (2001) 175–189

Cov[Pi,…,m (Xi, …, Xp, Xr, …, Xm )

The decomposition permits assessing the effect on the model output variability of each input variable through the variances of the main effect, and of any group of input variables through higher order variances (see Section 2.4 of the main text).

, Qi,…,n (Xi, …, Xp, Xs, …, Xn )] =E[Pi,…,m (Xi, …, Xp, Xr, …, Xm )

& && && &

× Qi,…,n (Xi, …, Xp, Xs, …, Xn )]

=

i,… p

Pi,…,m (xi, …, xp, xr, …, xm )

r,… m s,… n

5

×Qi,…,n (xi, …, xp, xs, …, xn )

f(x)

i,…,p,r,…,m;s,…,n

& && && & dx

=

i,… p

!&

Qi,…,n (xi, …, xp, xs, …, xn ) f(xn ) dxn

"

1. The first order variance (the main effect):

n

f(x) dx.

(a7)

i,…,p,r,…,m,s,…,n − 1

&

A.3. The relati6e weight of the main and total effects may be estimated by Monte Carlo integrals As pointed in Section 2.4 of the main text, due to the high number of terms appearing in the variance decomposition formula, only two types of variances are usually quantified:

r,… m s,… n −

1 Pi,…,m (xi, …, xp, xr, …, xm )

5

187

Var[Ai (Xi )]= E[Ai (Xi )2] = E[(E[S Xi ])2]− 2vE[E[S Xi ]]+ v 2

As pointed out in Eq. (a6) one has

= E[(E[S Xi ])2]− v 2.

(a9)

Qi,…,n (xi, …, xp, xs, …, xn ) f(xn ) dxn

n

+··· +Var[H(Xi, Xj, …, Xn )]

An alternative interpretation for Var[Ai (Xi )] comes from the fact that it is equal to Var[E(S Xi )], which, in turn, is equivalent to Var[S]− E[Var(S Xi )] (see Parzen, 1962). Thus, the main effect can be thought as the average reduction of the model output variance that would be achieved by fixing Xi to its true value xi. The average is made by varying Xi over all possible values xi (as xi is unknown). 2. The sum of the variances of all the terms in the decomposition that include the variable Xi (the total effect): This is equivalent to computing E[Var[S X1, …, Xi − 1, Xi + 1, …, Xn ]], i.e. the average variance that would remain as long as Xi is allowed to vary.

= %E[Ai (Xi )2]+ % E[Bij (Xi, Xj )2]

E[Var[S X1, …, Xi − 1, Xi + 1, …, Xn ]]

= E[Qi,…,n (Xi, …, Xp, Xs, …, Xn ) Xi, …, Xn − 1] = 0, thus leading to a null covariance, i.e. orthogonality between any pair of variables. Consequently, the variance of S(X1, …, Xn ) can be split into the sum of the variances of the functions in Eq. (4): Var[S(X1, …, Xn )] =% Var[Ai (Xi )]+ % Var[Bij (Xi, Xj )] iBj

i

+ %

iBjBk

Var[Cijk (Xi, Xj, Xk )]

iBj

i

+ %

iBjBk

= E[E[S 2 X1, …, Xi − 1, Xi + 1, …, Xn ] 2

E[Cijk (Xi, Xj, Xk ) ]

−(E[S X1, …, Xi − 1, Xi + 1, …, Xn ])2]

+ ··· +E[H(Xi, Xj, …, Xn )2] (due to the nullity of expectations).

(a8)

=E[S 2]− E[(E[S X1, …, Xi − 1, Xi + 1, …, Xn ])2]. (a10)

188

R. Sal6ador et al. / Ecological Modelling 136 (2001) 175–189

The joint computation of main and total effects requires evaluating v 2, E[S 2] (which are both computed once), E[(E[S Xi ])2] and E[(E[S X1, …, Xi − 1, Xi + 1, …, Xn ])2] (which are computed for each variable). As said in Section 2.4, due to their analytical complexity, numerical methods (Monte Carlo integrals) are required to estimate them. Since both E[(E[S Xi ])2] and E[(E[S X1, …, Xi − 1, Xi + 1, …, Xn ])2] are conditional expectations their estimation from simulations have a high computational cost (a whole set of simulations has to be carried out for each simulated value of Xi ). In order to minimise this problem Saltelli et al. (1993) give an improved computational scheme.

References Albini, F.A., 1976. Estimating wildfire behavior and effects. Tech. Rep. INT-30. Forest Service, Intermountain Forest and Range Experiment Station, Ogden (UT), pp. 91. Albini, F.A., 1985. A model for fire spread in wildland fuels by radiation. Combustion Science and Technology 42, 229– 258. Anderson, H.E., 1982. Aids to determining fuel models for estimating fire behavior. Tech. Rep. INT-30. Forest Service, Intermountain Forest and Range Experiment Station, Ogden (UT), 22 pp. Andrews, P.L., 1986. BEHAVE: Fire behaviour prediction and fuel modelling system. Burn Subsystem. Part 1. PMS-4392. USA. Dept. of Interior, Boise (Idaho), 130 pp. Beer, T., 1991. Bushfire rate of spread forecasting: deterministic and statistical approaches to fire modelling. Journal of Forecasting 10, 301–307. Burgan, R.E., 1987. Concepts and interpreted examples in advanced fuel modelling. Tech. Rep. INT-238. Forest Service, Rocky Mountain Research Station, Ogden (UT), 40 pp. Catchpole, E.A., Catchpole, W.R., 1991. Modelling moisture damping for fire spread in a mixture of live and dead fuels. International Journal of Wildland Fire 1, 101–106. Cheney, N.P., Gould, J.S., Catchpole, W.R., 1993. The influence of fuel, weather, and fire shape variables on firespread in grasslands. International Journal of Wildland Fire 3, 31 – 44. Cukier, R.I., Schaibly, J.H., Shuler, K.E., 1975. Study of the sensitivity of coupled reaction systems to uncertainties in rate coefficients. III. Analysis of the approximations. The Journal of Chemical Physics 63, 1140–1149. Efron, B., Stein, C., 1981. The jackknife estimate of variance. The Annals of Statistics 9, 586–596.

Elvira, L.M., Hernando, C., 1989. Inflamabilidad y energı´a de las especies de sotobosque, INIA, Madrid, 73 pp. Ferson, S., Kuhn, R., 1992. Propagating uncertainty in ecological risk analysis using interval and fuzzy arithmetic. In: Zannetti, P. (Ed.), Computer Techniques in Environmental Studies IV. Elsevier, London, pp. 387– 401. Finney, M.A., 1998. FARSITE: Fire area simulator — model development and evaluation. Res. ap. RMRS-RP-4. Forest Service, Rocky Mountain Research Station, Ogden (UT), 47 pp. Frandsen, W.H., 1971. Fire spread through porous fuels through the conservation of energy. Combustion and Flame 16, 9 – 16. Green, D.G., Gill, A.M., Noble, I.R., 1983. Shapes of similar fires in discrete fuels. Ecological Modelling 20, 21 – 32. Haje´k, J., 1968. Asymptotic normality of simple linear rank statistics under alternatives. The Annals of Mathematical Statistics 39, 325– 346. Hamby, D.M., 1994. A review of techniques for parameter sensitivity analysis of environmental models. Environmental Modelling and Assessment 32, 135– 154. Helton, J.C., 1993. Uncertainty and sensitivity analysis techniques for use in performance assessment for radioactive waste disposal. Reliability Engineering and System Safety 42, 327– 367. Homma, T., Saltelli, A., 1996. Importance measures in global sensitivity analysis of model output. Reliability Engineering and System Safety 52, 1 – 17. Lopes, A.M.G., Cruz, M.G., Viegas, D.X., 1998. FIRESTATION — An integrated system for the simulation of wind flow and fire spread over complex topography. III International Conference on Forest Fire Research — 14th Conference on Fire and Forest Meteorology, Luso (Portugal) 16/20 November, Vol. I, pp. 741– 754. Martı´nez-Milla´n, J., Saura, S., 1998. CARDIN 3.2. Forest fires spread and fighting simulation system. III International Conference on Forest Fire Research — 14th Conference on Fire and Forest Meteorology, Luso (Portugal) 16 – 20 November, Vol. I, pp. 907– 914. McAlpine, R.S., 1995. Testing the effect of fuel consumption on fire spread rate. International Journal of Wildland Fire 5, 143– 152. McArthur, A.G., 1966. Weather and grassland fire behaviour. Australian Forestry and Timber Bureau Leaflet 100, Camberra, AFTB. McKay, M.D., 1996. Variance-Based Methods for Assessing Uncertainty Importance in NUREG-1150 Analyses. LAUR-96-2695, Los Alamos National Laboratory, Los Alamos (NM), pp. 1 – 27. Panov, P., 1994. Analysis of pyric properties of Mediterranean type vegetation. MSc thesis. Mediterranean Agronomic Institute of Chania, Chania. Papio´, C., 1990. Ecologia del foc i regeneracio´ en garrigues i pinedes mediterranies. Phd thesis. Universitat Auto`noma de Barcelona, Barcelona, 314 pp. Parzen, E., 1962. Stochastic processes. Holden day, San Francisco, 324 pp.

R. Sal6ador et al. / Ecological Modelling 136 (2001) 175–189 Perry, G.L.W., 1998. Current approaches to modelling the spread of wildland fire: a review. Progress in Physical Geography 22, 222– 245. Pla, E., Roda`, F., 1999. Aproximacio´ a la dina`mica successional de combstible en brolles mediterra`nies. Orsis 14, 79– 103. Pook, E.W., Gill, A.M., 1993. Variation of live and dead fine fuel moisture in Pinus radiata plantations of the Australian Capital territory. International Journal of Wildland Fire 3, 155– 168. Rothermel, R.C., 1972. A mathematical model for predicting fire spread in wildland fuels. Res. Pap. INT-115. Forest Service, Intermountain Forest and Range Experiment Station, Ogden (UT), 40 pp. Rothermel, R.C., 1983. How to predict the spread and intensity of forest and range fires. PMS-436-1, USA. Department of interior, Boise (Idaho), 161 pp. Salazar, L.A., 1985. Sensitivity of fire behavior simulations to fuel model variations. Res. Pap. PSW-178. Forest Service, Pacific Southwest Forest and Range Experiment Station, Berkeley, (CA), 11 pp. Saltelli, A., Andres, T.H., Homma, T., 1993. Sensitivity analysis of the model output. An investigation of new techniques. Computational statistics and data analysis 15, 211– 238. Saltelli, A., Chan, K., Scott, M., 2000. Mathematical and statistical methods for sensitivity analysis. Wiley, New

.

189

York. Saltelli, A., Tarantola, S., Chan, K.P.-S., 1999. A quantitative model-independent method for global sensitivity analysis of model output. Technometrics 41, 39 – 56. Sobol’, I.M., 1990. Sensitivity estimates for nonlinear mathematical models. Matematicheskoe Modelirovanie, 2: 112– 118, translated as: Sobol’, I.M., 1993. Sensitivity analysis for nonlinear mathematical models, Mathematical Modelling and Computational Experiment 1, 407– 414. Stocks, B.J., Lawson, B.D., Alexander, M.E., van Wagner, C.E., McAlpine, R.S., Lynham, T.J., Dube´, D.E., 1989. The Canadian forest fire danger rating system: an overview. Forestry Chronicle 65, 450– 457. Viegas, D.X., Pin˜ol, J., Viegas, M.T., Ogaya, R. Estimating live fine moisture content using meteorologically-based indices, International Journal of Wildland Fire In press. Weber, R.O., 1989. Thermal theory for determining the burning velocity of a laminar flame using the inflection point in the temperature profile. Combustion Science Technology 64, 135– 139. Wilson, R.A., Jr., 1990. Reexamination of Rothermel’s fire spread equations in no-wind and no-slope conditions. Res. Pap. INT-434. Forest Service, Intermountain Forest and Range Experiment Station, Ogden (UT), 13 pp. Wu, Y., Sklar, F.H., Gopu, K., Rutchey, K., 1996. Fire simulations in the everglades landscape using parallel programming. Ecological Modelling 93, 113– 124.