Addressing Wildfire Risk in a Landscape-Level Scheduling Model: An Application in Portugal

The paper presents and discusses research aiming at the development of a forested landscape management scheduling model that may address the risk of wildfires. A general indicator is built from wildfire occurrence and damage probabilities to assess stand-level resistance to wildfires. This indicator is developed to further address the specificity of each stand configuration (e.g., shape and size) and spatial context (neighboring stands characteristics). The usefulness of the development of such an indicator is tested within a mixed integer programming (MIP) approach to find the location and timing of management options (e.g., fuel treatment, thinning, clearcut) that may maximize the forested landscape expected net revenues. The Leiria National Forest, a Portuguese forest in central Portugal, was used as a case study. Results suggest that the proposed approach may help integrate wildfire risk in forested landscape management planning and assess its impact on the optimal plan. Results further show that prescriptions that include fuel treatments are often chosen over others that do not include them, thus highlighting the importance of wildfire management efforts. Finally, they provide interesting insights about the role of thinnings and fuel treatment in mitigating wildfire risk.

wildfire risk, forest management scheduling, mixed integer programming

Wildfires are one of the main threats for forests in the Mediterranean countries, namely in Portugal (Alexandrian et al. 2000, Pereira et al. 2006). Large-scale forest fires throughout Mediterranean countries (e.g., Portugal, Spain, Italy, and Greece) have substantially increased during the last few decades (Velez 2006). The need to address wildfire risk in forest management planning is evident and yet fire and forest management are currently carried out mostly independently of each other in these countries (Borges and Uva 2006). Forest management planning addresses landscapewide goals and encompasses the scheduling of management options in each stand in the landscape mosaic. Stand management options are tied together by the landscapewide goals, e.g., even timber flows and ending inventory concerns (Borges et al. 1999, Hoganson et al. 2008), hence the importance of landscape-level models to develop stand management schedules.

The losses caused by wildfires are one of the major sources of uncertainty when projecting timber supply in forest management planning. They thus have a substantial impact on landscapewide objectives. Wildfires are a spatial phenomenon that may propagate from stand to stand due to several factors such as adjacency, wind, slope, and altitude, among others. Therefore, wildfire size, intensity, and postfire mortality in any given stand will be influenced by its characteristics as well as by the features of adjacent stands and their relative position. Thus, the much-needed incorporation of wildfire risk in forest management planning does require stand-level information (e.g., inventory and environmental features), as well as geographical information to characterize each stand spatial context (e.g., Agee et al. 2000, Acuna et al. 2010, Konoshima et al. 2010, Gonzáélez and Pukkala 2011). In this context, Kalabokidis et al. (2002) developed a conceptual framework to compute fire potential at landscape level, identifying external and internal factors that may influence it.

The development of decision trees (Cohan et al. 1983), mathematical programming models (e.g., Reed and Errico 1986, Gassmann 1989, Lilieholm et al. 1991, Boychuk and Martell 1996, Martins et al. 2005, Constantino et al. 2008, Acuna et al. 2010), simulation (Van Wagner 1983, Peter and Nelsson 2005), combination of Monte Carlo simulation and mathematical programming (Armstrong 2004), Markov chains (Zhou and Buongiorno 2006, Campbell and Dewhurst 2007), or heuristics (e.g., Hoganson and Rose 1987, Borges and Hoganson 1999, Falcão and Borges 2002, Caro et al. 2003, Gonzáélez et al. 2005) have proved to provide a framework that facilitates the integration of considerations about wildfire risk in forested landscape management planning models. Van Wagner (1983) pioneered research of modeling approaches to assess the impact of the proportion of burned area on timber supply. Reed and Errico (1986) developed a linear programming harvest scheduling model that took as input the fraction of burned forest in each year, considering the revenue losses due to wildfires. Hoganson and Rose (1987) developed a multistage stochastic linear programming harvest scheduling approach where the short-term stage was considered as risk free and the long-term stage was associated with several probabilistic scenarios. Gassmann (1989) developed a nonspatial stochastic programming model with several stages to optimize harvest scheduling in the presence of fire. Boychuk and Martell (1996) also developed a stochastic programming model where the first stage involved decisions made in a risk-free framework and where, at the beginning of each period, the forest is characterized by the area in each age class. Gonzáélez et al. (2005) integrated stand-level wildfire risk models into forestwide management planning models to optimize landscape metrics reflecting the spatial distribution of forest fuels and fire hazards. Despite the difficulty of ad- dressing forested landscape management planning with dynamic programming, due to the dimensionality curse, some authors developed formulations to overcome this constraint (e.g., Borges and Hoganson 2000, Spring and Kennedy 2005, Hoganson et al. 2008, Forsell et al. 2001, Konoshima et al. 2010).

Since forest management impacts the fuel load as well as its spatial distribution over the landscape, the most promising modeling approaches consider wildfire risk as endogenous to the management planning problem. Palma et al. (2007) and Acuna et al. (2010) developed approaches to estimate the contribution of each stand to reduce wildfire spread and losses. The latter developed a mixed integer programming model (MIP) to extend the approach by Reed and Errico (1986) and specify where and when to harvest blocks of stands so that clearcuts may lead to the reduction of wildfire risk and to maximize financial returns. Wei et al. (2008) expanded MIP models to locate fuel treatments in management plans based on ignition risk, probabilities of spread, fire intensity, and values at risk. Wei (2012) further demonstrated the potential of MIP to examine the location of fuel treatments to improve the efficiency of controlling future catastrophic fires. Moreover, several authors addressed fuel treatment scheduling to further integrate wildfire considerations in forested landscape management planning (e.g., Agee et al. 2000, Finney 2000, Gonzáélez et al. 2005, Palma et al. 2007, Kim et al. 2008, Konoshima et al. 2010, Acuna et al. 2010). Other authors have explicitly incorporated the use of wildfire spread simulators in forest planning models (Johnson et al. 1998, Bettinger 2009, Young-Hwan et al. 2009). However, in most cases, the use of wildfire spread simulators in forest planning is somehow limited as they require input data that are hardly available in a management planning framework characterized by extended temporal horizons (i.e., weather conditions, fuel moisture).

This approach to integrate wildfire risk in landscape-level forest management planning involved the development of an indicator to characterize a stand resistance to fire. This indicator builds from stand-level wildfire occurrence and damage probabilities (Garcia-Gonzalo et al. 2012, Marques et al. 2012). It takes further into account stand geometric (e.g., shape and size) and topological (e.g., spatial context) features so that it may address wildfire spread concerns. The emphasis was in adjusting a stand intrinsic resistance to wildfire according to its neighboring relations and the features of adjacent stands. The resistance indicator was incorporated in a harvest scheduling mixed integer model. This model seeks the combination of management alternatives that maximizes the forest value while targeting an aggregated value of wildfire resistance in each planning period. After presenting the modeling approach, we discuss results of an application to Leiria National Forest, a public forest in central Portugal.

The Modeling Approach

The forest management scheduling problem will be designed as a linear programming Model I (Johnson and Sheurman 1977). As such, its solution encompasses the assignment of a prescription to each stand i, i^aI = {1, …, I} in the forest so that landscapewide goals may be approximated. The landscape-level model should thus reflect spatial and temporal relations between stand-level decisions. In this research, the modeling approach recognizes the impact of decisions, made in a stand, on wildfire occurrence and severity probabilities in the set of its neighboring stands (e.g., the implementation of fuel treatments in a stand may help constrain the spread of fires to adjacent stands) (e.g., Finney 2000, Gonzáélez et al. 2005, Konoshima et al. 2010, Acuna et al. 2010). This entailed the development of a periodic stand-level wildfire resistance indicator. A periodic landscapewide wildfire resistance indicator was computed further as the weighted average of the stand-level indicators where the weights correspond to the stand area.

Adjusted Resistance Indicator

The wildfire resistance indicator of any stand i was thus built to reflect the environmental and biometrical features (e.g., stand density, shrubs biomass) of both stand i and its adjacent stands. The former contribute to the stand-i-specific resistance, while the latter are combined with the risk of wildfire spread to stand i to compute an adjusted resistance indicator. The specific resistance of stand i in period t is termed R_it. It reflects the proportion of trees in stand i that are prone to survive a wildfire occurring in planning period t, according to wildfire occurrence and postfire mortality probability models (Garcia-Gonzalo et al. 2011, Marques et al. 2012).

The adjusted resistance level of stand i in period t is termed RA_it. It incorporates the wildfire resistance levels of its adjacent stands and the risk of fire spread. It is defined as the solution of the following equations  

with

R_it = specific resistance of stand i in period t; R_it^a [0,1]; (1 − w_i) = weight to reflect the impact of neighboring stands on the wildfire resistance of stand i; (1 − w_i) [0,1]; w_i depends on parameters that reflect the size and shape of the stand i; we assume w_i > 0.

α_is = parameter reflecting the likelihood of a fire that occurs in stand s to spread to stand i; 0 ≤ α_is ≤ 1 and , where ν {i} is the set of neighbors of stand i. With these assumptions, it can be shown Equations 1 have exactly one solution (Ferreira 2011) (see Appendix).

The impact of neighboring stands, (1 − w_i), on wildfire resistance of stand i is assumed to decrease (increase) with the area (perimeter) of the latter, while the likelihood of a wildfire that occurs in stand s to spread to stand i, α_is, is assumed to increase with the length of the edge between the two stands. The difference RA_st − R_it thus corresponds to the adjustment of the specific wildfire resistance of stand i(R_it) might undergo to address the impact of its spatial context. For example, the wildfire resistance of stand i is prone to increase if the adjusted resistance of a neighboring stand s is greater than the specific resistance of stand i. A stand-specific wildfire resistance indicator might also be adjusted by considering the specific wildfire resistance of its neighbors s, rather than their adjusted wildfire resistance  

Nevertheless, in this case, the spatial context is interpreted more limitedly, as it is assumed that the wildfire resistance of a stand is impacted only by its adjacent stands rather than by a wider neighborhood. Thus, in this research, we assumed that a stand-specific wildfire resistance indicator is to be adjusted according to the adjusted wildfire resistance of its neighbors. In both cases (Equations 1 and 2), the wildfire resistance indicator value ranges from 0 and 1 (Ferreira 2011) (for more details, see Appendix).

The MIP Approach

The approach to address wildfire risk concerns in forest management scheduling integrates the stand-adjusted resistance index into an MIP model that may be described as follows  

 

 

 

 

 

 

 

 

 

where

Sets and indices:

• i identifies the stand, with iI={1, …, I} and I corresponds to the total number of stands; k identifies a prescription that may be assigned to stand i, k = {1, …, K_i}; t identifies the planning period, with t = {1, …, T}; and T corresponds to the number of planning periods.

• {i} = set of stand i's neighbors.

Parameters:

• c_ik = net present value, euros per hectare, resulting from the use of prescription k to manage stand i. It includes the value of ending inventory and it is adjusted to reflect the impact of wildfires as described in the Case Study section, under Results.

• A_i = area of stand i, in hectares; A_F = total forest area, in hectares, A_F = ∑_{i = 1}^IA_i; r_ikt = specific wildfire resistance of stand i in period t if it is managed according to the prescription k; RES_t = landscapewide wildfire resistance in period t; RES_t [0,1]; (1 − w_i) = as in Equation 1; α_is = as in Equation 1; ν_ikt = volume per hectare harvested in period t if prescription k is assigned to stand i; F_ik = age of stand i at the end of the planning horizon if prescription k is assigned to it; f = minimum required average age of the forest at the end of the planning horizon; λ = maximum allowed variation of volume harvested in two consecutive periods.

Variables:  

• R_it = specific resistance of stand i, in period t.

• RA_it = adjusted wildfire resistance of stand i, in period t.

• V_t = volume harvested in period t.

Equation 3 expresses the aim of maximizing the forest value. Equations 4 ensure that one and only one alternative is assigned to each stand.

Equation 5 defines the volume harvested in each planning period, while Equations 6 define the volume wood flow constraints. The latter state that the volume harvested in consecutive periods should not fluctuate over more than λ. The concern about the state of the forest at the end of planning horizon is reflected by Equation 7 that ensures that the average age of the ending inventory is greater or equal than a minimum value f.

Equations 8 compute the specific wildfire resistance of each stand in each period, according to the prescription assigned to it. Equations 9 compute the adjusted wildfire resistance of each stand, in each period, as described in Equation 1. To address wildfire risk concerns, harvest scheduling is constrained by setting a minimum level of forestwide adjusted wildfire resistance, computed as the weighted average of its stands' adjusted wildfire resistance, where the weights correspond to the stand area (Equation 10).

The MIP model is deterministic and yet it may address wildfire fire risk by taking advantage of the information provided by the adjusted wildfire resistance. This indicator sets a protection level target and is used further to adjust the cost coefficients of the decision variables. In the The Adjusted Cost Coefficient in the MIP Model section, we show, through a case study, how this adjustment is made.

Case Study

Leiria National Forest

To test the modeling approach we will consider Leiria National Forest (LNF), a public even-aged maritime pine forest in Portugal. LNF extends over an area of approximately 11,000 hectares, of which about 8,600 have been allocated to timber production. This area has been classified into 393 stands very regular in shape. A typical pine prescription includes a plantation of about 2,000 seedling per ha, a mandatory precommercial thinning at 15 years of age that leaves about 1,500 trees per ha, thinnings every 5 years between 20 and 50 years of age, and rotation ages that may extend from 50 to 80 years (Falcão and Borges 2001). In the instance considered in this case study, stands with age lower than 20 years extend over 13% of the forest area, while stands with age between 20 and 50 years extend over 65% of LNF productive area. The remaining 22% of the forest area is occupied by mature stands with age greater than 50 years.

The planting cost in LNF includes a fixed and a variable component (Table 1). The former includes the soil preparation and the fuel treatment costs. The latter includes the seedling costs. In LNF, timber is typically sold before harvesting and stumpage prices are a function of age (Table 2). After a wildfire, some timber may be salvaged and sold. Generally, salvage prices correspond to 75% of the original prices. In the case of stand ages for which timber prices were not available, these were computed by linear interpolation of available prices (Table 1). A 4% rate was used to discount costs and revenues.

Stumpage prices used for forest management scheduling in Leiria National Forest.

Table 1.

Stumpage prices used for forest management scheduling in Leiria National Forest.

Values in Borges and Falcão (1999) were updated according to pine stumpage prices variation since 1998. Source: Borges and Falcão (1999).

Operational costs.

Table 2.

Operational costs.

Source: CAOF's Database of ANEFA—The National Association of Forest, Agriculture and Environment Enterprises.

Maritime pine growth was estimated using a simplified version of the stand-level growth-and-yield model DUNAS (Falcão 1998). This model is based on difference equations and characterizes the existence of a stand at each age, according to the basal area, using the McDill and Amateis model (McDill and Amateis 1992). Understory growth was estimated according to a model developed by Botequim et al. (2009). According to this model, the understory biomass level depends on its age and on the percentage of sprouters.

To estimate the percentage of live trees in stand i in period t, if it is managed according to prescription k, this research used wildfire occurrence and damage models developed for pine stands in Portugal (Garcia-Gonzalo et al. 2011, Marques et al. 2012). The wildfire occurrence model estimates the probability of wildfire occurrence in a stand, according to its state (e.g., age, basal area, number of trees, biomass level). The damage model estimates the proportion of dead trees whenever mortality occurs as a consequence of a wildfire.

Specific Resistance Indicator

Forest structure and composition evolve in time and a stand wildfire resistance changes over the planning horizon. The indicator to estimate the specific wildfire resistance of a stand in LNF, in a planning period, should consider the stand state in this period, as well as the annual probabilities of wildfire occurrence, the probability of a wildfire to cause mortality, and the proportion of dead trees if mortality indeed occurs (12)  

where

• Pf_u^ikt = probability of wildfire occurrence in stand i in year u of period t if it is managed according to prescription k.

• Pmort_u^ikt = probability of mortality to occur if there is a wildfire in stand i, being managed according to prescription k, in year u of period t.

• Pam_u^ikt = proportion of dead trees in stand i, being managed according to prescription k, in year u of period t if mortality occurs as a consequence of a wildfire.

The computation of a LNF stand-specific wildfire resistance in period t, when it is managed according to prescription k, r_ikt, assumes independence of wildfire occurrence in different years. The proportion of surviving trees, in each year of a period t, builds from the number of live trees in the previous year. The indicator thus computes the proportion of surviving trees at the end of period t. r_ikt may also be interpreted as the expected percentage of surviving trees in stand i, in period t, if prescription k is adopted.

The Adjusted Wildfire Resistance Indicator

The specific wildfire resistance indicator is adjusted using two parameters (Equation 9). The first, (1 − w_i), reflects the overall impact of neighboring stands on the wildfire resistance of stand i, while the second, α_is =, explains the contribution of each individual neighboring stand to that impact. The assignment of values to w_i, iI is based on the relation between the area and the perimeter of stand i, as it reflects its shape and size. The greater the perimeter, the greater will be the impact of neighboring stands on the wildfire resistance of stand i. Conversely, the greater the area, the lower will be this impact. This research uses the ratio to compute w_i. This ratio has been introduced (e.g., Thomas 1979, Ohman and Lamas 2005) to address spatial considerations in forest management. It reaches its maximum value, one, when the stand i is a perfect circle. Specifically, the weight w_i will be calculated according to the expression w_i = , with 0 < θ ≤ 1. If stands have the same area, the weight w_i will decrease with the perimeter, while in the case of stands with the same perimeter, the weight will increase with the area. To ensure further that different weights are assigned to stands that have the same value for the ratio , but that have different areas, the computation of weights w_i considers the parameter θ = , where ρ_i ≥ 0 may be another scale-related parameter (e.g., the stands average area). For example, in the case of stands with very large areas, the parameter θ is closer to 1, while in the case of small stands it approaches 0.

To assess the impact of a neighboring stand s on the wildfire resistance of stand i, the parameter α_is takes into account the percentage of border that stand i shares with stand s, the slope, the relative position of stands i and s, and the existence of barriers between the stands (e.g., roads). Accordingly, it may be computed as  

where

• f_is = fraction of the border of stand i that is shared with stand s; f_is [0,1] and .

• m_is = parameter that reflects the likelihood of wildfire spread from stand s to stand i; it takes into account the relative position of i and s and the existence of barriers between the two stands; m_is [0,1].  

  

In this case study, the factor m_is is decomposed into two components—u and p—i.e., m_is = u_isp_is and α_is is thus computed as  

where

u_is = likelihood of fire spread from stand s to stand i, based on the relative position of i and s; u_is [0,1].

p_is = likelihood of fire spread from stand s to stand i, based on the existence of barriers between the two stands; p_is [0,1].

The relative altitude of both stands s and i influences the possibility of wildfire spread from s to i. If stand s is at an altitude lower than stand i, the more likely it is that a fire that crosses s spreads to i, as, in general, wildfire spreads uphill (e.g., Van Wagner 1988). In this case, the impact of s on the wildfire resistance of i is greater. Wildfire spread is also influenced by slope. In general, the rate of wildfire spread increases with the slope (e.g., Van Wagner 1988, Gonzáélez and Pukkala 2011), so if the stand s is at a lower elevation than i, the higher the slope the greater the likelihood of wildfire spread from s to i. Conversely, if s is at a higher elevation than i, the higher the slope the lower the likelihood of wildfire spread from s to i.

The aspect of neighboring stands also influences the spread of fire. If the aspect of neighboring stands s matches the direction from where wind blows and if s is at a lower elevation than i, then the higher is the likelihood of fire spread from s to i. Finally, the existence of barriers hinders wildfire spread and leads to a smaller influence of neighboring stands on the resistance of the stand i.

LNF is at a low altitude and the slope of its stands varies between 0 and 26 degrees. In general, the wind is toward southeast. The estimation of the values of parameters needed to compute the stands' adjusted wildfire resistance took advantage of expert knowledge available (Table 3) and met two conditions. Firstly, u_is is larger when s is at an altitude lower than i, the slope is high, and the s aspect is northwest. Secondly, u_is is smaller when s is at an altitude higher than i, the slope is high, and the s aspect is northwest.

Fire spread likelihood u considering the relative position between i and s, the orientation and the slope of s.

Table 3.

Fire spread likelihood u considering the relative position between i and s, the orientation and the slope of s.

In LNF, the road network includes roads perpendicular to the coast that are 10 m wide and roads that run parallel to the coast and are 5 m wide. The likelihood of wildfire spread from a stand to another thus depends on the type of road (if any) between them (Table 4). The weight p_is is equal to 1 if there is no road between the stands and it decreases with road width.

Weight to fire spread likelihood, considering the type of road between s and i.

Table 4.

Weight to fire spread likelihood, considering the type of road between s and i.

The LNF area that is allocated to timber production shares borders with other land uses. Thus, this research used normalized values of f_is (f_{is_n}) to consider only neighboring land units within that area, f_{is_n} = , with s being a stand within that area. Thus, it is assumed implicitly that the contribution from other land units to the resistance of stand i is similar to the contribution of neighboring stands in the LNF area that is allocated to timber production. An alternative would be to assume a fixed value for the adjusted resistance of those land units.

The Adjusted Cost Coefficient in the MIP Model

The MIP approach is deterministic. Nevertheless, it may effectively address wildfire risk concerns by targeting forestwide wildfire resistance conditions. These concerns are addressed further by adjusting the cost coefficients of the decision variables in the objective function. Decision variables consist of prescriptions, i.e., sequences of management options, such as clearcutting (CCP), thinning (TC), cleaning (C), or no intervention (NI) (Table 5) to be implemented in each stand of LNF in each planning period. All prescriptions involve at most a clearcut.

Possible prescriptions, considering a planning horizon of four periods.

Table 5.

Possible prescriptions, considering a planning horizon of four periods.

CCP—clearcutting involving also a fuel treatment and a new plantation; TC—thinning, involving also a fuel treatment; C—fuel treatment; NI—no intervention.

The value of the cost coefficient c_ik of a prescription k is typically computed by adding the net present value of its management options (NPV_{ikt_int}, with int = {CCP, NI, C, TC}) and the present value of ending inventory. These values should be interpreted as expected values, as they depend on the specific wildfire resistance indicators r_ikt, which in turn depend on the probabilities of wildfire occurrence, mortality, and proportion of dead trees. To reflect the impact of wildfires, the net present value of management options is adjusted so that it includes the present value of revenues resulting from the sale of salvaged burned wood (SAL_{ikt_int})  

Moreover, the net present value of other management options is also adjusted according to the probability of wildfire occurrence, the probability of mortality to occur after a wildfire, and the proportion of dead trees if mortality indeed occurs. This information is encapsulated in the specific wildfire resistance indicators (r_ikt).

Thus, in the case of thinning and clearcut management options, the net present value resulting from timber sales is computed as  

 

 

where

r_ikt = specific wildfire resistance of stand i, in period t, when managed according to prescription k.

v_ikt = volume harvested from stand i in period t, when managed according to prescription k.

P1_ikt = present value of stumpage price (€/m³).

P2_ikt = present value of salvage price (€/m³).

PC_t = present value of planting cost.

Equations 16 and 17 compute the net present value resulting from the sale of live trees, while Equation 18 computes the net present value of volume salvaged (1 − r_ikt)v_ikt.

Since in LNF timber is typically sold as stumpage, it is assumed that the cost of fuel treatments associated to harvesting operations, e.g., like CCP or TC, is supported by the buyer. Moreover, no fuel treatments are implemented whenever a wildfire occurs. Thus, in this research, the cost of a fuel treatment scheduled to occur every d₁ years is computed as follows: If d is the length of a planning period, whenever a wildfire occurs in the last d₁ years of period t, no fuel treatment is required. Assuming d > d₁, the probability of no wildfire occurrence in last d₁ years of period t is  

where

Pf_u^ikt = probability of wildfire occurrence in stand i, in year u of period t, when managed according to prescription k.

Therefore, the cost of a fuel treatment in period t will be  

where

ĉ_t = present value for fuel treatment cost.

The values for the probabilities take into consideration the possibility a fire in stand i to originate elsewhere. One of the limitations of the model is that those values depend only on the prescriptions assigned to stand i.

Results

The MIP model considered a four 10-year periods planning horizon and up to 40 prescriptions per stand in the area allocated to timber production in LNF. The minimum average age at the end of the planning horizon was set to 40 years (f = 40), while the timber flow constraints considered a maximum fluctuation of 25% (λ = 0.25) between consecutive periods. The minimum landscapewide wildfire resistance levels were set to range from RES_t = 0.7 to RES_t = 0.89, with t = 1, …, 4. The model encompassed 1,977 constraints and 4,069 variables, of which 2,496 are binary (Table 6).

Report of the MIP model size.

Table 6.

Report of the MIP model size.

The MIP model was implemented with the algebraic language ILOG OPL 6.3 and it was solved by CPLEX 12.1, in a computer-i7 2600 central processing unit, 3.40 GHz with 16 GB of random access memory. CPLEX returns the value of the linear relaxation of the MIP model, i.e., its optimal solution when the constraints that require integer values for decision variables are removed. CPLEX starts by generating cuts that lead to the first node of the branch and bound algorithm's tree. The output is provided when the integer optimal solution is reached or when CPLEX finds an integer solution that is within a 0.01% gap of the solution corresponding to the best node examined so far. When the minimum landscapewide wildfire resistance level was set to 0.89, CPLEX reported an optimal solution associated with a forest value of approximately 69, 127.4 × 10³. The gap to the linear relaxation was 0.015%, while the gap to the best node was 0.010%. This solution was obtained after 898.44 seconds and corresponds to landscapewide wildfire resistance levels of 0.9444, 0.8944, 0.8900, and 0.8910 in periods 1, 2, 3, and 4, respectively (Table 7).

Solutions for the MIP model, considering different levels for minimum resistance.

Table 7.

Solutions for the MIP model, considering different levels for minimum resistance.

Results of this computer run showed that, as expected, thinnings occur frequently. This management option is implemented in 270, 241, 281, and 333 stands in periods 1, 2, 3, and 4, respectively. The number of clearcuts is higher in the first period when 97 stands reach the rotation age. Fuel treatments not associated with harvest operations occur only in period 2 when 97 stands are treated. The solution to the scheduling problem reflects the initial age structure of the forest, as well as the traditional silviculture model in LNF encompassing compulsory thinnings between 20 and 50 year of age. The area of mature forest in period 1 explains the number of clearcuts in this period (Figures 1 and 2).

Figure 1.

Spatial distribution of management options when the minimum considering a minimum landscapewide wildfire resistance level was set to 0.89.

Figure 2.

Stand age in the area allocated to timber production in LNF at the beginning and at the end of the planning horizon.

The MIP model was used to assess the sensitivity of the solution to the weighting of the impacts of neighboring stands on a stand wildfire resistance (Table 8). When no impact was considered (α_is = 0), 30 stands were assigned a different prescription and the forest value decreased by 75.4 × 10³ euros. In this solution, some clearcuts previously scheduled to the fourth period are anticipated to periods 2 and 3. The number of stands where no management options are implemented in period 4 increased.

MIP solutions when the minimum landscapewide wildfire resistance level was set to 0.89 and the impacts of neighboring stands on a stand wildfire resistance are not considered α_is = 0.

Table 8.

MIP solutions when the minimum landscapewide wildfire resistance level was set to 0.89 and the impacts of neighboring stands on a stand wildfire resistance are not considered α_is = 0.

This result highlights that management options assigned to some stands to increase its specific wildfire resistance do impact too the wildfire resistance of their neighbors. Acknowledging the spatial context of stands and, thus, adjusting the specific wildfire resistance values may lead to savings in fuel treatment costs, which may even overcome the increase of present revenue resulting from the anticipation of some clearcuts.

The MIP model was also used to assess the sensitivity of solutions to the discount rate. In the case of a minimum landscapewide wildfire resistance level of 0.89, the decrease of the discount rate from 4 to 2% led to an increase of the forest value to 138,862.6 × 10³ euros (Table 9). The number of thinnings decreased and 29 prescriptions were changed. Conversely, when the discount rate was increased to 6%, the forest value decreased to 44,480.3 × 10³, while the number of thinnings increased.

MIP solutions for different discount rates when the minimum landscapewide wildfire resistance was set to 0.89.

Table 9.

MIP solutions for different discount rates when the minimum landscapewide wildfire resistance was set to 0.89.

The typical LNF silviculture model used to generate the MIP prescriptions was modified to assess whether it made sense to consider fuel treatments in the middle of periods when interventions as fuel treatment or thinning took place. For that purpose, in addition to clearcut (CCP), thinning (TC), and fuel treatment (C) options, this research considered a fuel treatment plus thinning (TCC) option and a two fuel treatments (CC) option. This change led to the increase of the number of prescriptions per stand from 40 to 233. Consequently, the number of variables increased to 20,242, of which 18,669 are binary variables (Table 10).

Report of the MIP model size, considering 233 possible prescriptions.

Table 10.

Report of the MIP model size, considering 233 possible prescriptions.

The cost of a fuel treatment planned to be implemented in the middle of a period is fully considered by the model if, and only if, the wildfire occurrence probability in the previous years in that period is zero. Else, the fuel treatment cost is multiplied by the nonoccurrence wildfire probability in the first half of the period. Computer runs, considering a minimum landscapewide wildfire resistance level of 0.89, involved the analysis of 4,242,379 nodes over 56,889 seconds. The solution reported a forest value of 69,166.2 × 10³ euros. The landscapewide wildfire resistance level reached 0.9444, 0.8948, 0.8900, and 0.8910, in periods 1, 2, 3, and 4, respectively (Table 11).

MIP solutions for different minimum landscapewide wildfire resistance target levels, considering 233 prescriptions.

Table 11.

MIP solutions for different minimum landscapewide wildfire resistance target levels, considering 233 prescriptions.

Results of this computer run showed again that, as expected, thinnings are the most frequently implemented management option in all periods. The number of clearcuts is again higher in the first period when 98 stands reach the rotation age and 26 other stands are not harvested or treated. Fuel treatments are more frequent in the second period when they are carried out in 98 stands. Very few stands are clearcut in the third and fourth periods (42 and 11, respectively). Therefore, the solution to the scheduling problem reflects again the structure of the current inventory.

The two fuel treatments (CC) option was never selected for this level of resistance, while the fuel treatment plus thinning (TCC) option was implemented in 30 stands in the third period. Apparently, these options do not impact substantially the landscapewide wildfire resistance (Tables 7 and 11), while they involve higher costs resulting from additional fuel treatments in the middle of the planning period. The selection of a limited number of these management options explains the similarity between the forest value when either 40 (Table 7) or 233 (Table 11) prescriptions per stand are considered. The fact that the MIP approach does not provide exact optimal solutions (deviation less or equal than 0.01%) may explain the differences.

Discussion and Conclusions

In this paper, we presented an approach to address wildfire risk in forest management scheduling. The approach built from wildfire occurrence and damage models (Garcia-Gonzalo et al. 2012, Marques et al. 2012) to quantify a stand specific wildfire resistance. Further, it adjusted this indicator to address landscapewide features of wildfire risk. The indicator developed in this article tries to incorporate the endogenous nature of fire risk and severity (Garcia-Gonzalo et al. 2012, Marques et al. 2012) and the relative position and elevation of neighboring stands (McArthur 1967, Van Wagner 1988). This helps take into account the neighborhood context of a stand and adjust the cost coefficients of the corresponding decision variables based on its neighbors wildfire occurrence and damage features. (e.g., Gonzáélez-Olabarria and Pukkala 2011). In fact, our approach does not try to mimic wildfire spread over landscapes, rather it demonstrates the potential of a methodology that may estimate an adjusted risk indicator at the landscape level (that is not based solely on the intrinsic risk associated to each of the stands). Such an indicator conveys information valuable to the decisionmaker. It quantifies the impact of neighboring stands on a stand wildfire resistance. This indicator may be easily updated if other wildfires occurrence and damage models are available. The models in Garcia-Gonzalo et al. (2012) and Marques et al. (2012) omit the effects of weather and surface fuels and yet the absence of reliable information about the future evolution of these factors prevents their use in long-term forest planning (Martell 2001).

An MIP model was developed to assess the potential of the adjusted wildfire resistance indicator. An innovation of this model is the attempt to address the risk of fire in forest management plans at the landscape level by including fuel treatment options in addition to the usual clearcut and thinning options. Another innovation lies in the development of a landscapewide wildfire resistance indicator and the setting of wildfire resistance targets. Previous research has focused on incorporating similar indicators in heuristics methods (Gonzáélez-Olabarria and Pukkala 2011). Still, another innovation consists of adjusting the cost coefficients of decision variables, according to the specific wildfire resistance indicators, so that the solution may take into account losses due to wildfires.

In the case study, the development of the wildfire resistance indicator was based on information from wildfire risk and damage models for maritime pine stands in Portugal. The results demonstrated that the model, although combinatorial, can be solved with reasonable computation costs. Contrary to approaches that ignore wildfire risk, it provides a more realistic estimate of forest outputs. Results further show that the approach may contribute to integrate wildfire management in the development of forest plans. It contributed to integrate wildfire management in the development of forest management plans. Specifically, it provided further information about the timing of management options such as fuel treatments to mitigate wildfires risk and damage and as such it contributes to decrease losses that might otherwise occur. It also contributed to overcome shortcomings of using solely a simple fire spread simulator where ignition points need to be simulated, as this approach may be very sensitive to the ignition points selected. The approach also provided the possibility of assessing tradeoffs between production and wildfire protection targets.

The approach highlighted that wildfire management must take into account the number and distribution of fuel treatment options. The model repeatedly chose prescriptions involving fuel treatments. These impact substantially the specific wildfire resistance and, consequently, also impact the adjusted wildfire resistance indicator. The model confirms the pertinence of the coordination of fuel treatments in a stand and in its neighboring stands, as an increased neighbor wildfire resistance may impact positively the resistance of the stand. The optimal distribution of fuel treatments helps increase the landscapewide resistance.

Hence, it may be possible to achieve an adequate landscapewide wildfire resistance by targeting key stands that may impact substantially the wildfire resistance of its neighbors, rather than implementing fuel treatments indiscriminately. This is consistent with the findings from Wei et al. (2008) and it may help save wildfire management costs. Wei et al. (2008) demonstrated that ignoring the contribution of neighbors when assessing the wildfire resistance of a stand may lead to increased costs resulting from the implementation of protection options for each stand alone.

The silviculture model used in LNF constrains the decision space, as thinnings are mandatory and must occur every 5 years between ages 20 and 25. In the future, it would be interesting to develop silviculture and growth-and-yield models to enable the expansion of the decision space. We think that the potential of this approach may then be further highlighted. Nevertheless, further research is being developed to address the decomposition of the forest-level problem into subproblems that may expand the immediate neighborhood context considered in this manuscript.

Other ongoing research to enhance this approach targets is its integration within a multiple criteria Pareto frontier method framework to facilitate the assessment of tradeoffs between wildfire resistance and production goals. It further focuses on the development of modeling features that may address adaptive management planning concerns.

Acknowledgments: This study was partially supported by the projects PEst-OE/MAT/UI0152, PTDC/AGR-CFL/64146/2006 (decision support tools for integrating fire and forest management planning) and by the scholarship SFRH/BD/37172/2007, funded by the Portuguese Foundation for Science and Technology (FCT—Fundação para a Ciência e a Tecnologia). It was also financed by the ERDF—European Regional Development Fund through the COMPETE Programme (operational program for competitiveness) within Project Flexible Design of Forest Fire Management Systems MIT/FSE/0064/2009, and by the project MOTIVE (Models for Adaptive Forest Management), funded by 7th European Framework Programme. Furthermore, this research has received funding from the European Union′s Seventh Programme for research, technological development, and demonstration under grant agreements: (i) Nr 2,82,887 INTEGRAL “Future-Oriented Integrated Management of European Forest Landscapes” and (iii) Nr PIRSES-GA-2010–2,69,257 (ForEAdapt, FP7-PEOPLE-2010-IRSES) knowledge exchange between Europe and America on forest growth models and optimization for adaptive forestry. It was also partially supported by Project PTDC/AGR-FOR/4526/2012 Models and Decision Support Systems for Addressing Risk and Uncertainty in Forest Planning (SADRI).

Appendix

Adjusted Resistance Index is between 0 and 1

It is intended to show that the adjusted resistance index varies between 0 and 1 both in the case of Equation 1 and in Equation 2. We use Banach fixed point theorem; for background see, e.g., Puterman (2005).

Theorem 1.Let 0 < w_i ≤ 1, for all i, α_is ≥ 0 for all i and s, ≤ 1 for all i, and 0 ≤ R_it ≤ 1 for all i and t. Thus, RA_it and RA_it^′ as defined by Equations 1 and 2, respectively, have their values in [0,1].

Proof

Equation 2 is equivalent having  

As 0 ≤ (1 − w_i) < 1, α_is ≥ 0 and ≤ 1, then (1 − w_i) α_is ≥ 0 for all s and 1 − (1 − w_i) ≥ 0.

Hence, RA_it^′ is a convex combination of R_it and R_st for s{i}. By assumption 0 ≤ R_it ≤ 1 and 0 ≤ R_st ≤ 1, thus 0 ≤ RA_it^′ ≤ 1.

We prove now that 0 ≤ RA_it ≤ 1. To simplify the notation we drop down the index t and rewrite Equation 1 as  

Let Q be the function defined by  

where each of its component functions is defined by  

Consider ℝ^I with the maximum norm ‖x‖ = max_iI{|x_i|}. ℝ^I is a Banach space.

Let x, yR^I. Then

 

 

 

 

 

 

 

as w_i > 0, ∀iI then 0 ≤ (1 − ) < 1. Thus, Q is a contraction. By Banach fixed point theorem, Q has a unique fixed point that corresponds to a solution of the Equation 1. Therefore, it is concluded that there is only one solution and that the sequence {r_i^k} defined by  

converges to that solution.

The aim is then to show that 0 ≤ r_i^k ≤ 1, ∀i I, ∀ k ℕ. In this case, it is proved that 0 ≤ RA_i ≤ 1.

The proof will be done by induction.

1. For k = 0, 0 ≤ r_i^k ≤ 1, ∀iI, as by assumption, r_i⁰ = R_i and 0 ≤ R_i ≤ 1, ∀iI.

2. Induction hypothesis: 0 ≤ r_i^k−1 ≤ 1, ∀iI.  

This equation is similar to Equation 21. Repeating the argument used to prove RA_it^′ [0,1], we conclude that r_i^k [0,1], ∀iI, ∀k ℕ.

Therefore, by the definition of {r_i^k}, the sequence converges to a limit RA_i and 0 ≤ RA_i ≤ 1, ∀iI.