In: Bevers, Michael; Barrett, Tara M., comps. 2005. Systems Analysis in Forest Resources: Proceedings of the 2003 Symposium. General Technical Report

PNW-GTR-656. Portland, OR: U.S. Department of Agriculture, Forest Service, Pacific Northwest Research Station.

1 M. Thor is a program leader at Skogforsk, Uppsala Science Park, S-751 83 Uppsala, Sweden; T. Möykkynen is a researcher at the Faculty of Forestry,

University of Joensuu, Finland; J.E. Pratt is a forestry officer (emeritus) at the Northern Research Station, Scotland; T. Pukkala is a professor at the Faculty

of Forestry, University of Joensuu, Finland, J. Rönnberg is a researcher at the Southern Swedish Forest Research Centre, Alnarp, Sweden; C.G. Shaw III is

a national program leader for pathology research, USDA Forest Service, Washington D.C. USA; J. Stenlid is a professor at the Swedish University of

Agricultural Sciences, Dept. Forest Mycology and Pathology, Uppsala, Sweden; G. Ståhl is a professor at Swedish University of Agricultural Sciences,

Dept. Forest Resource Management and Geomatics, Umeå, Sweden; S. Woodward is Head of Discipline Agriculture & Forestry, School of Biological

Sciences, Univerity of Aberdeen, Scotland.

Systems Analysis in Forest Resources: Proceedings of the 2003 Symposium

MODELING INFECTION AND SPREAD OF

HETEROBASIDION ANNOSUM IN CONIFEROUS 

FORESTS IN EUROPE

M. Thor, T. Möykkynen, J.E. Pratt, T. Pukkala, J. Rönnberg, 

C.G. Shaw III, J. Stenlid, G. Ståhl and S. Woodward 1

INTRODUCTION

The pathogenic fungus Heterobasidion annosum (Fr.)

Bref. causes severe problems for forestry throughout the

northern temperate zone, infecting mainly coniferous trees.

Annual losses in the European Union are estimated at €790

million. In Sweden annual revenue losses attributed to

decay and impacts on growth reach approximately €60 million. Apart from affecting income for the forest owner, H.

annosum impacts the industry: pulp and paper manufacturers and sawmills experience log quality problems because

of the fungus. 

The extent of decay is difficult to estimate in standing

trees; even so, some deduction for the risk of decay generABSTRACT

The pathogenic fungus Heterobasidion annosum causes severe problems for forestry throughout the northern temperate

zone, infecting mainly coniferous trees. Annual losses in the European Union are estimated at €790 million. The infection

biology of the fungus is well described, although variation is considerable. This paper describes two modeling projects in

Europe. The first, MOHIEF (Modeling of Heterobasidion infection in European forests), is a project within the European

Union aiming to produce a decision-support tool for forest managers, as well as a tool for scientists throughout Europe.

MOHIEF models infection and spread of the pathogen for various tree species over a range of forest conditions. The other

project, Heureka, aims to produce a system for forest planning on strategic, tactical and operative levels in Sweden. Although

national in scope, the model can be applied at regional and local levels. In Heureka the impact of H. annosum is one factor

amongst other components that include growth and yield, wood properties, biodiversity, biomass-carbon relationships and

forest-owner behavior. In this paper simulations of infection by Heterobasidion are demonstrated together with examples of

how the models can be exploited, from a forest manager’s and a scientist’s perspective.

105

ally is built into stumpage prices. The likelihood of infection can be minimized by active means such as selection 

of the right season for cutting, stump treatment, or stump

removal. These measures have a cost, however, and thus

must be integrated into forest planning to fit the “right”

operations to the “right” stands.

THE PATHOGEN

H. annosum forms perennial fruiting bodies that produce vast numbers of spores under suitable environmental

conditions. In Sweden, spore dispersal takes place during

the whole growing season (April to September) with a peak
in the summer months, whereas in warmer climates, the fall

and winter periods appear more favorable for spore production. Spores germinate on fresh woody tissue, favoring

newly cut stumps (Rishbeth, 1951a). The mycelium grows

into the stump and its roots, and infects neighboring trees

via root contacts and grafts. Damage occurs as decayed

wood on standing spruce and fir, and death of young pines.

Fungal mycelium may remain viable in old stumps for

decades, providing inoculum for the next generation of

trees. In summary, H. annosum enters a stand primarily via

stump surfaces and injuries, spreads from fresh stumps to

trees, from tree to tree and/or from old stumps to standing

trees (fig. 1). 

The incidence of H. annosum disease can be higher

under certain conditions; for example, high soil pH (Rishbeth

1951b), in first rotation stands (Rennerfelt, 1946), a high

number of stems per hectare (Venn & Solheim, 1994) and

on mineral soils with fluctuating water table (von Euler &

Johansson, 1983). On the other hand, the risks of infection

and spread appear to be lower on peat soil (Redfern, 1998),

during rainfall (for example, Sinclair, 1964) and in stands

with mixed tree species (Huse, 1983; Piri and others, 1990). 

An active way to control H. annosum is to prevent

stump infection by logging during a safe period when few

spores are present, or by treating stumps with an agent that

inhibits spore germination. In Europe protective stump

treatment, with biological or chemical agents, is carried out

on more than 200.000 hectares annually at an average cost

of €1.2 per m3 (Thor, 2002). A remedial method of control

is to remove old stumps prior to the establishment of a new

stand. This expensive action is carried out on especially

sensitive sites in the United Kingdom (UK) (Pratt, 1998).

Silvicultural measures such as planting of less susceptible

tree species, admixture of species, shortening the rotation

period, or minimizing stand entries are other options that

actively reduce spread of H. annosum. All in all, there is

ample scope and complexity to apply a modeling approach

to resolve this forest management problem. Models are

useful for forest managers as well as scientists in order to

understand the mechanisms of infection development and

impact of H. annosum.

MODELING

An empirical model of H. annosum infection in Norway

spruce exists for southern Sweden and Denmark (Vollbrecht

& Agestam 1995, Vollbrecht & Jorgensen, 1995). Mechanistic

models have been developed for H. annosum on Sitka spruce

(Picea sitchensis Bong Carr.) in the UK (Pratt and others,

1989) and for H. annosum on Norway spruce in Finland

(for example Möykkynen, Miina & Pukkala, 2000). Other

root diseases that have been modeled include Phellinus

weirii on Douglas fir (Pseudotsuga menziesii (Mirb.) Franco)

(Bloomberg, 1988). See also Pratt and others (1998) for

further details about modeling disease development in

forest stands.

The most comprehensive modeling effort to date, however, is the western root disease (WRD) model describing

the dynamics of H. annosum, Armillaria spp. and Phellinus

weirii in, for example, stands of fir (Abies spp.) and

Ponderosa pine (Pinus ponderosa Dougl. ex Laws.) in

western North America (Frankel, 1998). The WRD model

was developed by American and Canadian forest pathologists over a period of more than ten years. In addition to

the pathogens, the effect of windthrow and bark beetles can

be simulated in the model which is linked to growth and

yield models that respond to management actions. 

The WRD model has inspired a European Concerted

Action, MOHIEF (Woodward and others, 2003), in which 

a set of mechanistic models of H. annosum is to be developed for European conditions covering the range of, for

example, soil, climate, geography and hosts.

In MOHIEF (Modeling of Heterobasidion infection in

European forests) forest pathologists, forest managers and

modelers are preparing a prototype decision-support tool to

estimate potential losses from Heterobasidion attack. The

work is divided into three subtasks; 1) establish the output

106

Figure 1—Spores of H. annosum germinate on fresh woody tissue, for

example, stumps. Subsequently, fungal mycelium can spread to adjacent

root systems and stems, which are decayed. Illustration from Swedjemark

(1995).
requirements needed by the end-users of such a tool, 2)

clarify and collate the biological and management inputs

required, and 3) construct a simulation prototype model.

The project has been ongoing for 24 out of 36 months.

Draft models and an interface have been produced, and are

being tested. The range of conditions in the participating

countries – as regards climate, forestry conditions, soils

etc. – is challenging to the model. Also the basic growth

and yield models differ significantly across the participating countries. Nevertheless, significant progress has been

made: it is likely that the goals set for the project will be

achieved within the 3-year time frame.

Heureka (Lämås & Eriksson, 2002) is a research 

programme developing computer based tools for forest

analysis and planning in Sweden. Heureka uses a multidisciplinary approach that includes 13 subprojects (fig. 2).

Modeling of Heterobasidion is one of the 13 Heureka 

projects, with disease dynamics handled as in MOHIEF.

Integrated applications will be designed from the projects,

and associated systems will be developed. The project aims

at four applications to serve different user groups.

1. National and Regional Analysis 

This analyses modeling outcomes in terms of the state

of the forest ecosystem, and supply and demand for various

utilities. Basic data for this application could be the National

forest inventory, remote sensing data or information from

estate data bases. Possible users include the Board of

Forestry, the Environmental Protection Agency, forest owners’ associations and organizations for forest fuels.

2. Strategic Planning in Large Forest Enterprises

The main focus is on wood production, but other values

are also considered. Large forest areas can be handled.

Optimizing methods are used to generate a strategic plan

for the business. The application can generate future levels

of cutting, net revenues etc., and also a list of stands proposed for harvest over the next 6 months to 3 years.

3. Operational Planning

The basis for operational planning in the short term (1 to

6 months) is the list of stands available for harvest created

in the strategic planning application. The application will

optimize the net revenues from timber and forest fuels, and

allocate operations and transport.

4. Planning for Small-Scale Forestry

This application is aimed at smaller estates, and can

also handle multi-purpose forestry and environmental issues.

A problem area linked to this application is methods for

data capture that keep costs low. The users are private forest owners, forest owners’ associations, consultants and

producers of forest management plans. 

107

Figure 2—The principle of inter-disciplinary subprojects and syntheses in the multi-disciplinary Heureka project.
MODELING INITIAL CONDITIONS

Before modeling of disease dynamics, initiating data on

the incidence of root rot are needed. At present, such information in stand records is limited. Consequently, the normal situation would be to model the root rot incidence in a

stand. Current work (Thor and others 2005) uses a number

of stand and site variables to model the probability of decay

affecting a single tree. Data from the national forest inventory in Sweden were used to build the model. Records from

over 45,000 Norway spruce trees were evaluated. These

were analyzed for correlation between root rot incidence

and environmental conditions. In a stepwise logistic regression, sets of functions were developed to show significance for the variables: stand age, site class index,

temperature sum, elevation above sea level, tree diameter

at breast height, soil moisture and texture, proportion of

spruce in the stand, the occurrence of peat, and longitude.

These parameters are commonly noted in most stand

records in larger forest enterprises.

So far, the model only estimates the incidence of decay

at breast height. A calibrating function is being developed

to transform this to provide an estimate of rot frequency at

stump level. 

MODELING DISEASE DYNAMICS

For Fennoscandian conditions a model will be published

within the MOHIEF structure (Pukkala and others, 2005).

Dynamics of the disease are modeled on the basis of the

development of the stand and by the changes in fungal

dynamics over time. These sub-models are, however, not

independent of each other.

To simulate tree growth, distance-dependent or distance–independent growth models for individual trees can

be used. For cases in which there are no individual tree

measurements available, a rectangular plot is generated from

stand-level variables. Tree growth is calculated in 5 year

time steps. The timing for thinning and clear-fell operations

is triggered by basal area.

Modeling of the disease is based on the biology of H.

annosum (Woodward and others 2003, Pratt and others

1998). The probability of spore infection is estimated for

Scots pine and Norway spruce, depending on time of the

year and the use of any stump treatment. Not all stumps

that are spore infected transmit disease to adjacent trees.

Growth and infectivity of the root system depends on its

size, age and vitality. The growth rate of H. annosum is

much higher in dying roots of a stump than in a living tree. 

Once colonized, the root system is assigned a probability

for transferring the pathogen to adjacent root systems. This

process depends on soil conditions. When the fungus has

reached a neighboring tree, the growth rate in the root system is simulated, as well as the spread of decay into the

stem.

In cases of advanced decay in Norway spruce, tree

growth is reduced. Young Scots-pine trees die more frequently after attack by H. annosum than do spruce trees.

MODELING ECONOMIC OUTCOMES

The effect of disease on the economic performance of a

stand depends on growth reduction of infected trees, reallocation of produce assortments and their change in price, the

costs of prophylactic treatment, interest rate, etc.

Once the models of the initial conditions and the dynamics of the disease are in place, the economic outcome for a

number of feasible management schemes can be estimated.

One tool for modeling the outcome of assortments is

TimAn, developed at Skogforsk. TimAn is a software

package for advanced analysis of bucking-to-order, and it

covers stand inventory, construction of price lists, analysis

and follow-up of the actual outcome. TimAn uses stem

profiles for individual trees, which will be useful in these

analyses.

APPLICATIONS OF THE MODELS

One of the main objectives of the disease model is to

incorporate knowledge of root rot problems into forest

planning by means of functions built into the planning systems of various forest enterprises. Stand-alone applications

could be useful teaching tools, but will not help foresters or

forest production researchers appreciate the significance of

root rot problems.

Forester’s point-of-view

When making strategic decisions it is crucial to include

the impact of root rot in the analysis. For example, actions

by large forest enterprises will definitely affect and may

exacerbate future disease development. Issues to consider

include altering the length of the rotation, modifying thinning schedules, exercising choice of tree species and use 

of protective treatments on stumps.

In turn, these factors will affect the potential harvest

volumes and the assortment mix. One approach could be

108
scenario-handling, where different conceivable scenarios

are hypothesized and the consequences are subsequently

analyzed with the models. The models should predict average outcomes for large areas of forest: the high endemic

variation requires caution when interpreting estimates on

relatively small units, such as single stands or forests.

Strategic decisions and policies based on this type of 

analysis are a practical way to use the models.

In the short term, down to a planning horizon of one

year, the root rot situation should be considered as part of

the ordering priority in which stands will be harvested. What

are the consequences – for customers, for the supply organization, for the forest – of harvesting a Norway-spruce

stand during the low-risk period instead of the high-risk

period? What is the cost/benefit of applying stump treatment in thinning and clear-fell operations? 

Supply organizations are likely to apply the models 

differently than companies with large forest holdings. Nevertheless, the basic need to incorporate root rot modeling into

the planning process remains important for both groups.

Researcher’s Point of View

A good model will describe the disease impact and

improve understanding of the mechanisms involved in disease development. A sensitivity analysis can identify parameters critical to the outcome, and thus help to identify needs

for research. Current work suggests that some issues are

crucial, and need to be further investigated. They include

the rate of expansion of disease centers and the probability

of disease transfer in various soil types.

The type of modeling described above enables us to

have a more complete understanding of the dynamics of 

H. annosum. In the typical sequence a primary descriptive

model is developed (which has already been completed in

several countries for root disease). After a period of research,

data become available that allow for development of a

more complex mechanistic model. Ideally this model

improves our understanding so that critical processes can

be described in further detail, enabling an iterative mechanistic approach to the process of modeling to be adopted.

Used in this manner, modeling is a powerful means to

achieve a more complete understanding of the processes

involved in complex systems, such as the interactions

between pathogen, host, soil, climate and human activity.

PRESENT CONDITIONS 

AND TIME-TABLE

With regard to modeling of H. annosum, there is a welldeveloped mechanistic model for Western North America

(WRD model). In Europe a number of mechanistic models

for H. annosum dynamics in a stand are under development.

These models could be available, at least as prototypes,

before the end of 2004.

Models to describe initial conditions in Sweden will be

ready for practical use in 2004.

Models and simulation tools for detailed analyses of

assortments and cost/revenue calculations exist, but they

need to be modified to accommodate root disease dynamics.

COSTS AND REVENUES

In Sweden the primary problem with H. annosum

occurs in Norway-spruce trees. However, the variation in

conditions (rotation periods, management schemes, time of

thinning and clear-felling) suggests a clear need for reliable

quantification of the benefits of varying management

schemes and measures taken in order to minimize the

impact of root rot. Intuitively, stump treatment appears

profitable in stands growing on productive soils with short

rotations. The economic benefit is not equally obvious in

stands growing on poorer sites with longer rotation periods.

To handle this problem, tools are required; tools that the

current projects could help to develop.

In Sweden, the difference between conifer pulpwood (in

which decay is allowed to occupy 50 percent of the crosssection of a single log) and saw timber (in which no decay

allowed) is about US$ 19 (150 SEK) per cubic meter

(Brunberg, 2003). With severe decay, a fuel assortment must

be prepared, which doubles the value difference. Degradation

due to decay has a significantly higher economic impact

than a change from a higher timber class to a lower. Because

large timber volumes are involved, even a small improvement in decay reduction results in large savings. A simple

example from Sweden suffices: The annual saw-timber

volume derived from Norway spruce is about 16 million

cubic meters. Every one-percent of timber which could be

reallocated from decayed to sound wood is worth US$ 3 to

6 million (25-50 million SEK) to forest owners.

A typical stand with average Fennoscandian conditions

would produce model outputs similar to those presented in

109
Figures 3 and 4. Note that the models used are somewhat

preliminary. The output should therefore be regarded as

illustrating the magnitude of the differences, rather than

providing exact value differences.

CONCLUSIONS

Decay caused by H. annosum costs money for European

forestry, but there are means of minimizing the effects.

It is possible to include dynamics of H. annosum in

forestry planning, on strategic, tactical and operational 

levels.

Functions, models and tools for modeling H. annosum

exist in part, and further models are under construction.

Integrating functions and models of H. annosum

dynamics into existing forest growth and yield models are

likely to result in more accurate and useful results than

stand-alone applications.

LITERATURE CITED

Bloomberg, W.J. 1988. Modelling control strategies for 

laminated root rot on managed Douglas-fir stands:

model development. Phytopathology 78: 403-409.

Brunberg, T. 2003. Forestry costs & revenue in Sweden, 

2002. Skogforsk Results No. 5.

Euler, F. von and Johansson, M. 1983. Granens rotröta frodas i Norden [Incidence of root and butt rot in the

Nordic countries]. Swedish University of Agricultural

Sciences, Skogsfakta 10, 6 pp. [in Swedish].

Frankel, S. 1998. User’s guide to the western root disease 

model, Version 3.0. General Technical Report PSWGTR-165. Pacific Southwest Research Station, USDA,

Forest Service. 160 pp.

Huse, K.J. 1983. Forekomst av råte i norsk granskog. 

[Frequency of butt rot in stands of Picea abies L.

Karst]. Rapport Norsk Institutt for skogforskning 83(3),

39 pp. [In Norwegian with English summary].

Lämås, T. & Eriksson, L.O. 2002. Analysis and planning 

systems for multi-resource, sustainable forestry – the

Heureka research programme at SLU. Can. J. For. Res.

33: 500-508.

Möykkynen, T., Miina, J. and Pukkala, T. 2000.

Optimizing the management of a Picea abies stand

under risk of butt rot. Forest Pathology 30(2): 65-76.

Piri, T., Korhonen, K. and Sairanen, A. 1990.

Occurrence of Heterobasidion annosum in pure and

mixed spruce stands in southern Finland. Scand. J. 

For. Res. 5: 113-125. 

110

Figure 4—Example of results from a simulation of an average

Fennoscandian Norway-spruce stand with two thinnings (at stand age 42

and 54). Estimated decay-affected volume for two alternatives: Management

with stump treatment in previous clear fell and both thinnings (“Treatment”),

and management with logging during high-risk season and no stump treatment (“No treatment”).

Figure 3—Example of results from a simulation of an average

Fennoscandian Norway-spruce stand with two thinnings (at stand age 42

and 54). Estimated percentage of trees expressing some decay for two

alternatives: Management with stump treatment in previous clear fell and

both thinnings (“Treatment”), and management with logging during highrisk season and no stump treatment (“No treatment”).
Pratt, J.E. 1998. Economic appraisal of the benefits of 

control treatments. In: Woodward, S., Stenlid, J.,

Karjalainen, R and Hütterman, A. (eds.). Heterobasidion

annosum: Biology, Ecology, Impact and Control. CAB

International, ISBN 0-85199-275-7, pp. 315-332.

Pratt, J.E., Redfern, D.B. and Burnand, A.C. 1989. 

Modelling the spread of Heterobasidion annosum in

Sitka spruce plantations in Britain. In: Morrison, D.J.

(ed.). Proceedings of the Seventh International

Conference on Root and Butt rots, Canada, August

1988. Forestry Canada, Victoria, British Columbia, 

pp. 308-319.

Pratt, J.E., Shaw, C.G. III and Vollbrecht, G. 1998.

Modelling Disease Development in Forest Stands. 

In: Woodward, S., Stenlid, J., Karjalainen, R and

Hütterman, A. (eds.). Heterobasidion annosum: Biology,

Ecology, Impact and Control. CAB International, ISBN

0-85199-275-7, pp. 213-234.

Pukkala, T., Möykkynen, T., Thor, M., Rönnberg, J. 

and Stenlid, J. 2005. Modelling infection and spread 

of Heterobasidion annosum in Fennoscandian conifer

stands. Can. J. For. Res. 35(1): 74-84.

Redfern, D.B. 1998. The effect of soil on root infection 

and spread by Heterobasidion annosum. In: Delatour,

C., Guillaumin, J.J., Lung-Escarmant, B. and Marçais,

B. (eds.). Proceedings of the Ninth International

Conference on Root and Butt Rots, Carcans, France,

September 1997. INRA Colloques No. 89. ISBN 2-

7380-0821-6, pp. 267-273.

Rennerfelt, E. 1946. Om rotrötan (Polyporus annosus Fr.) 

i Sverige. Dess utbredning och sätt att uppträda.

Meddelande från Statens Skogsforskningsinstitut 35: 

1-88. [In Swedish with English summary].

Rishbeth, J. 1951. a. Observations on the biology of 

Fomes annosus with particular reference to East

Anglian pine plantations. II. Spore production, stump

infection, and saprophytic activity in stumps. Ann. Bot.

15: 1-21.

Rishbeth, J. 1951b. Observations on the biology of Fomes

annosus with particular reference to East Anglian pine

plantations. III. Natural and experimental infection of

Pines and some factors affecting severity of the disease.

Ann. Bot. 15: 221-246.

Swedjemark, G. 1995. Heterobasidion annosum root rot 

in Picea abies: Variability in aggressiveness and resistance. Swedish University of Agricultural Sciences,

Dept. For. Mycol. Path. Dissertation. ISBN 91-576-

5053-5.

Sinclair, W.A. 1964. Root and butt-rot of conifers caused 

by Fomes annosus, with special reference to inoculum

and control of the disease in New York. Memoir No.

391, Cornell University Agriculture Experiment Station,

New York State College of Agriculture, Ithaca, New

York, USA, 54 pp.

Thor, M. 2002. Operational stump treatment against root 

rot – European survey. Skogforsk Results No. 1.

Thor, M., Ståhl, G. and Stenlid, J. 2005. Modelling root 

rot incidence in Sweden using tree, site and stand variables. Scand. J. For. Res. 20(2): 165-176.

Venn, K. and Solheim, H. 1994. Root and butt rot in first 

generation of Norway spruce affected by spacing and

thinning. In: Johansson, M. and Stenlid, J. (eds.).

Proceedings of the Eighth International Conference on

Root and Butt Rots. Swedish University of Agricultural

Sciences, Uppsala. pp. 642-645.

Vollbrecht, G. & Agestam, E. 1995. Modelling incidence 

of root rot in Picea abies plantations in southern

Sweden. Scand. J. For. Res. 10(1): 74-81.

Vollbrecht, G. and Jorgensen, B.B. 1995. Modelling the 

incidence of butt rot in plantations of Picea abies in

Denmark. Can. J. For. Res. 25: 1887-1896. 

Woodward, S., Stenlid, J., Karjalainen, R. and 

Hüttermann, A. (eds.). 1998. Heterobasidion annosum:

Biology, ecology, impact and control. CAB International,

ISBN 0 85199 275 7.

Woodward, S., Pratt, J.E., Pukkala, T., Spanos, K.A., 

Nicolotti, G., Tomiczek, C., Stenlid, J. Marçais, B.

and Lakomy, P. 2003. MOHIEF: Modelling of

Heterobasidion annosum in European forests, a EUfunded research program. In: LaFlamme, G., Bérubé, J.

and Bussières, G. (eds.). Root and Butt Rot of Forest

Trees - Proceedings of the IUFRO Working Party

7.02.01 Quebec City, Canada, September 16-22, 2001.

Canadian Forest Service, Information Report LAU-X126. pp. 423-427.

111
This page is intentionally left blank.
FOREST

ASSESSMENT AND

PLANNING CASE

STUDIES
This page is intentionally left blank.
In: Bevers, Michael; Barrett, Tara M., comps. 2005. Systems Analysis in Forest Resources: Proceedings of the 2003 Symposium. General Technical Report

PNW-GTR-656. Portland, OR: U.S. Department of Agriculture, Forest Service, Pacific Northwest Research Station.

1 H.M. Hoganson is an Associate Professor, Dept of Forest Resources, Univ. of MN, North Central Research and Outreach Center, 1861 Highway 

169 East, Grand Rapids, Minnesota, 55744; 2 Y. Wei is a Research Associate, Dept. of Forest Resources, Univ. of MN, 115 Green Hall, 1530 North

Cleveland Avenue, St Paul, MN 55108, and 3 R.H. Hokans is Regional Analyst, USDA Forest Service Eastern Region, 626 East Wisconsin Avenue,

Milwaukee, WI 53202.

Systems Analysis in Forest Resources: Proceedings of the 2003 Symposium

INTEGRATING SPATIAL OBJECTIVES INTO FOREST

PLANS FOR MINNESOTA’S NATIONAL FORESTS

Howard M. Hoganson1, Yu Wei2, and Rickard H. Hokans3

INTRODUCTION

The Chippewa and Superior National Forests in Minnesota

released their draft forest plans in April 2003. Analyses for

the draft plans used the University of Minnesota’s Dualplan

forest management-scheduling model. Additional modeling

work is refining schedules using DPspace, a spatial model

based on dynamic programming. DPspace addresses

explicitly the core area of mature forest produced over time

in each major landscape ecosystem. Integrating the spatial

model with the Dualplan model has been a critical step in

the analysis process. 

First, this paper describes the two models and an approach

to link them. Then, the Minnesota situation is described

and test results are presented for the Chippewa National

Forest. Results presented are draft and deliberative test

runs for the National Forests. They are not linked directly

ABSTRACT

National Forests in Minnesota are currently developing new management plans. Analyses for the plans are integrating

the Dualplan forest management-scheduling model with DPspace, a dynamic programming model to schedule core area of

mature forest over time. Applications have been successful in addressing a wide-range of forest-wide constraints involving

60,000 to 100,000 analysis areas, each with potentially thousands of treatment options. A key aspect for practical application has been trimming the list of possible treatment options for each analysis area without impacting optimality characteristics of the model formulations. Draft and deliberative results for the Chippewa National Forest suggest that the core area of

mature forest can be increased substantially with little reduction in sustainable timber production levels. This is somewhat

contrary to results from an aspatial model. Results also indicate that core area concerns are more of an immediate nature

because past plans have not addressed spatial objectives. Careful planning is needed because existing core area of mature

forest cannot be replaced rapidly once it is harvested. Given more lead-time, planning can increase core area over time. 

115

with a specific forest-wide alternative presented in the draft

environmental impact statement (USDA Forest Service

2003). 

DUALPLAN

Dualplan is a forest management-scheduling model

similar to linear programming (LP) models like the USDA

Forest Service Spectrum model. Key differences with

Dualplan are that it focuses on the dual formulation of the

LP problem rather than the primal formulation (Hoganson

and Rose 1984), and it subdivides the dual problem into

many small subproblems that are each easy to analyze.

Each subproblem of the dual is a simple stand-level optimization problem for one analysis area (stand). Estimates

of the marginal costs of achieving the forest-wide constraints of the LP primal formulation are key to tying the

subproblem analyses together. These estimates are used
like market prices in the dual to value the effects of each

stand-level treatment option on each forest-wide constraint.

Initially the user estimates these marginal costs for Dualplan.

Dualplan iteratively re-estimates them after using them to

develop a forest-wide schedule. This initial schedule is

optimal mathematically, but likely infeasible because of

errors in the estimates of the marginal costs. Dualplan then

uses information about the infeasibilities of the forest-wide

constraints to re-estimate the marginal costs of the forestwide constraints. For example, if the harvest volume for 

a specific time period is constrained to meet a minimum

level and the resulting management schedule falls short 

of that level, then the estimate of the marginal cost for that

constraint is likely too low, so the estimate is increased for

the next iteration in Dualplan. In early applications of

Dualplan, like those for the Minnesota Generic Impact

Statement on Timber Harvesting in Minnesota (Jaakko

Pöyry Consulting 1994), emphasis was on satisfying forestwide constraints associated with timber production. Recently,

many options have been added to track and constrain forest

ecological conditions over time. Currently the USDA

Forest Service is using Dualplan to analyze a range of 

forest-wide alternatives for the two National Forests in

Minnesota. Besides constraints on timber production levels,

these applications have recognized constraints that help

describe desired future ecological conditions for the forest.

Constraints have addressed: (1) the age distribution of the

forest in each landscape ecosystem, (2) area targets for

selected forest cover types over time in each landscape

ecosystem, and (3) biodiversity values associated with 

having a mix of forest cover types and ages. This last set 

of constraints is similar to the concept of downward sloping

demand curves for timber where less of a forest condition

in a given time period implies a higher value for that condition in that period. 

Dualplan can use either a model I or model II type 

of formulation for addressing the stand level problems

(Johnson and Scheurman 1977). A model II formulation 

is usually more efficient because it eliminates the need to

enumerate many combinations of first-rotation treatment

types and rotation lengths with similar options for future

rotations. Dynamic programming can be used to solve the

stand level decision trees. An example decision tree (fig. 1)

helps show how the size of such trees depends more on the

number of options considered for future rotations than on

the number of options considered for the current rotation.

Options to convert to other forest cover types at the end of

the first rotation can increase substantially the size of the

decision trees for individual stands. 

DPSPACE

DPspace is a forest management-scheduling model

designed to address the spatial arrangement of the forest.

Specifically, it can track and value the amount of core

area produced. It uses a dynamic programming (DP) formulation of the problem similar to the formulation used by

Hoganson and Borges (1998) to address adjacency constraints. Also similar to Hoganson and Borges (1998), a

series of overlapping subproblems are used to address large

problems. These overlapping subproblems are similar to

moving windows used in geographic information systems

where a separate DP formulation of the problem is solved

for each window. For each window in the moving windows

process only a portion of the solution will be accepted. The

portion accepted is the portion of the window not to be

included in the next window. It is the portion farthest from

the remainder of the forest yet to be analyzed in the moving windows.

Core area is assumed to be area of the forest that meets

specific age requirements and is surrounded by a protective

116

Figure 1—A stand-level decision tree linking first rotation options with

regeneration options.
buffer. Öhman and Eriksson (1998) recognized core area as

an important spatial factor to address in forest management

scheduling models. Core area of mature forest was identified

as an important spatial measure to address in planning 

in Minnesota by the interdisciplinary planning team for

Minnesota’s National Forest (USDA Forest Service 2003).

For DPspace applications for both the Chippewa and

Superior National Forests in Minnesota focus has been on

the production of core area of mature forest. A 328-foot

(100-meter) protective buffer distance has been assumed

with minimum ages for core area varying by forest cover

type and approximately equal to the age of financial maturity for each forest cover type. Minimum age for protective

buffer also vary by forest cover type and are approximately

20 years younger than the age requirement for the core area.

DPspace uses “influence zones” to define the spatial

interrelationships that impact the production of core area.

The set of influence zones for a forest can be thought of as

simply a map of the forest that indicates, for each point in

the forest, those stands that influence that point in terms of

its potential to produce core area. Each influence zone is a

region of the forest that is influenced by the same set of

stands. Each influence zone may not be contiguous. In most

existing large stands there will be a center core area that is

influenced only by that stand itself. Other areas in the same

stand are influenced by additional stands if they are within

the core area buffer distance of the stand boundary. Generally

for a forest as a whole there are many more influence zones

than stands. Each influence zone can be labeled by the

stands that influence it, and the number of stands influencing it can be considered the influence zone’s dimension.

For example an influence zone ABC is influenced by stand

A, stand B, and stand C because all of that influence zone

is within the core area buffer distance from each of those

stands. It has a dimension of three because it is influenced

by three stands. A frequency distribution of influence zone

numbers by influence zone dimension for the Minnesota’s

National Forests (fig. 2) shows that the dimensions can 

be large and that the distribution is quite sensitive to the

assumed core area buffer distance. Many influence zones

involve three or more stands. Comparing the frequency

distributions (fig. 2) with corresponding distributions of

influence zone area (fig. 3) shows that the lower dimensioned influence zones are substantially larger in area on

average. They contain a much larger portion of the forest

117

Figure 2—Breakdown of the Minnesota National Forests in terms of the

number of influence zones by influence zone dimension.

Figure 3—Breakdown of the Minnesota National Forests in terms of the

area of influence zones by influence zone dimension. 
than the frequency distributions might suggest. With a 328-

foot core area buffer distance, relatively little area is in

influence zones involving more than six stands. As one

would expect, with larger buffer widths, more of the forest

involves more area in higher-dimensioned influence zones

(fig. 3). 

Each stage in the DP formulation used for solving each

subproblem (moving window) corresponds with a standlevel decision. The decision addressed at each stage in the

network is the choice of the management treatment option

to assign to the corresponding stand. Treatment options are

defined using a model I format (Johnson and Scheurman

1977) with a separate arc for each option. State variables

for each stage describe the treatment options selected for

some stands addressed in earlier stages of the network

(much like a simple decision tree). Each stand remains as 

a state dimension of the DP tree in later stages of the network until all other stands it “influences” have also been

represented by a stage decision. A stand influences another

stand if there exists an influence zone that both stands

influence. Each influence zone is incorporated into the

decision tree at the last stage in tree that corresponds with

one of the stands that makes up its influence zone dimension. For example, assuming influence zone ABC is influenced by stand A, stand B and stand C, and stand C is the

last of the three stands addressed in the DP network, then

the condition of the influence zone can be evaluated for all

arcs associated with stage C nodes. This condition is based

on: (1) the treatment option for the corresponding stand C

arc and (2) the corresponding conditions for stand A and

stand B as defined by the stand A and stand B state variables for the corresponding node at the start of stage C. 

Addressing the production of core area in the DP network is complicated somewhat by the assumption that the

buffer area surrounding core area need not meet the required

age of core area itself. In effect, each subcomponent of

each influence zone (subcomponents defined in terms of

the stand in which the it resides) must be examined separately. For all stands making up the influence zone, its portion of the influence zone will produce core area if the

stand meets the age requirements of core area and all other

stands in the influence zone meet the minimum buffer age

requirements for core area. When addressing each influence

zone in the DP network, these checks must be made for all

subcomponents of the influence zone. It is not simply an all

or none answer as to whether an influence zone produces

core area during each period. 

INTEGRATING COMPONENTS

A shortcoming of the DPspace model is its inability to

address forest-wide constraints directly. Its focus is on integrating spatial and aspatial considerations assuming that

the aspatial value of each stand-level treatment option is

known. Net present value (NPV) estimates for the aspatial

aspects of each treatment option for each stand is a key input

to the model. These input values need not be based strictly

on timber returns. They provide an opportunity for a direct

linkage with the Dualplan model to consider forest-wide

constraints. In the Dualplan solution process, stand-level

treatment options are valued taking into account forest-wide

constraints. Rather than view the Dualplan process as selecting one option for each stand, one can use the same approach

to estimate the aspatial value of various treatment options

to be recognized in DPspace. 

The linkage described above is straightforward. But

several factors make the problem challenging. First, the

DPspace model will be impractical to use for large problems if most stands require a large number of treatment

options to be considered in the spatial model. DPspace uses

a model I format. Converting stand-level decision trees,

like those shown in figure 1, to a model I format results in

many treatment options per stand. Second, if spatial values

addressed in DPspace are large, then the aspatial solution

may be impacted substantially causing some violations of

the forest-wide constraints in the aspatial model. In other

words, the shadow price estimates from Dualplan alone are

not likely accurate estimates of the marginal costs of the

modeled, aspatial forest-wide constraints if they don’t take

also into account the spatial aspects of the problem.

To help keep the number of treatment options small for

the spatial model without eliminating potentially optimal

treatment options, a treatment trimming model was developed that uses both the results from Dualplan and a detailed

analysis of the influence zones that define the potential for

each stand to contribute to the production of core area. Two

basic factors are considered for dropping a treatment option

in the trimming process. First, treatments are dropped if

crediting the treatment option for spatial benefits from all

stands it influences spatially cannot raise its total value above

the maximum aspatial value of the stand based on its best

aspatial treatment option. Second, a treatment option can be

dropped if there exists another treatment option that has a

greater aspatial net present value and that other option spatially dominates the treatment option to be dropped. A treatment option spatially dominates another if it meets core area

condition requirements in all time periods that the other

treatment option meets core area condition requirements.

118
To recognize the impact of spatial considerations on

achieving aspatial forest-wide constraints, Dualplan and

DPspace were linked such that the DPspace solution becomes

an intermediate solution in the Dualplan iterative process.

Essentially the subproblems of the Dualplan model are no

longer independent, single stand-level problems. Instead,

the overlapping subproblems from DPspace are the Dualplan

subproblems. Dualplan uses a summary of the DPspace

solutions in terms of the measures of the forest-wide constraints to re-estimate the shadow prices for the forest-wide

constraints. The process is repeated until schedules are

found that are “near-feasible” in terms of the forest wide

constraints. For each repetition (iteration), stand-level treatment options are re-evaluated and trimmed for use in

DPspace based on the updated forest-wide shadow price

estimates from Dualplan. 

Multiple runs of the entire system can be done to examine a range of assumptions about the value for core area.

Management schedules are easy to link with GIS systems

for further analysis because schedules produced are spatially

explicit. 

MINNESOTA APPLICATIONS

Dualplan has been used to analyze seven forest-wide

alternatives for each National Forest in Minnesota. These

alternatives have focused on sustaining timber harvest levels

over time while moving the forest closer to desired future

conditions. Desired future conditions are defined specifically

for each of seven landscape ecosystems in the Chippewa

National Forest and each of eight landscape ecosystems in

the Superior National Forest. Forest-wide alternatives differ

in terms of assumptions defining desired future conditions

and the rate at which the forest is moved towards the

assumed desired future conditions. 

For both forests the planning horizon consists of ten

ten-year periods. Analysis areas (AA’s) are spatially explicit

with over 67,000 AA’s for the 549,000 acres modeled for

the Chippewa National Forest. Over 101,000 AA’s are used

for the 1,212,000 acres modeled for the Superior National

Forest. Many of the AA’s are substands representing riparian

areas along lakes and streams. Seventeen silvicultural treatment types were recognized, representing various types of

even-aged and uneven-aged management. For each treatment type, a range of harvest timings was considered. Twenty

forest cover types were recognized for each forest with five

timber-productivity site-quality classes for most all types.

Desired future conditions for each landscape ecosystem were

defined in terms of a desired mix of forest cover types and

a desired stand age distribution. Natural succession was

modeled recognizing that some stands change forest cover

type over time without treatment. Estimates of pre-settlement age distributions and stand replacement fire intervals

served as a benchmark for defining desired future conditions and the rate at which harvesting is applied to mimic

stand replacement fires.

Forest cover type conversion (restoration) options were

an important consideration for stand-level decisions because

desired future forest-wide conditions were defined in terms

of a desired forest cover type mix. This added substantially

to the number of treatment options modeled for most stands.

Five key map layers influenced the set of treatment options

considered for each stand. Besides the landscape ecosystem

and riparian map layers described earlier, other layers

included a map of management areas, a map of scenic classes

and a map identifying sensitive areas for a number of plant

and animal species. The map of management areas included

18 classifications used by the Forest Service to help define

the theme of each alternative. For some forest-wide alternatives that emphasize older forest conditions, management

area classifications limited treatment options substantially. 

Tests of the spatial model for the National Forest applications focused on a forest–wide alternative somewhat similar to the alternative with the highest timber harvest level

in the draft forest plan. That alternative was identified in

the draft environmental impact statement (USDA Forest

Service 2003) as the forest-wide alternative with the most

fragmentation of mature forest over time. That alternative

would also likely be the most complicated alternative to

model simply because, on average, it is the least restrictive

and thus involves more possible treatment options for each

analysis area. Multiple model runs of the integrated modeling system were applied, varying the value of mature forest

core area between applications. For each application, multiple iterations of the integrated model were needed to account

for the impact of valuing core area on the marginal cost

estimates for the many aspatial constraints. These aspatial

constraints defined desired future conditions for each landscape ecosystem and sustained timber harvest levels over

time. For all the values of mature forest core area examined,

the aspatial constraints were satisfied after multiple iterations of the integrated modeling system. The resulting management schedules are potentially improved schedules in that

they are near-feasible in terms of the aspatial constraints

and produce more core area of mature forest.

Core area is valued in terms of dollars per acre per

decade with values assumed to occur at the end of each

decade. Assuming a four percent annual discount rate, a

price of $100 per acre per decade for core area, and that the

midpoint of the first decade is used as the base year for

119
NPV calculations, then an acre that produces core area in

every decade over an infinite horizon has a NPV of

approximately $253 per acre.

Results of the test applications were similar for both

National Forests. Here, results are presented only for

Chippewa National Forest. Two types of core area were

valued and tracked for both forests: upland forest core area

and lowland forest core area. Both types had core area

minimum stand age requirements that varied by forest

cover types with a minimum age of at least 40 years for all

forest cover types. Core area value assumptions impacted

core area production levels quite similarly for both uplands

(fig. 4) and lowlands (fig. 5). Results for the $0 per acre

per decade value are equivalent to the Dualplan results

where core area was not valued. Without valuing core area,

core area levels decline substantially over the first two

decades (fig. 4 and fig. 5). It should not be surprising that

the draft environmental impact statement for the Minnesota

Forest Plan (USDA Forest Service 2003) raises concern

about potential declines in core area of mature forest. Raising

the core area values to $100 per acre per decade raises the

amount of mature area for both uplands (fig. 4) and lowlands (fig. 5) in the later decades, but there is still a decline

in the early decades from the starting condition. A core area

price of approximately $300 per acre per decade is needed

before such a drop is not present in the early decades. For

the core area values modeled, higher values for core area

led to higher core area output levels for both uplands and

lowlands, with levels substantially higher in the long-term

than in the short-term. Increases are less in the short-term

because of the time required to produce larger blocks of

mature forest that are relatively well suited for producing

core area. In past planning efforts producing core area of

mature forest was not a primary objective, so it is not surprising that the intermediate-aged stands are not arranged

spatially such that they can produce large quantities of core

area when they soon reach the minimum age requirements

for core area of mature forest. With more lead-time and

planning, substantially more core area can be produced in

the longer-term (fig. 4 and fig. 5). 

120

Figure 4—Impact of core area values on upland core area production for

five model runs for the Chippewa National Forest. 

Figure 5—Impact of core area values on lowland core area production for

five model runs for the Chippewa National Forest.
Recognizing higher core area values tended to shift

optimal management strategies from even-aged management

to uneven-aged management (fig. 6). However, even-aged

management was still predominant even at the highest core

area values examined. Higher core area values also tended

to lengthen rotation ages. As one might expect, higher values for core area resulted in larger patches of mature forest.

Results did not suggest that current harvest blocks should

be large so as to later produce more core area. With levels

of core area more difficult to increase in the short- term,

designing harvests to focus more on what remains on the

landscape in the short- term seems of more importance

than the size and shape of harvests. Focus is more on what

remains on the landscape rather than on what is harvested.

Model applications were performed using a Dell workstation with two 1.8 GHz processors and 512 MB of memory. Each iteration of DPspace took approximately one

hour for the Chippewa National Forest and two hours for

the Superior National Forest. Several heuristics have been

used to help reduce the number of iterations needed for

DPspace model. In effect, more iterations of the much faster

Dualplan model are applied with heuristics incorporated to

tie solutions to the DPspace results. These heuristics are

not applied in the last iteration of the process so as not to

impact optimality characteristics of the final solution.

Solutions that are close to satisfying all of the aspatial constraints have generally been found in approximately 10 to

20 iterations of DPspace. Publications to describe and test

the specific details of the modeling process are in progress.

The overall process for integrating the models shows enormous potential for improvement as intermediate iterations

of the DPspace model unlikely need to use all stands or all

influence zones. Early and intermediate iterations of the

process can also likely use smaller subproblems (moving

windows). The ability to use current shadow price estimates to screen and trim stand level alternatives is a key

component of the modeling system. Like with Dualplan,

intermediate iterations of the integrated model are done

simply to help find better estimates of the shadow prices

associated with the forest-wide aspatial constraints. Once

good estimates are found, the DPspace model seems quite

effective at finding near-optimal spatial solutions for the

problem as it is formulated mathematically

Future work is needed to better automate system details

for use by others. Improvements in the system can likely

be achieved by better understanding how heuristics might

be best used to help speed the solution process. Clearly, a

key in application is having the ability to apply many model

runs to learn about key model assumptions and trade-offs.

In future research it seems desirable to consider options for

recognizing multiple types of core area and other spatial

measures. Linkages to forest plan implementation are also

important, as issues are complex without simple guides for

managers. The ability to decompose large problems into

parts makes the system especially appealing. It suggests

opportunities for large-scale collaborative planning involving multiple ownerships and model expansion opportunities

to consider additional facets of the problem.

CONCLUSIONS

A modeling system closely tied to optimization techniques has been developed to better integrate timber harvest

scheduling with ecological objectives of moving forests

closer to desired future conditions. Desired future conditions are complex conditions that are difficult to define,

varying by ecological land classes and involving species

and age composition of the forest as well as its spatial

arrangement. Core area of mature forest has been the spatial

measure of focus in modeling applications by the National

Forests in Minnesota. Test applications have found that relatively high timber harvest levels can be sustained while

increasing the amount of core area of mature forest over

121

Figure 6—Impact of core area values on the type of stand-level management options selected for five model runs for the Chippewa National Forest.
time. Spatial arrangement of the forest has not received

attention in past planning efforts so it is not surprising that

considerable time is needed to improve spatial conditions

over time. Without recognizing spatial detail in planning,

timber harvesting can reduce the existing core area of

mature forest quite rapidly. 

LITERATURE CITED

Hoganson, H.M.; Borges, J.G. 1998. Using dynamic 

programming and overlapping subproblems to address

adjacency in large harvest scheduling problems. Forest

Science. 44(4): 526-538.

Hoganson, H.M.; Rose, D.W. 1984. A simulation 

approach for optimal timber management scheduling.

Forest Science. 30(1): 220-238.

Jaakko Pöyry Consulting, Inc. 1994. Generic 

Environmental Impact Statement on Timber Harvesting

and Forest Management in Minnesota. Tarrytown, NY:

Jaakko Pöyry Consulting, Inc. 813 p.

Johnson, K.N.; Scheurman, H.L. 1977. Techniques for 

prescribing optimal timber harvest and investment

under different objectives—discussion and synthesis.

Forest Science Monograph No. 18. 31 p.

Öhman, K.; Eriksson, L.O. 1998. The core area concept 

in forming contiguous areas for long-term forest planning. Canadian Journal of Forest Research. 28: 1032-

1039.

USDA Forest Service. 2003. Draft environmental impact 

statement: forest plan revision: Chippewa National

Forest and Superior National Forest. Milwaukee

Wisconsin. 695 p.

122
In: Bevers, Michael; Barrett, Tara M., comps. 2005. Systems Analysis in Forest Resources: Proceedings of the 2003 Symposium. General Technical Report

PNW-GTR-656. Portland, OR: U.S. Department of Agriculture, Forest Service, Pacific Northwest Research Station.

1 Brian McGinley and Gary Marsh are Resource Planners, 2 Allison Reger is a Forest Resource Analyst, and 3 Kirk Lunstrum is a Wildlife Biologist. 

Coauthors work on the Willamette National Forest, Oregon.

Systems Analysis in Forest Resources: Proceedings of the 2003 Symposium

LANDSCAPE CHANGES FROM 

MANAGING YOUNG STANDS 

A Fall Creek LSR Modeling Study

Brian McGinley1, Allison Reger2, Gary Marsh1, and Kirk Lunstrum3

INTRODUCTION

Public land agencies continually seek a better understanding of spatial and temporal changes to forested 

landscapes and to use this understanding in successive

management decisions. Recent focus has moved toward

complex assessments of forest habitat changes across

large landscapes. This shift has been partially prompted 

in the Pacific Northwest by desires to protect old growth

forests and its dependent species. 

The 1994 Northwest Forest Plan, with a network of

Late-Successional Reserves (LSRs), is an example of this

focus on managing forest habitat changes across large-scale

landscapes. One option in the LSR network strategy is

123

using silviculture treatments to improve late-successional

habitat conditions within young stands less than 80 years

old. 

Thinning is accepted as a useful tool for accelerating

vegetative development toward late seral conditions

(USDA,USDI 1994; USDA,USDI 1998). Although much

debate still remains on the best treatments or stand-level

designs to apply, thinnings can effectively influence species

composition, stem growth, canopy complexity, and understory communities in residual stands. 

With most thinning in LSRs occurring in established

plantations, managers should anticipate the temporal and

spatial influence these plantations can have in shaping the

ABSTRACT

Public land agencies continually seek a better understanding of spatial and temporal changes in forested landscapes to guide

successive management decisions. The 1994 Northwest Forest Plan developed a network of Late-Successional Reserves (LSR)

with the primary goal of protecting and improving late-succesional habitat. Improving LSR habitat can involve silviculture

treatments in stands less than 80 years old, which would improve stem growth, canopy height complexity and understory community development. The Fall Creek Modeling Study, using TELSA (Tool for Exploratory Landscape Scenario Analysis),

explores connections among young stand management strategies and resulting landscape habitat patterns. TELSA is a spatially

explicit planning tool that simulates succession, natural disturbances, and management actions on a landscape over time. The

Fall Creek LSR is a typical mosaic of young stands, fragmented late successional habitat, and high road densities created by

intense forest management. The impetus for this project was insufficient funding to promote late-successional habitat by the

thinning of all available hectares over the next 50 years. This modeling study evaluated a number of resource-driven scenarios

for selecting prescriptions and prioritizing stands for treatment. The team used spatial mapping and stand simulation tools to

assess the effects of these thinning scenarios on LSR habitat patterns over a 200 year period. For specific landscape objectives,

such as interior habitat or mimicking natural fire patterns, location of habitat improvements from thinning are shown to be

more valuable than total hectares of improved habitat.
landscape patterns of late-successional habitat. Plantations

comprise around 25 to 35 percent of most large LSRs in

the Mid-Willamette network, making their landscape role

significant. 

The Fall Creek Modeling Study explores the connections among young stand management strategies and

resulting landscape habitat patterns.

STUDY AREA

The Fall Creek watershed and LSR were combined to

create a 36,338 hectares analysis area on the Willamette

National Forest (fig. 1). While 20 percent of the area is

classified in the Forest Plan as outside the LSR boundary,

this analysis was simplified to focus only on LSR objectives.

The Fall Creek analysis area is well fragmented with

about 45 percent (16,498 hectares) occurring in plantations. Only one-third of plantation hectares have been

thinned prior to age 20, with another one-third beyond 

the optimal age range for precommercial thinning.

Over the next 50 years, roughly 22,719 hectares will

potentially be available for thinning, assuming multiple

entries are prescribed in many candidate stands. This

includes 2,328 hectares of 40 to 80 year old stands currently

available for thinning in the next 10 years.

The Fall Creek LSR is well roaded with some of the

highest road densities in the Mid-Willamette LSR network. About 89 percent of Fall Creek hectares have a

road density of greater than 1.2 km/sq. km. A minimum

124

Figure 1—Existing Stand Conditions within the Fall Creek LSR.
road network for this LSR includes only 23 percent (226

kilometers) of its total road system.

MANAGEMENT SITUATION

The Forest Plan identifies two principal goals for thinning stands less than 80 years old in LSR’s: 1) to develop

old growth forest characteristics, and 2) to reduce the risk

of large-scale disturbance events like wildfire (USDA,

USDI 1994, pp. B-5-6).

The Mid-Willamette LSR Assessment identifies more

interior habitat and lower road densities as key management objectives for the Fall Creek LSR (USDA, USDI

1998). Enhancing late seral connective corridors along its

eastern boundary is another focus for the Fall Creek LSR.

While plenty of thinning opportunities exist in the analysis

area, the Willamette National Forest does not expect to

have sufficient funding or staffing to thin all available

hectares over the next 50 years. 

Current thinning efforts are using the same oldest/largest/

closest-first priority strategy used in placing clearcuts on

the landscape. A thinning strategy built from resource

objectives for prioritizing stands is likely to create greater

landscape habitat benefits.

MANAGEMENT SCENARIOS

This analysis evaluated a number of resource-driven

management scenarios for selecting thinning prescriptions

and prioritizing young stands for treatment. The following

briefly describes the scenarios examined. 

Management Scenarios

Scenarios are comprised of management constraints, a

stand priority strategy, and thinning prescriptions. Commercial

thinning prescriptions could potentially occur at one or

more of three stand ages (40, 60 and 80 years) and at one

of three treatment intensities (light, medium, heavy). There

were six scenarios studied, each meeting a single resource

objective or providing a benchmark for comparison.

1) The No Action scenario was modeled to provide a

point of comparison. In this scenario no action was taken

on any stands permitting the managed stands to now follow

natural succession.

2) The Owl Sites scenario focuses on young stands

within 2.0 kilometers of nesting activity centers and uses

nesting success to prioritize activity centers. The objective

of this scenario was to improve habitat near northern spotted owl nest sites.

3) The Roads scenario focuses thinning on habitat

improvement and reducing road densities in land blocks

(roadsheds) possessing high aquatic risk ratings (USDA,

1999). Thinning entries into roadsheds were spaced at

least 20 years apart to allow for road closures and decommissioning. The objective of this scenario was to reduce

aquatic risks through habitat improvement and lower road

densities.

4) The Natural Fire scenario uses thinning to create

habitat conditions predicted from defined natural fire

regimes. Fire frequency, fire intensity and slope position

were used to define resulting habitat conditions and

assign treatments. This scenario was developed to mimic

vegetative patterns created by natural fire regimes.

5) The Interior Habitat scenario focuses on young

stands surrounded by or next to large interior habitat blocks.

Thinning intensity was reduced to minimize edge effects

on existing late-successional habitat and to expand large

interior habitat blocks.

6) Finally, a Two-Thin Benchmark was applied to 

all young stands to provide an upper limit of thinning for

comparing with other scenarios. This scenario was the

most aggressive thinning strategy.

Several iterations of each scenario were tested to look at 

a range of possibilities.

Analysis Tools

The analysis team used spatial mapping and stand simulation modeling tools to assess the effects of thinning treatments on LSR habitat patterns. TELSA (Tool for Exploratory

Landscape Scenario Analysis) is a strategic, landscapelevel planning tool that can project planning strategies onto

a landscape and assess long-term spatial and temporal

effects using a range of performance indicators (ESSA,

1999).

PNW_GAP (originally ZELIG-PNW) is a stand growth

simulator used to track stand attributes along development

curves in response to thinning treatments (Garman, 1999).

Garman (1999) tested a variety of thinning prescriptions

for accelerating the development of late-successional stand

attributes. These tools were used to apply thinning schedules for each scenario to the landscape and project late-successional attribute development over a 200-year time period.

Late-Successional Stand Attributes 

with Minimum Thresholds

Garman (1999) simulated a wide range of thinning

125
treatments on a 40-year old Douglas-fir stand that had been

precommercially thinned at age 12 years. He tracked thinning effects on five late-successional attributes and the

time needed to reach threshold values. (See attributes with

minimum thresholds below). 

Garman also used a composite Late-Successional

Index (LSI) for tracking when stands meet all three live

stand attribute thresholds. The Fall Creek analysis team

focused only on the three live stand attributes and LSI to

evaluate scenario influences on late-successional habitat

and other management goals. The snag and downed wood

attributes were not used in this study, because they could

be readily created (topping or felling) once minimum

large bole densities were achieved. 

Thinning prescriptions selected from Garman’s simulations for this analysis quickly met thresholds for the live

stand attributes and matched other scenario objectives.

Late-successional Criterion

The TELSA spatial model was used to map stand attributes across the landscape at different time points for each

scenario. Other key landscape attributes were used to evaluate scenarios (see below).

BUILDING RESOURCE SCENARIOS

TELSA, a spatially explicit model, requires data in a

vector-based format. Several GIS layers used in this analysis included vegetation and planning area boundaries. 

Using GIS technology the landscape was stratified into

planning zones. Planning zones signify areas with discrete

management goals or areas requiring separate reporting of

results. Planning zones common to all scenarios were the

project area boundary and transportation routes.

Transportation routes entered into TELSA tracked which

roads were used during thinning operations and how long

each road was left inactive (unused by thinning activities).

A combination of unique planning zones and thinning

prescriptions customized each scenario to achieve key

resource objectives highlighted by the scenario. In response

to limited budget assumptions, the Fall Creek analysis team

also created priority management goals for thinning within

planning zones. These unique planning zones defined the

key difference between resource scenarios. Each scenario

had 1 to 2 additional planning zones added for stratifying

the landscape based on differing management goals.

Planning zones specific to each scenario are discussed

below.

1) No Action Scenario: This scenario used no unique

planning zones and no management was applied to the

landscape. The scenario provided a point from which to

measure the effects of thinning stands.

2) Owl Scenario: This scenario focuses habitat improvement in young stands near known spotted owl (Strix occidentalis caurina) activity centers. Based on monitoring data

on the nest site occupancy and reproductive success, owl

habitat centers were ranked into three site categories. A

1.1-km buffer around nest sites and the sites corresponding

ranking became planning zones in the TELSA model. A

two-thin prescription (a moderate thin at age 60 years, light

thin at 80 years) was selected for all stands.

This prescription was selected by the analysis team 

to maintain forest cover and canopy closure for potential

owl dispersal and foraging opportunities, while attaining

threshold values of late-successional attributes quickly.

Owl scenario iterations differed in which owl circles

were treated and which were left to natural succession.

3) Roads Scenario: This scenario divides the landscape

into management blocks called “roadsheds”—lands that

would feasibly log to a local road system. Roadsheds varied in size from 1,500 to 4,000 hectares. Two criteria were

used to prioritize roadsheds for thinning and resulting road

closures:

a. Aquatic risk rating

b. Plantation hectares available for thinning

The roadsheds and priority ranking were integrated into

the TELSA model as a planning zone.

126

1. Density of large boles >100 cm

• 10 boles per hectare

2. Canopy Height Diversity Index

• Index = 8.0

3. Density of shade tolerant boles >40 cm

• 10 boles per hectare

4. Density of snags >50 cm

• 10 boles per hectare

5. Mass of logs >10 cm at large end

• 30 Mg/ha
4) Natural Fire Scenario: Fire regimes defined in the

Integrated Natural Fuels Management Strategy (USDA,

USDI 2000) served as a template for the Natural Fire scenario. These fire regimes as identified for the study area

became a key planning zone map in the model. These three

regimes can be briefly described as:

1. moderate frequency/mixed severity; 

2. low frequency/mixed severity; and 

3. low frequency/high severity. 

Slope position (lower 50 percent of the slope and the

upper 50 percent of the slope) was a second planning zone

added to include the variation of fire effects due to slope

position. Thinning prescriptions were selected to match the

stand conditions likely to result from representative fire

events under these regimes. Moderate frequency areas

received a two thin prescription (upper slopes – heavy at

age 40 years, moderate at 80 years; and lower slopes –

moderate at 40 years, light at 80 years). Low frequency

areas received one heavy thin at age 60 years in the upper

slopes only.

5) Interior Habitat Scenario: Old growth dependent

species thrive best in late succesional interior habitat,

which provides higher relative humidity, moderated temperatures, and desirable physical structure. Interior habitat

(IH) is also generally free of invasive and competing nonnative species. As a rule, interior habitat conditions in latesuccessional stands improve as you move further from

edges with early seral conditions.

This scenario was designed to focus thinning in young

stands that help expand existing interior habitat (IH) blocks.

Four landscape blocks were identified for prioritizing treatments and are listed by priority:

1. Stands within largest IH blocks

2. Stands creating sinuosity on the edge of the largest IH

block

3. Stands connecting large IH blocks

4. All other stands

Because thinning may increase the risk of edge effects

(edaphic changes and blowdown) on adjacent late-successional habitat, prescriptions were varied based on distance

from late-successional edge: 

• Within 120 meters of habitat edges—Moderate thin at 

60 years, and light thin at 80 years.

• Beyond 120 meters from habitat edges—Heavy thin at 

60 years, and moderate thin at 80 years.

Landscape blocks and distance from edge both became 

key planning zones in the TELSA model.

6) The Two Thin scenario prescribed a two-thin prescription (heavy thin at age 40 years, medium thin at 60

years) to all candidate stands with a constraint that reentry

into any roadshed should not occur for at least 20 years.

This constraint reduced open road densities by increasing

opportunities for road decommissioning work.

Comparing Resource Scenarios 

The following management goals were selected to 

compare resource scenarios:

• Total hectares of late-successional habitat (an LSI rating 

above 75 percent)

• Total hectares of interior late-successional habitat

• Total hectares of late-successional habitat within owl 

nesting circles.

• Total hectares of late-successional habitat within 

moderate fire intensity areas.

• Total kilometers of roads remaining inactive for 10 and 

20 years.

Management goals were also expressed as relative measures of total hectares thinned. Resource scenarios were also

evaluated and compared visually, as seen in figures 2 and 3,

to generally understand how landscape habitat patterns

change under each scenario. Finally the flow of thinning

hectares across the analysis timeline was considered for

operational reasonability.

ANALYSIS RESULTS

After 70 years, resource scenarios had thinned between

8,600 and 13,000 hectares. By contrast, the Two-Thin

Benchmark thinned 31,000 hectares in just 50 years and the

No Action scenario 0 hectares. Harvest levels by decade

varied noticeably between scenarios. The Roads scenario

created a tri-modal harvest pattern due to the 20-year harvest constraint for roadsheds, while other scenarios more

equitably spaced harvest over time.

By analysis year 80, scenarios had increased late-successional habitat (hectares with an LSI rating above 75 

percent) by 2 to 17 percent above the No Action scenario.

Scenarios had increased LS habitat within owl circles by 

3 to 13 percent above the No Action scenario.

127
128

Figure 2—Description of Habitat using LSI values by Scenario at Year 80.
129

Figure 3—Description of Habitat using LSI values by Scenario at Year 120.
By analysis year 80, scenarios had increased interior

habitat by 2 to 40 percent above the No Action scenario.

The Roads and Interior Habitat scenarios achieved highest

gains in interior habitat.

By analysis year 40, scenarios were able to keep 327 

to 395 km of roads inactive (or unused) for more than 20

years. These roads are considered good candidates for closure. The Roads scenario created the most kilometers of

inactive roads, largely due to the 20-year re-entry scheduling constraint and its narrow focus on only high-risk roadsheds.

Habitat patterns varied notably between scenarios,

though total late-successional habitat differences were not

great (figures 2 and 3). The Owl Site scenario concentrates

habitat improvement efforts around only 37 nest sites spread

across the landscape. The Natural Fire scenario also tended

to spread thinning treatments across the landscape, in high

fire frequency areas.

By contrast, the Roads scenario focused thinnings on

the north side of the landscape, where aquatic risk ratings

were highest. Similarly, the Interior Habitat scenario concentrated treatments within the core of the landscape along

Fall Creek.

Efforts to include precommercial thinning options in 

the analysis were stymied by incompatible data sets of

stand development curves, and were eventually abandoned.

CONCLUSIONS

Modeling helps to visualize the cumulative benefits/

impacts of thinning strategies on LSR landscapes and 

compare management strategies. The Fall Creek analysis

should provide useful reference points for future project

analysis or creation of a 5-year action plan of thinning

projects.

Results should prove useful for building public and

agency confidence in using thinning to improve habitat

conditions within LSR landscapes. For interior habitat values, location of habitat improvements on the landscape can

be more important than total hectares of habitat improved.

Landscape objectives help drive this placement of thinning

efforts. Landscape attributes or habitat indices used to compare management scenarios are more useful when made

relative to total treatment hectares.

NEXT STEPS IN THIS INVESTIGATION

• Model other resource scenarios of interest to managers.

¡ Connectivity (travel corridors) between LSRs.

¡ Combine several resource objectives into one scenario.

¡ Allocate resource scenarios to different sections of the 

LSR and pursue simultaneously.

• Strengthen the simulations by testing complicated 

prescriptions involving precommercial thinning.

• Introduce implementation constraints (planning, budget, 

resource issues) into the scenarios.

• Use modeling results to create a five-year action plan 

of thinnings.

LITERATURE CITED

Garman, S.L. 1999. Accelerating development of latesuccessional conditions in young managed Douglas-fir

stands: a simulation study. USDA, Forest Service,

Willamette National Forest. On file with: Willamette

National Forest, 211 E. 7th Avenue, Eugene, OR 97401.

ESSA Technologies LTD. 1999. TELSA: Tool for 

Exploratory Landscape Scenario Analyses, User’s Guide

Version 1.0. Vancouver, BC: ESSA Technologies Ltd. 

164 p.

USDA Forest Service and USDI Bureau of Land 

Management. 1994. Record of Decision for Amendments

to Forest Service and Bureau of Land Management

Planning Documents within the Range of the Northern

Spotted Owl. http:/www.oregon.gov/library/

USDA Forest Service and USDI Bureau of Land 

Management. 1998. Mid-Willamette Late-Successional

Reserve Assessment. On file with: Willamette National

Forest, 211 E. 7th Avenue, Eugene, OR 97401.

USDA Forest Service.1999. Fall Creek/Little Fall Creek 

Access and Travel Management Environmental Assessment. Willamette National Forest. On file with:

Willamette National Forest, 211 E. 7th Avenue, Eugene,

OR 97401.

130
In: Bevers, Michael; Barrett, Tara M., comps. 2005. Systems Analysis in Forest Resources: Proceedings of the 2003 Symposium. General Technical Report

PNW-GTR-656. Portland, OR: U.S. Department of Agriculture, Forest Service, Pacific Northwest Research Station.

1 P. Bettinger is an Associate Professor, Warnell School of Forest Resources, University of Georgia, Athens, GA 30602. 2 J. Drake is an Assistant Professor,

Warnell School of Forest Resources, University of Georgia, Athens, GA 30602.

Systems Analysis in Forest Resources: Proceedings of the 2003 Symposium

EXAMINATION OF TRADEOFFS BETWEEN 

TIMBER PRODUCTION AND SPOTTED OWL

HABITAT IN COASTAL OREGON

Pete Bettinger1, and Jason Drake2

INTRODUCTION

As forest and wildlife sciences evolve, researchers and

practitioners become increasingly exposed to new hypotheses concerning the link between management (and hence

landscape condition) and ecological and social values. The

development of new hypotheses also provides opportunities

for their integration with management models and subsequent evaluations of impact, effect, and association. Advances

in wildlife habitat relationships research have provided forest planning the opportunity to more closely incorporate

habitat goals into forest planning processes. The type of

mathematical programming approach used to accomplish

this integration will depend on the wildlife species under

consideration and type of relationships used to evaluate the

capability of a landscape to provide habitat. One example

of new research on wildlife-habitat relationships is that

provided by McComb and others (2002), who developed 

a wildlife capability model for the northern spotted owl

(Strix occidentalis caurina).

ABSTRACT

A forest planning model was developed to assess the trade-offs among timber production and spotted owl (Strix occidentalis caurina) habitat in coastal Oregon. The model utilized a threshold accepting heuristic technique to schedule, via 1-opt

moves, timber harvests subject to spotted owl, adjacency, and green-up goals. The forest planning model was designed as a

Model I problem, and only clearcut harvest decisions were made. The trade-offs among timber production and habitat goals

are presented as a set of production possibility frontiers. With this preliminary analysis, it seems that there is a linear relationship between even-flow timber harvest levels and an allowable decrease in habitat levels (from projections of habitat

under a “grow-only” scenario), up to a point. When habitat levels were allowed to decrease beyond 20% of the “grow-only”

levels, no effect on timber harvest levels was noted. While this model helps managers and policy makers understand the

trade-offs among timber production and spotted owl habitat levels, results may be limited to the landscape examined and

the range of constraints tested.

131

The McComb and others (2002) owl habitat capability

index (HCI) model includes both a within-stand nesting

component and a home-range analysis, and thus requires

measuring the characteristics of forest vegetation up to

2,400 m from each area of interest. Since a management

activity scheduled to a management unit can affect the

habitat quality of other management units up to 2,400 m

away, this suggests the need for integer decision variables.

A number of mixed integer and integer formulations of forest planning problems have been addressed in the forestry

literature, from the pioneering work of Hof and Joyce (1992),

Hof and Raphael (1993), and Hof and others (1994) using

the exact techniques, to approaches provided by Bettinger

and others (1997, 2002) and Calkin and others (2002)

using heuristic techniques. This paper presents a process

that integrates McComb and others (2002) HCI model with

commodity production goals in a heuristic forest planning

model.
METHODS

We selected a large landscape in western Oregon to use

in modeling the integration of spotted owl HCI and commodity production goals in a forest planning context. The

landscape (fig. 1) represents approximately 235,000 ha of

land, mostly forested, and mainly managed by state and

industrial land managers. Pixels in figure 1 are colored

based on their vegetation class, and while not completely

informative of the vegetative structure across the landscape,

the figure allows one to locate the agricultural, urban, and

cleared areas (the light colors), and to discern changes in

topography based on the variations in the gray scale. The

vegetation database was acquired from the CLAMS project

(CLAMS 2002), and consists of a raster GIS database containing 3,752,714 pixels with a spatial resolution of 25 m.

Each pixel in the vegetation GIS database has associated

with it a tree list, as the database was created using a gradient nearest neighbor approach to classification that facilitates this association (see Ohmann and Gregory 2002). A

tree list is simply a list of trees used to represent the forested

conditions within each pixel. The lists used were derived

from the U.S. Forest Service’s Forest Inventory and

Analysis program, the U.S. Forest Service’s Continuous

Vegetation Survey program, and research plots. Ohmann

and Gregory (2002) provide more detailed information

regarding the classification approach and assignment of

tree lists to pixels in the GIS database.

The vegetation data associated with each pixel includes

standing timber volume, average tree age, quadratic mean

diameter, and several measures of coarse woody debris.

This data facilitates measuring the impact of activities on

both commodity production and habitat goals. The age

class distribution of this area (fig. 2) shows that most of the

forest vegetation is between 20 and 50 years old. A GIS

database describing management units was acquired from

the CLAMS project. Management units average 6-8 ha in

size. The spatial relationship of each management unit in

relation to others of which it shares a side was also developed, to facilitate the modeling of adjacency and green-up

restrictions. Pixel data within management units are aggregated to arrive at average stand ages, volumes, and so on,

for each management unit.

Objective function

The forest planning problem we solved attempts to

maximize even-flow of timber harvest volume subject to

adjacency, green-up, minimum harvest age, and spotted

owl habitat constraints. The objective function was designed

to measure and minimize the deviations of actual scheduled

harvest volume from some target harvest volume.

minimize ∑ ∑(xit vit) –TV

2

132

Figure 1—A 235,000 ha landscape in the Coast

Range of Oregon.

  

 

T 

t=1

N

i=1
Where:

t = a time period

T = total number of time periods

i = a management unit

N = total number of management units

TV = target timber harvest volume per time period 

(a constant)

xit = a binary (0,1) variable indicating whether or not management unit i is clearcut during time period t

vit = the volume available in management unit i during

time period t

In this research, the time horizon was assumed to be 

40 years, and consisted of eight 5-year time periods. There

were 39,500 management units in the study area, each containing approximately 100 pixels. It is possible that not all

of the pixels in a management unit could be harvested in

each time period. A number of factors may prevent harvest

of a pixel, including the following: it is represented as having a non-forest vegetation condition, its age is below the

minimum harvest age, it is located within a riparian area, 

it is reserved as a leave-tree clump. Therefore, vit for each

management unit only represents the available volume 

during that time period.

Constraints

There are five main constraints contained within this

planning process. First, we assume that each management

unit could only be scheduled for harvest once during the

time horizon.

∑xit ≤ 1 i

Second, each management unit, based on the average

age of its forested pixels, had to be above a minimum harvest age before a clearcutting activity could be scheduled.

Thus the average harvest age of each management unit during each time period was computed, then the set of potential

choices for each management unit was determined.

∑ (aj Agejt)

∑ (aj)

If AHAit ≥ MHA Then xit ∈ {0,1}

Else xit = 0

Where:

j = a pixel within a management unit

Ji = the set of forested pixels contained within management

unit i

aj = area of a pixel

Agejt = age of trees contained within pixel j during time

period t

AHAit = average potential harvest age of management unit 

i during time period t

MHA = the minimum harvest age for a set of trees within a

management unit

The third constraint indicates that clearcut sizes are limited to a maximum of 48.6 ha (120 acres). The scheduling

process used a recursive function to determine the actual

size of all proposed clearcuts by sensing the size of potential clearcuts from neighbors of proposed sales, and their

neighbors, and so on.

xit CAit + ∑xzt CAzt ≤ MCA i,t

Where:

CAit = potential clearcut area of management unit i during

time period t

z = an adjacent management unit

Ni = the set of all management units adjacent to management unit i

Si = the set of all management units adjacent to those management units adjacent to management unit i

MCA = the maximum clearcut area assumed (48.6 ha)

This definition is very similar to that proposed by

Murray (1999) for area restriction problems. Si is essentially

the set of clearcut units that contain all units adjacent to the

neighbors of unit i, and all units adjacent to the neighbors

of neighbors, etc. It is referred to as a recursive function

133

Figure 2—Age class distribution of forest land in a 235,000 ha landscape

in the Coast Range of Oregon.

Ji

Ji

T

j=1

j=1

t=1

A

A



 



z∈Ni∪Si

= AHAit i,t A
because of the potentially sprawling cluster group that must

be assessed depending on the direct or indirect contiguity

(of other clearcut units) with management unit i (Murray

1999).

The fourth constraint is related to the heterogeneity of

vegetation contained in the initial raster vegetation GIS

database. Given the salt and pepper nature of resulting

classifications of satellite images, we wanted to ensure that

a significant proportion of the forested area of each management unit was actually available for clearcut harvest

prior to allowing a clearcut to be scheduled. While this

constraint may seem redundant, given the minimum harvest

age constraint noted above, it does prevent the scheduling

of small (and very old) portions of management units. Here,

the percentage of forest pixels that are available for clearcut

harvest are summed each time period, and if this percentage is higher than a minimum percentage, the management

unit may be available for clearcut harvest. A minimum harvest area of 65% was chosen based on numerous evaluations

conducted within the CLAMS project. Higher minimum percentages significantly constrain harvest levels, while lower

percentages allow small parts of management units to be

harvested.

If γit ≥ MHP Then xit ∈ {0,1}

Else xit = 0

Where:

γit = the percentage of management unit i available for

clearcut harvest during time period t

MHP = the minimum harvest percent (minimum percentage area available for clearcut harvest within each

management unit)

Finally, the fifth constraint allows the scheduling model

to constrain activities to those that are possible without

reducing the “grow only” spotted owl HCI levels “too

much.” Here, the spotted owl HCI levels were first computed for the landscape under the assumption that no harvests will occur over the 40 year time frame. Then, a

maximum allowable decrease was determined for each

time period.

LHCIt ≥ w (“grow only” LHCIt ) t

Where:

LHCIt = the potential spotted owl HCI level for the entire

landscape during time period t

w = the maximum percentage decrease (0 ≤ w ≤ 1)

The landscape HCI levels are computed from individual

management unit HCI estimates

∑ (Ait HCIit)

∑ Ait

Where:

Ait = forested area of management unit i during time period

t

HCIit = the HCI estimate for management unit i during

time period t

We utilized the “model 2” spotted owl HCI model

developed by McComb and others (2002) to make the HCI

evaluations across the landscape. The spotted owl HCI

model was developed to measure the capability of a patch

and its surrounding neighbors to provide conditions important to survival and reproduction of spotted owls. HCI levels

were estimated by computing a nest stand capability index

and a landscape capability index:

HCIit = (NCIit2 LCIit )(1 / 3)

Where:

NCIit = the nest stand capability index component of the

HCI for management unit i during time period t

LCIit = the landscape capability index component of the

HCI for management unit i during time period t

NCIit is a function of the density of trees 10-25 cm dbh,

trees 25-50 cm dbh, and trees 75+ cm dbh, as well as a

diameter diversity index. LCIit is a function of the percentage of large and medium sized trees within three zones

around each management unit: 28 ha, 202 ha, and 1,810 ha.

These spatial zones were translated to radii around each

management unit: 300 m, 800 m, and 2,400 m. The centroid of each management unit was used as the base point

of each evaluation, and if any other management unit’s

centroid was within these search distances, they contributed to LCIit. This represents a simplification of the

process described in McComb and others (2002), where

LCI was computed as a moving window around each pixel,

more precisely measuring the zonal areas. Here, we computed management unit HCI levels, and by basing the LCI

computations on this, gain significant processing speed

during the scheduling process. The disadvantage of this

approach is that the difference in the evaluation of owl

134

A

A

N

N

i=1

j=1



 



= LHCIit t
habitat (by pixel or by management unit) reduces the heterogeneity possible in habitat values from 3,752,714 possible distinct values to 39,500. In addition, the spatial zone

evaluated using the original approach would resemble a

smooth circular radii evaluated around each pixel, while

the approach implemented here results in a rough approximation of the radii around each management unit. The

impact of these deviations from the original model, in

terms of the quality of habitat evaluations, is uncertain,

however, maps of the initial habitat condition values were

compared to those provided by McComb and others

(2002), and through visual inspection we determined that

the difference in the spatial distribution of habitat (within

the ranges mapped) was negligible even though it was not

entirely clear whether McComb and others (2002) used the

same data as used here.

Heuristic scheduling process

A threshold accepting heuristic technique called Habitat

Optimization for Owls (HOO) was developed and utilized

to schedule management units for clearcut harvest over the

40-year time horizon. Threshold accepting was introduced

by Deuck and Scheuer (1990), and has been applied to

forestry problems in Bettinger and others (2002, 2003).

Our implementation of threshold accepting to this problem

used a Model I (Johnson and Scheurmann 1977) problem

formulation, and 1-opt moves (i.e., a change to the timing

of harvest of a single management unit) to arrange efficient

harvest schedules for the landscape of interest. Given a feasible solution, a management unit and a harvest time period

were chosen at random, as shown in the flow chart describing the search process (fig. 3). This “potential” change to

the solution was then evaluated with respect to the constraints. If a constraint is violated, the process reverted to

the previous feasible solution, and a new potential choice

of management unit and harvest time period was made. If

no constraints are violated, the objective function was evaluated. If the proposed solution results in the overall best

solution found, the proposed solution became both the best

and the current solutions. If the proposed solution is not the

best solution, yet was within a certain range of the overall

best (as defined by the best solution value - the threshold

value), the proposed solution was retained as the current

solution, and the process continued. If the proposed solution

was not within a certain range of the overall best solution,

the process reverted to the previous feasible solution. The

solution process was run a number of times with several

target even-flow harvest volumes. The preliminary results

provided here illustrate the potential trade-offs among

even-flow harvest volume and spotted owl habitat capability index levels.

RESULTS

When average harvest volumes for the eight planning

periods are plotted against the average spotted owl habitat

capability index levels, one can begin to visualize the

trade-offs among timber production and habitat levels for

this landscape (fig. 4). In this case, the spotted owl HCI

levels were unconstrained. If we were to constrain spotted

owl HCI levels to so that they may decrease x% from the

grow-only projections of HCI levels (for the entire landscape), a clearer picture of the interaction between the

commodity production and ecological goals arises (fig. 5).

Here, we see that allowable decreases in HCI of 20% or

135

Figure 3—A flow chart of the threshold accepting scheduling process.
more had no effect on timber harvest levels. When HCI

was constrained to an allowable decrease of less than 20%,

there was almost a linear relationship between the allowable decrease and the resulting even-flow timber harvest

levels. Similarly, average HCI levels decreased in a nearly

linear fashion as the allowable decrease in HCI from the

grow-only condition increased, up to 20% (fig. 6). Allowable

decreases in HCI beyond 20% seemed to have no effect on

average HCI levels on this landscape. The HCI level at this

point (0.168) represents the level most likely under full 

utilization of the landscape under economic objectives.

Allowing small deviations in owl habitat levels from the

“grow only” case constrains the achievement of economic

objectives to a point (up to a 20% allowable decrease), but

beyond this point this constraint has no impact on solution

values - HCI levels are at their lowest given the age class

distribution of the landscape, the management direction,

and the other constraints recognized (e.g., adjacency).

Figures 5 and 6 contribute to this insight. 

While we have been relating timber harvest levels to

average HCI levels over the 40-year planning horizon, HCI

levels were not static, and generally increased over time

(fig. 7). The trend in the increase of HCI levels over time

seemed to be consistent no matter what target even-flow

harvest level was specified in the HOO model. The ‘nocut’ HCI levels in figure 7 were the levels obtained when

the landscape was projected under a grow-only policy. It is

from these levels that the allowable decrease in HCI was

measured when HCI is constrained.

136

Figure 4—Production possibility frontier relating average HCI level to 

average even-flow timber harvest level, when HCI levels are unconstrained.

Figure 5—Even-flow harvest levels when HCI is allowed to decrease from

1 to 40 percent from grow-only projections of the landscape.

Figure 7—Actual HCI levels over the planning horizon under three different policies regarding the desired even-flow timber harvest levels (i.e., the

target harvest levels).

Figure 6—Average HCI levels when HCI is allowed to decrease from 1 to

40 percent from grow-only projections of the landscape.
Spotted owl HCI levels can also be mapped, as the HCI

for each management unit is reported by HOO. It is our

intent to perform some spatial analyses of these maps to

investigate the spatial effects of various timber harvest

level goals, although those analyses are not provided here.

Figure 8 illustrates the change in HCI levels from the initial condition of the landscape to the end of the time horizon, when the target even-flow harvest level was 20

million board feet per year.

DISCUSSION

The HOO model can be viewed as an optimization

model that allows the evaluation of different policies containing both economic and ecological goals, and one that

can be applied across a broad area and over a long time

frame. This ability, however, is not unique, as other heuristic models have also been described in the literature where

complex habitat relationships are incorporated into a mathematical programming context. HOO merges strategic and

tactical planning considerations, so that processes at a variety of temporal and spatial scales can be represented. We

have shown here the use of the area restriction model for

handling adjacency relationships of almost 40,000 management units (resulting potentially 320,000 integer variables),

which may be intractable in other mathematical optimization processes.

The quality of solutions generated by a model such 

as HOO can be assessed in three general ways. First, one

could compare the results provided here to solutions generated by formulating and solving the problems with linear

programming (LP) techniques. These LP results would necessarily should be viewed as “relaxed” results, as one or

more constraints (i.e., adjacency and owl habitat) would

need to be ignored to allow the LP model to solve the problem. However, if produced, the LP results could be viewed

as an upper bound on potential solutions to the problem.

Second, the problems could be formulated and solved with

integer programming (IP) techniques, resulting in direct

comparisons with the heuristic solutions, since the constraints may be formulated in such as way that they can all

be recognized. Given the number of potential integer variables, however, this seems an unreasonable goal. Finally,

one could generate an estimated global optimum solution

from a set of heuristic solutions, as proposed by Bettinger

137

Figure 8—Spatial representation of the HCI levels for the initial condition of the landscape, and the ending condition

(after the 40-year planning horizon) when the target 
et al. (1998) and Boston and Bettinger (1999). This latter

approach has its drawbacks, including the time required 

to generate a sufficient set of solutions, the assumptions

required to confidently use extreme value theory, and the

fact that a bad heuristic can produce a misleading estimate

of the global optimum solution to a problem. For these 

reasons, we leave this work for the future.

Other capabilities are being integrated into the HOO

model beyond those described here. For example, since

vegetation is represented at the pixel level, the ability to

restrict management within certain distances of the stream

system is possible. In addition, the ability to model leaving

“leave tree clumps” in clearcut areas is possible. We plan

to, given a desired level of retention specified by a user of

HOO, randomly leave x% of potential clearcut pixels uncut

in harvest units. This would facilitate analyses of potential

legacy policies as well as voluntary ecological-oriented

contributions of landowners. HOO also facilitates an examination of a variety of green-up and adjacency policies, as

the length of the green-up period and the size of the maximum clearcut is not fixed, but variable, and adjustable by

users of the model.

Other management activities have been suggested as

ways to accelerate the development of spotted owl habitat

in the Pacific Northwest. Thinnings, in particular, may

result in faster growth rates of residual trees, thus accelerated development of spotted owl habitat. However, given

the short tenure of the planning horizon (40 years), it remains

to be seen whether the potential gains through thinning can

be realized via this modeling process. Thinnings could,

however, have an effect on total timber harvest volume

produced, thus could potentially increase the maximum

even-flow harvest levels. Thinnings are, however, less

widely used in steep terrain, as the cost of logging may

exceed the revenues generated, so it remains to be seen

whether the impact of thinnings is significant on timber

harvest levels (i.e., more than 10 percent, which is a

common assumption about how much thinning volume

contributes to overall harvest levels in the Pacific Northwest).

LITERATURE CITED

Bettinger, P.; Graetz, D.; Boston, K.; Sessions, J.; 

Chung, W. 2002. Eight heuristic planning techniques

applied to three increasingly difficult wildlife planning

problems. Silva Fennica. 36: 561-584.

Bettinger, P.; Johnson, D.; Johnson, K.N. 2003. Spatial 

forest plan development with ecological and economic

goals. Ecological Modelling. 169: 215-236.

Bettinger, P.; Sessions, J.; Johnson, K.N. 1998. Ensuring 

the compatibility of aquatic habitat and commodity production goals in eastern Oregon with a Tabu search procedure. Forest Science. 44: 96-112.

Bettinger, P.; Sessions, J.; Boston, K. 1997. Using Tabu 

search to schedule timber harvests subject to spatial

wildlife goals for big game. Ecological Modelling. 94:

111-123.

Boston, K.; Bettinger, P. 1999. An analysis of Monte 

Carlo integer programming, simulated annealing, and

tabu search heuristics for solving spatial harvest scheduling problems. Forest Science. 45: 292-301.

Calkin, D.; Montgomery, C.A.; Schumaker, N.H.; 

Polasky, S.; Arthur, J.L.; Nalle, D.J. 2002.

Developing a production possibility set of wildlife

species persistence and timber harvest value using simulated annealing. Canadian Journal of Forest Research.

32: 1329-1342.

CLAMS (Coastal Landscape Analysis and Modeling 

Study). 2003. Home page. http://www.fsl.orst.

edu/clams/ (accessed 11/5/03).

Deuck, G.; Scheuer, T. 1990. Threshold accepting: A

general purpose optimization algorithm appearing superior to simulated annealing. Journal of Computational

Physics. 90: 161-175.

Hof, J.; Bevers, M.; Joyce, L.; Kent, B. 1994. An integer 

programming approach for spatially and temporally

optimizing wildlife populations. Forest Science. 40:

177-191.

Hof, J.G.; Joyce, L.A. 1992. Spatial optimization for 

wildlife and timber in managed forest ecosystems.

Forest Science. 38: 489-508.

Hof, J.G.; Raphael, M.G. 1993. Some mathematical 

programming approaches for optimizing timber ageclass distributions to meet multispecies wildlife population objectives. Canadian Journal of Forest Research.

23: 828-834.

Johnson, K.N.; Scheurman, H.L. 1977. Techniques for 

prescribing optimal timber harvest and investment

under different objectives - discussion and synthesis.

Forest Science. Monograph 18. 31 pp.

138
McComb, W.C.; McGrath, M.T.; Spies, T.A.; Vesely, D. 

2002. Models for mapping potential habitat at landscape

scales: An example using northern spotted owls. Forest

Science. 48: 203-216.

Murray, A.T. 1999. Spatial restrictions in harvest scheduling. Forest Science. 45: 45-52.

Ohmann, J.L.; Gregory, M.J. 2002. Predictive mapping 

of forest composition and structure with direct gradient

analysis and nearest neighbor imputation in coastal

Oregon, USA. Canadian Journal of Forest Research. 

32: 725-741.

139
This page is intentionally left blank.
In: Bevers, Michael; Barrett, Tara M., comps. 2005. Systems Analysis in Forest Resources: Proceedings of the 2003 Symposium. General Technical Report

PNW-GTR-656. Portland, OR: U.S. Department of Agriculture, Forest Service, Pacific Northwest Research Station.

1 Bruce Meneghin is GIS Analyst with USDA Forest Service, Inventory and Monitoring Institute, Fort Collins, CO, USA. 2 Allison Reger and Neal 

Forrester are operations research analyst and planner, on the Willamette National Forest, Eugene, OR, USA.

Systems Analysis in Forest Resources: Proceedings of the 2003 Symposium

ASSESSING SHORT- AND LONG-TERM 

THINNING OPPORTUNITIES ON THE 

WILLAMETTE NATIONAL FOREST

Bruce Meneghin1, Allison Reger2, and Neal Forrester2

INTRODUCTION

The Willamette National Forest is comprised of 1.6 

million acres of land on the western side of the Cascade

mountains in western Oregon, USA. From the 1960’s

through the early 1990’s the Willamette was a significant

timber producing forest, harvesting between 500 and 800

million board feet of timber annually. Since that time, management objectives of the forest have shifted dramatically.

The current management direction, established by the

Willamette Forest Plan (USDA Forest Service 1990) and

the Northwest Forest Plan (USDA Forest Service 1994),

emphasizes the establishment and maintenance of a system

of Late Seral Reserves (LSRs), and implementing an aquatic

conservation strategy for protecting riparian habitat and

native fish stocks. This management direction has divided

the Forest into three broad land allocation classes, shown

in figure 1. One third of the forest is allocated to LSRs

while another 14 percent is designated as riparian areas

critical to the implementation of the aquatic conservation

strategy. The remainder of the lands are referred to as

ABSTRACT

The Willamette National Forest conducted an analytical exercise to examine thinning treatments in the short and long

term. Objectives of the analysis were to verify short-term thinning estimates; to gain an understanding of long-term implications and capabilities of commercial thinning on the forest; to develop a model to be used for analyzing vegetation management proposals and policies in the future; and to provide useful information about thinning that could be linked with

short, mid and long-term management strategies. A linear programming model was constructed, and scenarios with varying

mixes of operational and allocation constraints were examined. Where practical, multiple entry prescriptions offered the

best prospect of achieving a sustained flow of forest products. Limiting thinning to Late Seral Reserve (LSR) lands or

stands previously managed greatly reduces available harvest volumes in the short term. Issues that need further investigation

are fuels treatments and management objectives on lands outside of LSRs.

141

matrix lands – a mix of vegetation types and land allocations that surround the LSRs and riparian areas. 

Under the current management direction the harvest

levels are projected to be 116 million board feet annually.

However, timber harvests in mature and old-growth forests

using clearcutting are very controversial on the Willamette.

Therefore, the most promising opportunity for achieving

the planned harvest levels is to conduct thinnings in the

350,000 acres of young forest that were created during the

last forty years. 

STUDY OBJECTIVES

This analysis was driven by four key questions from

managers:

• What would the maximum harvest level be if thinning

was the only harvest activity?

• What would the maximum harvest level be if thinnings were limited to lands that had been previously

managed?
• How long could a “thinning only” program be 

sustained?

• What would be the long-term implications of limiting 

harvest activity to thinning?

In order to answer these questions an analysis project

was designed with these analytical objectives:

• Verify thinning estimates made by District offices.

Crude estimates based on acreages had been developed at the field level. We wanted to improve the

accuracy of these estimates.

• Determine the long-term implications of a commercial 

thinning program. We needed to be able to quickly

and consistently estimate future forest age class distributions under different management scenarios.

• Provide information useful to short-term and longterm thinning strategies.

• Develop a model for future analyses. We wanted to 

develop a model that could be easily adapted to analyze questions posed by future administrative and

political proposals.

MODEL DEVELOPMENT

To meet these analysis objectives we decided to construct a fairly traditional linear programming harvest 

schedule model with a Model 1 structure as presented by

Johnson and Scheurman (1977). The model was formulated

using Spectrum, the LP based planning model supported by

the USDA Forest Service (1998) and commonly used in

forest plan revisions (Rupe 1995). The contents and structure of the model are described briefly below.

Time Frame

The planning horizon was delineated by thirty, five-year

planning periods. Five-year periods was chosen to synchronize with the yield simulator, and a 150 year planning horizon was chosen to observe long term implications of

management strategies.

Land Stratification

Land attributes used for stratification, illustrated in figure 2, were chosen to differentiate current timber inventory,

growth rates and response to management, and to allow

expression of current and potential management objectives

in the model. Watersheds were included to provide geographic specificity and to control activities within watersheds. The

land status layer was a reclassification of management designations that identified suitability for thinning prescriptions. The previous management layer was used to assign

appropriate thinning prescriptions from a silvicultural point

of view. Finally, species and age class stratifications were

carried to differentiate inventory, growth and response to

treatment. The combination of these layers resulted in 3900

analysis units with unique attributes. In many cases, analysis units were comprised of non-contiguous land units.

142

Figure 1—Proportions of the Willamette National Forest by major land

status groups.

Figure 2—Analysis units were created by overlaying five layers of information. A total of 3900 unique analysis units were used in the analysis.
Management options

A variety of thinning treatment options were developed

for each analysis unit. Within LSR areas, these treatments are intended to encourage the development of key

old-growth characteristics desired in LSRs (Spies 1991;

Willamette National Forest 1998). These included regimes

with and without pre-commercial thinning, and regimes

that had one, two or three commercial thinnings during 

the planning horizon.

Outputs and Effects

To keep the model simple, a small number of outputs

and effects were modeled. Timber growth and yield projections were developed using the Forest Vegetation Simulator

(FVS). These simulations provided estimates of inventory

and harvest amounts as well as estimates of fuels produced

from thinning operations. The per-acre cost of thinning

operations was also estimated. These costs varied by combinations of land status, species group and thinning intensity.

Management Objectives

Management objectives that were likely to intersect

with harvesting objectives were represented in the model.

These objectives and their representation in the model are

listed by resource area.

Fish and Watershed—

• Limit potential ground disturbance and possible sedimentation or impacts to streamside vegetation, amphibians,

and arthropod populations. 

• Ensure not all closed canopy dispersal corridors in a 

given watershed are impacted at the same time and limit

potential adverse cumulative effects within desired levels.

• Maintain or minimize disruption of microclimates in 

riparian reserves and minimize edge impact on streamside

stands.

• Minimize stream course impacts during peak flow events.

These objectives were represented in two sets of constraints in the model. The first controlled the number of

acres within riparian areas that could be harvested in a

period. Willamette biologists felt that limiting the harvesting to less than seven percent of the riparian acres per period would keep cumulative impacts to an acceptable level

(equation 6, below). The second set of constraints aimed 

at protecting the fish resource was the enforcement of standards established in the Northwest Forest Plan (USDA Forest

Service 1994). This standard limits the acres harvested in

selected watersheds to less than fifteen percent of the acres

in that watershed per period (equation 8).

Wildlife—

• Conduct silvicultural treatments in young stands within

LSRs to promote recruitment into suitable old-growth

habitat.

• Provide dispersal habitat for spotted owls as well as 

maintaining some level of landscape continuity and dispersal habitat for other species.

These objectives were modeled with decision variables

representing management choices designed to produce

desired habitat and with constraints that limited the number

of acres treated within LSRs. Specialists felt that limiting

the treatments to less than twenty five percent of an LSR in

a ten-year period would ensure that habitat requirements

could be met (equation 7).

SCENARIO DESIGN

In order to meet the analysis objectives, three broad

level scenarios were studied. Scenario 1, designed as a

benchmark, sought to maximize harvest levels with few

constraints. Scenario 2 was designed to represent the management objectives for fish, wildlife and watershed that

might interact with timber production objectives. Finally,

Scenario 3 was similar to Scenario 2 but harvest activities

were limited to stands that had been previously managed.

The mathematical formulation for Scenario 2, stated as in

Clutter et al. (1983), is shown below.

Maximize:

30 R

(1) Z = ∑ ∑ ∑ Vijk Xjk i=1 j=1 k∈LSR

2 R N

(2) Z′ = ∑ ∑ ∑ Vijk Xjk i=1 j=1 k=1

30 R N

(3) Z′′ = ∑ ∑ ∑ Vijk Xjk i=1 j=1 k=1

Subject to:

R N RN

(4) ∑ ∑ Vi+1jk Xjk ≤ ∑ ∑ 1.25Vijk Xjk (i = 1,2...29)

j=1 k=1 i=1 k=1

R N RN

(5) ∑ ∑ Vi+1jk Xjk ≥ ∑ ∑ 0.75Vijk Xjk (i = 1,2...29)

j=1 k=1 j=1 k=1

R R

(6) ∑ ∑ Hijk Xjk ≤ ∑ ∑ 0.05Xjk (i = 1,2...30)

j=1 k∈1 j=1 k∈RZ

143
R R

(7) ∑ ∑ Hijk Xjk ≤ ∑ ∑ 0.25Xjk (i = 1,2...30)

j=1 k∈LSR j=1 k∈LSR

R R

(8) ∑ ∑ Hijk Xjk ≤ ∑ ∑ 0.15Xjk (i = 1,2...30)

j=1 k∈WR j=1 k∈WR

where,

R = number of management regimes

N = number of analysis units

Vijk = volume per acre harvested from analysis unit k under

regime j in period i

Xjk = acres of analysis unit k allocated to management

regime j

LSR = set of analysis units residing in LSRs

Hijk = 1 if harvest occurs in period i on analysis unit k

under regime j, 0 otherwise

RZ = set of analysis units residing in riparian zone

WR = set of analysis units residing in watersheds with 

special restrictions

The model was solved using a lexicographic goal programming approach described by Mendoza (1987) and

commonly used in Forest Service planning models (Rupe

1995). The objective functions (1), (2) and (3) were solved

iteratively. Objective function (1) seeks to maximize harvest volumes produced in designated LSRs. After maximizing using equation (1), the model was solved to maximize

the volume harvested in the first two periods, (equation

(2)) subject to (4) through (8) and equation (1) reformulated as a constraint:

30 R

(9) ∑ ∑ ∑ Vijk Xjk>=0.99Z

i=1 j=1 k∈LSR

where Z is the objective function value achieved during the

first solution.

This process was repeated for objective function (3)

which maximizes volume produced over the entire planning horizon.

Equations (4) and (5) represent the maximum periodic

increase and decrease allowed for timber volumes. Equation

(6) meets management objectives for fish by limiting harvests in riparian areas, while equation (7) enforces wildlife

objectives by limiting harvests within LSRs. A series of

equations represented by (8) limited harvesting in select

watersheds.

RESULTS

All three scenarios showed similar trends in harvest volume through time as shown in figure 3. All scenarios were

characterized by steady or increasing harvest volume for

the first five decades followed by a steady decrease after

decade five. The lower volumes in the early periods of scenario 2 were due to the increased restrictions on percentage

of acres harvested in riparian and LSR areas in a given

period. These restrictions were compounded by a smaller

available land base in scenario 3 which resulted in an even

lower harvest in the early periods.

Scenarios 2 and 3 were able to match the volume production of the benchmark in decade five because of two

factors. First, the solutions for these scenarios chose regimes

that called for two and three entries on many acres. Secondly,

these solutions harvested more acres from matrix lands

than did the benchmark, as shown in figure 4. Matrix lands

allow greater use of the multiple entry regimes and thus

play an important role in maintaining harvest levels while

meeting other resource objectives. The slight increase in

harvest volume in decade eight for scenarios 2 and 3 comes

from a third thinning on matrix lands.

144

Figure 4—Proportion of thinned acres from different land classes in the

three scenarios.

Figure 3—Annual harvest volumes through time for the three scenarios.
All three scenarios show a transition in the long term 

to a forest that is greater than 95% late successional stands.

For all scenarios, the decline in acres available for thinning

begins in decade five and continues. Natural disturbance,

not accounted for in the model, will be the only means for

regenerating the forest.

CONCLUSIONS

This analysis shows a potential for the Willamette

National Forest to sustain a “thinning only” harvest program for seventy to eighty years. During that time harvest

levels climb from 65 mmbf on 4500 acres to 91 mmbf on

6200 acres annually for scenario 2. If harvests are limited

to managed stands, then early harvest levels fall to 40 mmbf

on 3200 acres annually. If LSR lands are maintained and

the aquatic conservation strategy objectives are met, the

greatest opportunities for harvest come from multiple-entry

thinning regimes on matrix lands. 

In addition to answering the questions outlined at the

beginning of the study, this analysis also raised new questions. The large proportion of thinning from matrix lands

called for by scenarios 2 and 3 will require that clear silvicultural objectives be established for these treatments. Some

mid-level analysis may also be needed to determine thinning priorities and entry schedules in sensitive watersheds.

Finally, additional analysis will surely be needed to examine the potential fire risks associated with fuels generated

from thinning operations.

LITERATURE CITED

Clutter, J. L.; et al. 1983. Timber management: a quantitative approach. New York: John Wiley & Sons. 333 p.

Johnson, K.N.; Scheurman, H.L. 1977. Techniques for 

prescribing optimal timber harvest and investment

under different objectives – discussion and synthesis.

Forest Science Monograph 18. 31 p.

Mendoza, G.A. 1987. Goal programming formulations 

and extensions: An overview and analysis. Can. J. For.

Res., 17: 575 – 581.

Rupe, J.B. 1995. The analysis of competing policy choices 

for the revised Black Hills forest plan: Ecosystem

restoration or sustained timber harvest levels? In:

Thompson, J.E. (compiler). Analysis in Support of

Ecosystem Management; Analysis Workshop III.; 1995

April 10 – 13; Fort Collins, CO: USDA Forest Service:

241 - 250.

Spies, T.A.; Franklin, J.F. 1991. The structure of natural 

young, mature, and old-growth Douglas-fir forests in

Oregon and Washington. In: Wildlife and vegetation of

unmanaged Douglas-fir forests. Ruggiero, L.F., KB.

Aubry, A.B. Carey, M.H. Huff, Tech Coordinators.

USDA For. Serv. Gen. Tech. Rep. PNW-GTR-285. 

pp. 91-109.

USDA Forest Service, Inventory and Monitoring 

Institute. 1998. Spectrum Users Guide.

http://www.fs.fed.us/institute/PAG/SPECTRUM.shtml

USDA Forest Service and USDI Bureau of Land 

Management. 1994.Record of Decision for Amendments

to Forest Service and Bureau of Land Management

Planning Documents within the Range of the Northern

Spotted Owl. USDA Forest Service PNW Region, P.O.

Box 3623, Portland, OR 97208.

USDA Forest Service, Willamette National Forest. 1990. 

Land and Resource Management Plan. http:/www.fs.

fed.us/rc/willamette/manage/forestplan/ index.html

145
This page is intentionally left blank.
In: Bevers, Michael; Barrett, Tara M., comps. 2005. Systems Analysis in Forest Resources: Proceedings of the 2003 Symposium. General Technical Report

PNW-GTR-656. Portland, OR: U.S. Department of Agriculture, Forest Service, Pacific Northwest Research Station.

1 Susan King is an Operations Research Analyst and 2 Brett Butler is a Research Forester at the USDA Forest Service; Northeastern Research Station; 

11 Campus Blvd., Suite 200; Newtown Square, PA 19073.

Systems Analysis in Forest Resources: Proceedings of the 2003 Symposium

GENERATING A FOREST PARCELIZATION

MAP FOR MADISON COUNTY, NY

Susan L. King1 and Brett J. Butler2

INTRODUCTION

A parcel is a contiguous polygon of land owned by a

single ownership. An ownership is an individual, group of

individuals, business, or public agency that has legal title to

the polygon of land. Between 1994 and 2002, the number

of private landowners in the United States increased by 7

percent (Butler 2004). At the same time, the area of forest

land collectively owned by these families, individuals, companies, and other groups increased by less than 3 percent

(www.fs.fed.us/woodlandowners). This continuing trend of

forest land ownership increasing faster than forest land is

resulting in the fragmentation of forest ownerships or forest

parcelization. An alternative way to define parcelization is

a decrease in the average landowner’s landholdings.

ABSTRACT

It is well documented that the size of forested landholdings directly affects the management of forest resources. Traditionally, this information has been presented in tabular format. A spatial format will provide a more useful mechanism for

summarizing and utilizing this important information and can provide new insight into the causes, or at least correlates,

with differences in sizes of forest landholdings. The objective of this study was to investigate techniques for modeling forest parcel size in order to facilitate future efforts to generate a national map depicting the distribution of forest parcel sizes.

County tax records and land-cover data for Madison County, New York, were used to generate a dataset containing forest

parcel size by owner. This area was selected because, unlike most of the country, the land ownership parcel data are readily

available in a format conducive to spatial analysis. Ancillary data included both social and biophysical predictor variables.

Multiple linear regression, neural network, and decision trees were used to model these data. The models were tested using

validation datasets and accuracy was assessed with enhanced kappa statistics. A decision tree was selected as the final

method because it provided superior validation results, fewer parameters to adjust, and ease of explanation and visualization

to non-experts.

KEYWORDS: Parcelization, decision trees, simple kriging with a varying local mean, neural networks.

147

The consequences of parcelization are not fully understood, but they have important implications for the preservation of working forest land (Sampson 2000). A basic

limitation imposed by parcelization is related to economies

of scale (Row 1978). For working forest land to be sustained, it is necessary for the land to be profitable, and timber is the most profitable market commodity from such

tracts. The profitability of a timber harvest is a function of

the realized stumpage value and fixed and variable costs.

While the variable costs can be adjusted, at least partially, 

to reflect the size of the parcel being treated, the fixed costs

(e.g., transportation of equipment) by definition are invariant to treatment area. 

Forest parcelization, itself is not detrimental to ecosystem functions, though it often is a catalyst for change.
Specifically, parcelization often precedes the conversion of

forest land to other land uses, e.g., residential development,

which, in turn, leads to fragmentation of the physical forest

resource. Forest fragmentation has been shown to be detrimental to ecosystem functions ranging from water quality

to the availability of nesting habitat (Rochelle et al. 1999;

Saunders et al. 1991).

Likewise, parcelization may not directly affect the

resources available to people desiring access (e.g., hunters),

but access will be hampered by the need to contact an

increasing number of land owners. This access also affects

programs designed to educate or influence landowners. As

the number of owners increases, the cost of implementing

cost sharing or extension programs will be affected. 

The characteristics of the owners also are correlated

with the size of landholdings (Butler 2004; Birch 1996).

Landowners of smaller forested tracts are more likely to

have objectives related to home ownership and aesthetics,

whereas owners of larger forested tracts are likely to

actively manage their forests for timber production.

To accurately predict forest parcelization patterns, the

processes or patterns of the underlying forces must be

understood. To date, most of the information on forest

parcelization has been obtained from coarse-scale ownership surveys that allow only tabular summaries, limited

spatial resolution, and limited ability to combine with other

data sources. Mehmood and Zhang (2001) built a multiple

linear regression model to describe parcelization for the

United States. Their dependent variable, a measure of

parcelization, is the change in average nonindustrial private

forest size between 1978 and 1994 on a per-state basis.

Their independent variables are suited for a state-level

model and reflect changes over time. Our model is a geographical approach; that is, we investigated methods that

model forest parcelization both on a continuous and categorical scale. The purpose of this case study was to gain

experience using different modeling techniques and independent variables which will enable us to build a national

map that depicts the size of forest parcels. This map will

148

Figure 2—Agricultural and forest land are the predominate land uses. An

urban corridor exists along Interstate 90 at the edge of the Lake plains

physiographic region. Syracuse is a major population center to the west.

Figure 1—Madison County, located in

central New York, is bordered to the

north by Oneida Lake. The Lake Plains

physiographic region surrounds Oneida

Lake; the southern portion of the county

is part of the Appalachian Uplands physiographic region. Numerous lakes are

interspersed among the rolling hills.
allow us to increase our understanding of the processes

underlying forest parcelization and its consequences.

METHODS

Study Area

The study area was Madison County in the center of New

York State (fig. 1). The county contains numerous lakes

interspersed among the flat plains and rolling hills of the

Lake Plains and Appalachian Uplands physiographic regions.

Agricultural and forest lands are the major land uses (fig.

2), but, as in many parts of the county, these industries

have been superseded by service and retail sectors in the

local economies. Currently, 69,000 people live in Madison

County (U.S. Census Bureau 2001). The county covers 661

mi2 giving it an average population density of 104 people

per mi2.

This area was selected because it is one of relatively

few counties in the United States where tax records 

have been digitized in a format that is amenable to further 

analyses and are publicly available for a nominal fee

(www.madisoncounty.org). All legal property boundaries

for this county have been digitized with polygon topology

and information about the ownership of each parcel is

available.

Variables

The dependent variable, size of forest parcel, was

derived by combining county’s tax records and remotely

sensed imagery. The (total landholdings) parcel map was

created by assigning every ownership in the county a

unique numerical identifier and then combining adjacent

parcels to remove “artificial” boundaries. A forest/nonforest map was created by reclassifying the 30m pixels in a

1992 National Land Cover Data (Vogelmann 2001) image

for Madison County as 1 for forest or 0 for all other landcover categories. The total number of forest pixels within a

parcel (based on pixel centroid) were calculated and multiplied by 0.2 ac to obtain the total forest area per parcel.

The independent variables for this analysis (table 1)

were limited to data that were publicly available and have

national coverage. Thus, the results can be used to create a

product that covers a larger geographic area. The physical

variables originated from remotely sensed imagery, digital

elevation models (DEM), and census data. Percent forest,

agriculture, and urban were derived from the 1992 National

Land Cover Data. Slope and elevation were estimated from

the DEM (Defense Mapping Agency 1997).

The Census Bureau (U.S. Census Bureau 2001) was the

source of the socioeconomic variables available at the census block or block group level. Data for the finest spatial

resolution available were incorporated. Variables related to

the number of people across the landscape included population, dwelling, and family density. These variables are all

highly correlated, but not completely redundant. The average income and education of the people within a census

block group were tested. A dummy variable indicating public ownership was derived from county tax records, though

this variable can also be derived from regional or national

data sources (McGhie 1996).

149

Table 1—Dependent and independent variables used in the study.

Variable Abbreviation Description

Forest parcel ForParc Area (acres) of forest in a parcel

Percent forest PerFor Percent forest within (1, 5, 25, or 100) acres of the parcel centroid

Percent agriculture PerAg Percent agriculture within (1, 5, 25, or 100) acres of the parcel centroid

Percent urban PerUrb Percent urban within (1, 5, 25, or 100) acres of the parcel centroid

Population density PopDen People per square mile within a census block

Dwelling density DwDen Dwellings or houses per square mile within a census block

Family density FamDen Families or households per square mile within a census block

Public ownership PubOwn Dummy variable indicating parcel centroid is within a public ownership

Distance to roads DistRoad Distance (mi) of parcel centroid to nearest improved road

Distance to water DistWater Distance (mi) of parcel centroid to nearest census water

Slope Slope Average slope within (5 or 100) acres of parcel centroid

Elevation Elev Average elevation within (5 or 100) acres of parcel centroid

Income Inc Average annual income (US$) within a block group

Education Ed Percent of population with a college degree within a census block group

Gravity Index Grav Weighted distance of the census block centroid to population centers
A gravity index (Kline 2001; Haynes 1984) was included

to represent the relative influence of population centers.

Parcels closer to an urban area are more likely to be subject

to urban land-use influences. It is hypothesized that this will

result in smaller parcels. The index was calculated as:

3 (Populationk)0.5

Gravity indexi = ∑

k=i Distanceik

where

i = U.S. Census block group centroid

k = a city having a large urban influence

The selection of three nearest cities was based on the

work of Kline and Alig (2001). Cities within an 50 mi

buffer of the county were included. The city with the

largest population (163,860 people) was Syracuse, located

west of the county. The remaining urban centers ranged in

population from 166 to 68,637. 

MODELS

Several procedures for mapping parcelization were

investigated. Both simple kriging with a varying local mean

and decision trees were investigated. In simple kriging with

a varying local mean, the large-scale variability is modeled

using a nonspatial approach such as multiple linear regression, neural networks, or generalized additive models. A

linear model, multiple linear regression, and a nonlinear

approach (neural networks) were investigated. Forest parcelization is predicted on a continuous scale with simple kriging and on a categorical scale with decision trees.

Multiple Linear Regression vs. Neural Networks

In multiple linear regression, the relationship between

independent and dependent variables is assumed to be linear in the unknown parameters and interactions between

the independent variables must be specified a priori by the

user. The parameters are estimated by least squares. In neural networks, there is no assumed relationship between the

independent and dependent variables. The relationship

between independent and dependent variables is learned

through an iterative process. 

A neural network requires no assumptions about the 

distributions, mean, or correlation of the errors. Neural 

networks may provide a superior solution over a traditional

statistical procedure when the distributions are unknown

and possibly nonlinear, outliers may exist, and noise is

present in the data (Burke 1991). In this research, a feedforward neural network with backpropagation, the most

widely used learning algorithm, is implemented. Briefly, a

backpropagation network consists of processing elements

analogous to the brain’s neuron. The backpropagation network is hierarchical, consisting of at least three layers. 

The first layer is an input layer with the same number of

processing elements as independent variables, which are

thought to relate to the output classification signal. This

layer typically operates as a simple input buffer and uses a

linear transfer function to distribute each of the input attributes to the second layer processing units. The second layer

is the hidden layer. There can be multiple hidden layers,

but usually only one is used. The number of processing

units in the hidden layer is undeterminable a priori. Each of

the hidden layers uses some form of a nonlinear sigmoidal

transfer function, which usually has the same number of

weights as the number of input attributes. The third and

usually final layer, the output layer, typically has the same

number of processing units as classification categories. 

The output signals from each of the hidden layer’s transfer

functions feed into another sigmoidal function at the output

layer. The difference between the observed and estimated

values of the dependent variable is formulated as a function

of weights. This is an unconstrained minimized problem

and there are numerous solution techniques available (King

1999). Selection of initial weights can influence the convergence to an optimum solution and ultimately the final

set of weights which are used on new data.

Simple Kriging with a Varying Local Mean

Because each datum has a geographical reference, 

geostatistical procedures are feasible. One geostatistical

approach for incorporating secondary information is simple

kriging with a varying local mean. Large-scale variability

is modeled using multiple linear regression or neural networks; small-scale variability is modeled using simple

kriging (Majure et al. 1996; Metzger 1997; Hunner et al.

1998; Goovaerts 2000; Lister et al. 2000). The two surfaces

are then added to produce the final map. If there is no smallscale variability, the large-scale map may stand alone.

Decision Tree

Decision trees partition the decision space into smaller

segments defined in terms of the input variables. These

partitions define the predictive relationship between the

inputs and the dependent variable. This partitioning continues until the subsets cannot be partitioned further using

user-defined stopping criteria. The homogenous groups

allow for greater accuracy in predicting the dependent variable from the independent variables. Chi-Square Automatic

Interaction Detector (CHAID), the decision-tree algorithm

used in this study, selects the most significant predictor

with the lowest adjusted probability at each level based on

the Bonferroni adjusted Chi-Square statistics (KASS 1980).

A small p-value indicates that the observed association

150
between the predictor and dependent variable is unlikely to

have occurred solely due to sampling variability. If a predictor has more than two categories, there may be many

ways to partition the data set based on the categories. A

combinatorial search algorithm finds a partition that has the

smallest p-value. CHAID partitions the data into multiway

splits, which result in more interpretable trees than binary

splits. The root of the tree is the entire data set. The subsets

and subsubsets form branches of the tree. Subsets that meet

a stopping criterion and cannot be further partitioned are

leaves. As with neural networks, over and underfitting are

modeling concerns in decision trees. A tree with too many

branches overfits the data and a tree with too few branches

underfits the data. 

Statistical Tests

The level of agreement between truth and the model

classification were evaluated using Cohen’s kappa (κ)

(Congalton and Green 1999). Some agreement between

truth and the model classification would be expected by

chance. The κ analysis gives the representation of the proportion of agreement beyond that expected by chance. A κ

score of 1 indicates perfect agreement and a score of 0

indicates no more agreement than would be expected by

chance. A negative score indicates less agreement than

would be expected by chance, that is, disagreement. The κ

statistic is based on an equal weighting of the classes, which

assumes that all of the errors are of equal importance.

Weighting the cells closer to the diagonal more heavily than

those cells farther away takes into account near disagreements. The Fleiss-Cohen weighting scheme was used.

Pontius (2000, 2002) has argued that the standard κ index

of agreement is inappropriate for map comparisons. The

standard κ was developed for studies when the marginal

distributions are unknown. The goal of a map is to obtain

similar marginal distributions as the population. Hence,

Pontius believes that the standard method to compute the

expected proportion correct classification by chance usually

is inappropriate for classification schemes that specify both

quantity and location. Two maps are best compared by tabulating the cross classifications between truth and prediction in a j x j contingency table, where j is the number of

categories in each map. Quantification error occurs when

the quantity of cells in a category in one map differs from

the quantity of cells in the other map. Even if there is no

quantification error, location error can occur where the

location of a category in one map differs from the location

of that category in the other map. The criteria of Landis

and Koch (1977) were used to interpret the κ scores.

Model and Validation Data Sets

Madison County has 35,370 land parcels of which

19,903 are forested. The dependent variable, forest parcelization, was divided into four classes (table. 2) based on

work by Birch (1996). The overwhelming number of small

parcels and the scarcity of large parcels pose a modeling

challenge. Equal number of observations in each class does

not represent all possible dependent and independent variable combinations, especially in the small parcel class. We

used different dependent variable class percentages. The

final model included 2 percent of the first class, (smallest

parcels), 10 percent of the second class, and 70 percent of

the third and fourth classes. Many trials were conducted

with randomly generated data sets with these percentages.

The final model and validation data sets had 791 and

19,112 observations, respectively. 

Results and Discussion

Models with a continuous dependent forest-parcelization variable were built using both multiple linear regression and neural networks. Variogram analysis showed that

there was a significant spatial dependence. The small-scale

structure was modeled using simple kriging. For both multiple linear regression and neural networks, adding the

small-scale variability improves the R2, 0.462 for multiple

linear regression and 0.892 for neural networks. The multiple linear regression model included percent forest within

100 acres, ownership, gravity, distance to roads, and distance to water, while neural networks substituted dwelling

density and slope within 100 acres for gravity. Unfortunately,

both models fail to generalize on new data. For the validation data set, multiple linear regression underpredicts parcel size and neural networks overpredicts parcel size,

particularly small parcels. After adding the small-scale spatial variability to the validation data sets, both procedures

predict numerous parcel values as negative. A neural network has many parameters that require adjustment (initial

weights, number of hidden layers and the number of hidden nodes within each hidden layer, type of sigmoidal

transfer function at the nodes in the hidden layers and output nodes, and optimization techniques). Additional adjustments could have been made, and possibly a better solution

obtained. However, it is not easy to interpret the output

from a neural network because it is not apparent which

input variables contribute heavily to the outcome. Also,

many parameters require adjusting in producing a variogram and in simple kriging. Because the ultimate objective

is a national map, these procedures were abandoned in

favor of a decision tree, which is highly interpretable has

fewer parameters to adjust, and is easily implemented on a

production scale. 

151
152

Figure 3—Final decision tree for classifying forest parcelization. Distance to roads is important in classifying parcel sizes. The larger parcels are farther

from the road than smaller parcels. Ownership is important is classifying large parcels. Slope and percent forest are important in classifying small parcels.

Figure 4—Map of a) observed and b) predicted parcel classes.

The large parcel sizes are on the public land in the higher elevations in the southern part of Madison County and on private land

in the plains near Lake Oneida in the northern part of the county.

The decision tree misses the higher parcels on private land. Most

of the misclassifications are classified as an adjacent parcel size.
The final decision tree is presented in figure 3. Good generalization is achieved by limiting the depth of the tree and

selecting a lower bound on the number of observations

required for a leaf. The final model was based on both the 

κ statistics and expert interpretation. In all trials, the most

significant predictor and the first to enter the model was

distance to roads, the most highly correlated variable with

forest parcelization (ρ = 0.34). The smaller distance to roads,

the more likely the forest parcel will be small. Slope and

percent forest within 100 acres further predict smaller

parcels. Public versus private ownership is the most significant predictor of large parcels. Dwelling density could be 

a replacement for percent forest within 100 acres of the

parcel in the fourth level of the tree. The final model completely misclassifies parcels over 99 acres on private land.

Figure 4 presents both the observed and predicted maps.

The over 99 acre parcels in the northern part of the county

are on private land and are misclassified in the lower class.

The confusion matrices for the model and validation data

sets are given in table 3. The percentage of perfectly

matched cells or overall accuracy is larger for the validation than the model data set (79.48 versus 71.94 percent)

due to the large number of small parcels. The κ statistic

decreases from 0.57 with a 95-percent confidence interval

of (0.53, 0.62) for the model to 0.37 with a 95-percent confidence interval of (0.36, 0.38) for the validation data set.

This indicates fair to moderate agreement according to

Landis and Koch (1977). Similarly, the weighted κ statistic

is 0.73 with a 95 percent confidence interval of (0.70, 0.77)

for the model, indicating substantial agreement, and 0.41

with a 95 percent confidence interval (0.40, 0.43) for the

validation data set, indicating moderate agreement. Pontius’

153

Table 2—Distribution of forest parcels by size class.

Parcel size (acres) Frequency Percent

0 - 9 17,088 85.86

9- 49 2,563 12.88

49 – 99 201 1.01

> 99 51 0.26

Table 3—Agreement between observed and predicted forest parcels 

for the model and validation data sets.

Model

Observed parcel class

Predicted 

parcel class 0 – 9 9 – 49 49 – 99 > 99 Total

0 – 9 311 42 9 0 362

9 - 49 31 160 58 2 251

49 – 99 11 41 83 6 141

> 99 0 11 11 15 37

Total 353 254 161 23 791

Validation

Observed parcel class

Predicted 

parcel class 0 – 9 9 – 49 49 - 99 > 99 Total

0 – 9 13,832 2,336 550 8 16,726

9 - 49 379 1,324 595 14 2,312

49 – 99 5 20 31 4 60

> 99 2 3 6 3 14

Total 14,218 3,683 1,1822 29 19,112
distribution of correct classification and errors percentages

is shown in figure 5. The validation data set is 79 percent

accurate.

CONCLUSION

It is difficult to build a model that captures both dominant and rare events. Of the methods tried, decision trees

gave the best results. There are fewer parameters to adjust,

it is easy for an expert to tweak the model, and it is visually explainable. Neural networks are superior to multiple

linear regression for this problem and newer neural network software might have more options to prevent overfitting. Ownership is the predominate variable in differentiating large parcels. Distance to roads, slope, and percent

forest best partition smaller parcels.

Currently, many of the consequences of forest parcelization are only hypotheses. Although, these hypotheses are

built on a plethora of antidotal evidence, more work is

needed to quantify the connections. The generation of forestparcelization maps across broad geographic areas will

allow these questions to be better explored. The methods

developed in this study are an initial step in this process

and will be expanded to cover a larger geographic area.

Future efforts will focus on using point-based estimates of

forest parcel sizes to generate a map depicting forest parcel

size. These estimates will come from the National Woodland

Owner Survey (www.fs.fed.us/woodlandowners). Due to

the sampling technique (probability proportional to size),

the inverse problem will be experienced: large parcels will

outnumber smaller parcels. The techniques explored in this

pilot study will still be applicable. From the map we will

be able to graphically display the distribution of forest parcel sizes across broad geographic areas and explore the

causes and consequences of observed patterns.

LITERATURE CITED

Birch, T.W. 1996. Private forest-land owners of the United 

States, 1994. Resour. Bull. NE-134. Radnor, PA: U.S.

Department of Agriculture, Forest Service, Northeastern

Forest Experiment Station. 183 p.

Burke, L.I. 1991. Introduction to artificial neural networks 

for pattern recognition. Computers and Operations

Research. 18(2): 211-220.

Butler, B.J.; Leatherberry, E.C. 2004. America’s Forest 

Owners. Journal of Forestry, 102(7): 4-9. 

Congalton, R.G.; Green, K. 1999. Assessing the accuracy 

of remotely sensed data: principles and practices. New

York: CRC Press. 137 p.

Defense Mapping Agency. 1997. Digital map atlas. 90 

meter digital terrain model. Washington, DC: Defenses

Mapping Agency.

Goovaerts, P. 2000. Geostatistical approaches for incorporating elevation into the spatial interpolation of rainfall.

Journal of Hydrology 228:113-129.

Haynes, K.E.; Fotheringham, A.S. 1984. Gravity and 

spatial interaction models. Beverly Hills, CA: Sage

Publications. 88 p.

Hunner G.; Reich, R.M.; Mower H.T. 1998. Modeling 

forest stand structure using spatial statistics. In: Whiffen

H. J-H.; Hubbard, W.C., eds. Proceedings of the 2nd

southern forestry GIS conference; Athens, GA. Athens,

GA: Univ. of Georgia, School of Forest Resources: 

103-120.

Kass, G.V. 1980. An exploratory technique for investigating large quantities of categorical data. Applied

Statistics. 29 (1980): 119-127.

154

Figure 5—Of the 19,112 parcels in the validation data set, 79 percent are

classified correctly. Because there are four classes, 25 percent are classified correctly due to chance, 42 percent classified correctly due to quantity, and 12 percent are classified correctly due to location. The location

error is 7 percent and the error due to quantity is 13 percent.
King, Susan L. 1999. Neural networks vs. multiple linear 

regression for estimating previous diameter. In Stringer,

Jeffrey W.; Loftis, David L., eds. Proceedings 12th

central hardwood forest conference; Lexington, KY.

Asheville, NC: Gen. Tech. Rep. SRS-24. US Department

of Agriculture, Forest Service, Southern Research

Station: 159-166.

Kline, J.D.; Alig, R.J. 2001. A spatial model of land use 

change in western Oregon and western Washington.

Res. Pap. PNW-RP-528. Portland, OR: U.S. Department

of Agriculture, Forest Service, Pacific Northwest

Research Station, 24 p.

Landis, J.; Koch, G. 1977. The measurement of observer 

agreement for categorical data. Biometrics. 33: 159-174.

Majure J,N.; Cressie, N.; Cook D.; Symanzik J. 1996.

GIS, spatial statistical graphics, and forest health. In:

Third international conference on integrating GIS 

and environmental modeling. National Center for

Geographic Information and Analysis. Sante Fe, NM,

Santa Barbara, CA: National Center for Geographic

Information and Analysis. CD.

McGhie, R.G. 1996. Creation of a comprehensive managed areas spatial database for the conterminous United

States: summary project technical report (NASANAGW-1743). Tech. Rep. 96-4. Santa Barbara, CA:

National Center for Geographic Information and

Analysis. 53 p.

Mehmood, S.R.; Zhang, D. 2001. Forest parcelization in 

the United States: a study of contributing factors.

Journal of Forestry. 99(4):30-34.

Metzger, K.L. 1997. Modeling forest stand structure to 

a ten meter resolution using Landsat TM data. Fort

Collins, CO: Colorado State University, Department 

of Forest Sciences. 123 p. M.S. thesis.

Pontius, R.G., Jr. 2000. Quantification error versus 

location error in comparison of categorical maps.

Photogrammetric Engineering & Remote Sensing.

66(8): 1011-1016.

Pontius, R.G., Jr. 2002. Statistical methods to partition 

effects of quantity and location during comparison of

categorical maps at multiple resolutions. Photogrammetric

Engineering & Remote Sensing, 68(10): 1041-1049.

Sampson, R.N.; DeCoster, L. 2000. Forest fragmentation: 

implications for sustainable private forests. Journal of

Forestry. 98(3): 4-8.

Smith, W.B.; Miles, P.D.; Vissage, J.S. 2004. Forest 

resources of the United States, 2002. Gen. Tech. Rep.

NC-241 St. Paul, MN: U.S. Department of Agriculture,

Forest Service, North Central Research Station.

Rochelle, J.A.; Lehmann, L.A.; Wisniewsk. J. 1999.

Forest fragmentation: wildlife and management implications. Leiden, The Netherlands: Koninklijke Brill NV

325 p.

Row, C. 1978. Economies of tract size in timber growing. 

Journal of Forestry. 76(9): 576-582.

Saunders, D.A.; Hobbs, R.J.; Margules, C.R. 1991. 

Biological consequences of ecosystem fragmentation: 

a review. Conservation Biology. 5:18-32.

U.S. Census Bureau. 2001. Census 2000 summary file 1 

(SFI). Washington, DC. U.S. Department of Commerce;

Census Bureau.

Vogelmann, J.E.; Howard, S.M.; Yang, L.; Larson, 

C.R.; Wylie, B.K.; Van Driel, N. 2001. Completion of

the 1990s national land cover data set for the conterminous United States from Landsat Thematic Mapper data

and ancillary data sources. Photogrammetric Engineering and Remote Sensing. 67: 650-652.

155
This page is intentionally left blank.
In: Bevers, Michael; Barrett, Tara M., comps. 2005. Systems Analysis in Forest Resources: Proceedings of the 2003 Symposium. General Technical Report

PNW-GTR-656. Portland, OR: U.S. Department of Agriculture, Forest Service, Pacific Northwest Research Station.

1 James Merzenich is an analyst, U.S. Forest Service, Pacific Northwest Region, P.O. Box 3623, Portland, OR 97208. 2 Leonardo Frid is a systems 

ecologist, ESSA Technologies Ltd., #300-1765 West 8th Ave., Vancouver, B.C. V6J 5C6.

Systems Analysis in Forest Resources: Proceedings of the 2003 Symposium

PROJECTING LANDSCAPE CONDITIONS 

IN SOUTHERN UTAH BY USING VDDT

James Merzenich1 and Leonardo Frid2

INTRODUCTION

The Vegetation Dynamic Development Tool (VDDT) is

a simulation model used to project landscape-scale vegetative conditions over long time frames (Beukema and others

2003; Merzenich and others 1999). VDDT models the effects

of alternative levels of management treatments on vegetation as influenced by stochastic disturbances such as wildfire. 

Discrete states are defined in VDDT on the basis of a

vegetative cover type and structure class. In each model

run VDDT uses up to fifty thousand simulation units to

project change. Simulation units are initially assigned to

states based on the proportion of area contained in those

states. These simulation units then progress along timedependent successional pathways or change states in

response to probabilistically applied disturbances or management treatments. VDDT projects the proportion of units

ABSTRACT

The Vegetation Dynamics Development Tool (VDDT) is a state-transition modeling system that integrates the effects of

succession and natural disturbances such as wildfire with management treatments. VDDT is used to evaluate the response

of non-spatial indicators to management actions. This study used VDDT to project changes in vegetation in response to different management scenarios for a mixed forest and range watershed in southwest Utah. The vegetative composition of this

landscape has been significantly altered because of the effects of fire suppression and grazing. Pinyon-juniper acreage has

increased from historic levels while aspen, sage, and oak habitats have declined. Key environmental factors projected by

VDDT include the percent of area dominated by open shrubs, pinyon-juniper, and aspen, and the percent of dry-site conifer

acres at high risk of uncharacteristically severe wildfire. The analysis shows that an aggressive vegetation treatment schedule is needed to reverse current trends and restore ecological conditions and processes in this area. 

KEYWORDS: Landscape modeling, natural disturbance, simulation, forest dynamics, vegetative succession.

157

(or area) contained in each state and the levels of disturbances that may be expected. 

The U.S. Forest Service, Pacific Southwest Station, is

conducting a study comparing landscape-scale fuel treatment models at seven major locations in the United States

(Weise and others 2001). Two reports describing the use of

VDDT in conjunction with this study are available. In a

retrospective analysis, VDDT models were used to predict

stand replacement fires between 1937 and 1996 in Yosemite

national park (Arbaugh and others, 2001). In the Bitterroot

front of western Montana VDDT was used to estimate the

level of fuel treatments needed to restore ecological conditions and reduce the risk of uncharacteristic stand-replacement wildfires (Merzenich and others, 2001). This paper

presents the results of applying VDDT to the Beaver River

drainage, an area of mixed range and forest, in southwest

Utah.
The study area contains 350,000 acres (142,000 ha.) located primarily on the Fishlake National Forest and adjoining

Bureau of Land Management land (fig. 1). Elevations vary

from approximately five to twelve thousand feet (1500 to

3700 meters). The study area contains acreage associated

with nineteen of twenty-four VDDT models developed 

for potential vegetation types in southern Utah (Long and

Merzenich, 2004). The major vegetation types represented

by these models are Wyoming big sage (Artemesia tridentata var. wyomingensis), mountain big sage (Artemesia 

tri-dentata var. vaseyana), Gambel oak (Quercus gambelii),

ponderosa pine (Pinus ponderosa), montane fir (mix of

Pseudotsuga menziesii and Abies concolor), and Spruce/fir/

aspen (Picea engelmannii, Abies lasiocarpa, and Populus

tremuloides). 

A combination of wildfire suppression and livestock

grazing has significantly altered the vegetative composition

of the study area over the past 150 to 200 years (Campbell

and others, 2003). At lower elevations, a mosaic of grass

and shrublands has been replaced by vast areas of dense

shrubs and non-native annuals such as cheat grass and red

brome (Bromus tectorum; Bromus rubra). On more mesic

sites pinyon pine (Pinus edulis) and juniper (Juniperus spp.)

are replacing sage and perennial grasses. Mid-elevation

mixed pine and fir forests, once maintained by frequent

ground fires, now have unnaturally high levels of easily

combustible fuel and are at high risk of catastrophic wildfire (Bradley and others 1992). Quaking aspen stands,

dependent upon fire disturbance for regeneration and detrimentally affected by livestock browsing, are being replaced

by shade tolerant spruce and subalpine fir (Bartos and

Campbell, 1998).

The objectives of this study were to track the response

of key indicators, such as aspen, ponderosa pine and pinyon-juniper acres, to management scenarios that represent:

1) no management; 2) management at levels commensurate

with current funding; and 3) management at a higher level

believed necessary to restore ecological conditions.

METHODS

Historically fires were common within most vegetative

types in the Beaver River drainage with estimated fire frequencies ranging from every 5-25 years for ponderosa pine

stands to 50-80 years for spruce/subalpine fir stands

(Campbell and others, 2003). 

Data on the acreage burned by individual wildfires for

the modern period are available for this drainage for the

years 1970 to 2002. Because of active fire suppression and

the reduction of fine fuels due to grazing, there has been a

decrease in the acres burned during this modern period, as

compared to historic. The average area burned per year in

wildfires was about two hundred acres for the period 1970

to 1991 and two thousand acres for the period 1992 to 2002.

The burn data for this latter period (1992-2002) are believed

to be more reflective of current fuel and climatic conditions

and were used to estimate annual wildfire probabilities for

VDDT. Based on these data, the average fire return interval

is now 175 years with 0.57 percent of the area burned per

year.

The wildfire acreage burned per year varies based on

ignitions, climatic conditions, and resources available to

suppress fires. Options available in VDDT to vary wildfire

disturbance probabilities temporally, by type of fire year,

were not used in this analysis.

Most wildfires in the last decade have occurred in the

grass/shrub and woodland classes. In this analysis average

fire return intervals are estimated to be 150 years (0.0067

annual probability) for grass/shrub and woodland, 200 years

for dry forest, and 300 years for high elevation spruce/fir

and aspen areas. These fire return intervals and annual

probabilities correspond to an expected acreage burned per

year for the entire study area of about 2,000 acres. This is

consistent with the burn data for the period 1992 to 2002. 

Wildfire probabilities also vary according to the succession class. For example, on areas dominated by annual

grasses or high-risk conifer stands wildfire probabilities are

increased, while on low-risk stands they are decreased. For

each vegetative class both stand-replacement and low intensity (underburns and mixed severity) fires are modeled. 

158

Figure 1—Map of Utah. Shaded area represents Beaver River drainage.
Three scenarios were developed to estimate the effect 

of alternative levels of management treatments over a 50-

year time frame. This time frame was chosen to make this

analysis consistent with other analyses comparing landscape-scale fuel treatment models (Weise and others 2001).

The VDDT models for this area are actually designed to

allow projections for 300 or more years.

The following assumptions apply to all scenarios: active

wildfire suppression will continue; regeneration harvest is

used to promote aspen; partial harvests are used to reduce

fuels and promote early seral species; and livestock grazing

will continue at current levels. The three scenarios are

described as follows.

1) No_mgt: No management except fire suppression

and grazing

2) Full_mgt: Management at a level designed to restore 

and maintain ecological function as quickly as practicable. 

3) Current_mgt: Management at a level commensurate 

with current funding. The treatment level is assumed

to be one-third that of the full-management scenario.

An explanation of the full-management scenario (scenario 2) as it applies to the major vegetation types follows.

Table 1 shows projected treatment levels associated with

this scenario along with the estimated historic and current

fire return intervals. 

Grass and Shrublands: Approximately 42 percent 

of the study area was historically dominated by perennial

grass and sagebrush. Two major management issues in

these vegetation types are the proportion of acres in which

pinyon-juniper is dominant and the amount of open versus

closed shrub. Pinyon pine and juniper develop slowly and

are susceptible to mortality by wildfire. Historically pinyon-juniper stands were mostly confined to rocky outcrops

or other areas that frequently escaped wildfire (Bradley 

and others 1992). The area dominated by pinyon-juniper is

presently at least three times higher than the historic level

(O’Brien 1999). Pinyon-juniper stands reduce water yield

and streamflow through evapo-transpiration, and retard the

growth of grasses, forbs and shrubs. Pinyon-juniper stands

are often fire resistant due to the lack of fine fuels in the

understory. 

Most sage stands are presently closed canopied. Natural

sagebrush areas maintained by fire would have nearly equal

acreage in perennial grass, open shrub, and closed shrub. 

Few large wildfires presently occur in this landscape.

An average fire return interval of 150 years was applied to

areas dominated by perennial grass, shrubs and pinyonjuniper. On areas dominated by annual grasses (e.g. red

brome or cheatgrass) the assumed fire return interval is 50

to 75 years. If large areas of contiguous annual grasses

were to develop wildfire probabilities could dramatically

increase. In cheatgrass dominated landscapes in southern

Idaho average fire-return intervals are now less than 10

years (Paysen and others, 2000).

The primary prescription is to apply mechanical treatment (Dixie harrow, brush choppers, cut and burn) and prescribed fire to closed shrub classes and classes containing

pinyon-juniper at a ratio of 75 percent mechanical treatment

to 25 percent prescribed fire. This is done at a 50-year return

interval for Wyoming sage and a 30-year interval for mountain sage areas. These treatments include seeding annual

grass areas to perennial grass where appropriate with an

assumed success rate of 50 percent. The purpose of these

treatments is to remove pinyon-juniper and to lower the

density of sage. The closed shrub and pinyon-juniper classes younger than 200 or 250 years are the target areas to be

treated. Some non-target areas are burned with prescribed

fire, however. When applying prescribed fire we assumed

that perennial grass, annual grass, and open shrub areas

159

Table 1—Wildfire frequencies and treatment levels for Full_mgt scenario (Campbell and others, 2003)

Vegetation Fire freq. Fire freq. Trt freq

type (historic) (current) Treatment Ann rate Targeted classes

Wyo. Sage 40-60 150 Mech/Pres burn 50 (2%) All stands with pj or closed shrub

Mtn. Sage 20-40 150 Mech/Pres burn 30 (3.3%) All stands with pj or closed shrub

Oak woodland 20-50 150 Pres burn 50 (2%) All but seedlings

Dry forest 5-25 200 Pres burn 30 (3.3%) Pole and larger open

Thin/burn 50 (2%) Mature dense

Sprucefir/aspen 50-80 300 Pres burn 50 (2%) All mature/old

Regen harv 100 (1%) Mature sprucefir/aspen
will be burned at a ratio that is approximately one-third

that of the other classes. Thus, if targeted states were

burned at a rate of 3 percent per year, these non-target

states would be burned at a rate of 1 percent per year.

Woodlands: Oak and associated woodlands comprise

about 29 percent of the Beaver River landscape. Historically

fires occurred in these areas every 20 to 50 years (Campbell

and others, 2003). Pinyon-juniper now dominates a majority

of this acreage. The major issue is controlling the amount

of pinyon-juniper and maintaining young oak stands commensurate with the historic fire return interval. The assumed

current fire return interval for these areas is 150 years. The

primary treatment is to prescribe burn these areas at a rate

of 2 percent per year. 

Dry Forest: Approximately 10 percent of the study area

is comprised of dry forest comprised mostly of ponderosa

pine, Douglas fir, and white fir. These stands often have an

aspen associate. The historic fire return interval is 5 to 25

years (Campbell and others, 2003). Wildfires in dry forests

are now uncharacteristically severe and most areas that

were dominated by ponderosa pine and aspen are now

overstocked with white fir and Douglas fir. Major issues

are the decline of aspen and ponderosa pine. Fuel treatments are needed to restore the natural role of fire. The

assumed current fire return interval for these areas is 200

years. The primary treatment is to apply a salvage treatment (partial harvest and burn) to the mature and old highrisk classes at a rate of 2 percent of these areas per year.

This treatment should transform these stands to the lowrisk class. Immature, mature, and old low-risk stands are

then prescribed burned every 30 years to maintain them in

an open-canopied low-risk state.

High Elevation SpruceFir/Aspen forests: Approximately 12 percent of the study area is comprised of higher

elevation forests. Dominant cover types are subalpine fir,

Engelmann spruce and aspen. Historically fires occurred 

in these areas every 50 to 80 years, maintaining about 80 

percent of the area in aspen (Campbell and others, 2003).

Presently most acres are dominated by spruce and fir. The

major management issue is the decline in the acreage and

vigor of aspen. The primary treatments are to prescribe

burn the mature aspen and spruce/fir stands capable of supporting aspen at a rate of 2 percent per year and to apply

regeneration harvest to mature spruce/fir stands still containing an aspen component at a rate of 1 percent per year.

These treatments restore aspen and reduce the amount of

overage spruce/fir. In southwest Utah all aspen reproduction

is by vegetative sprouting usually triggered by disturbance

(Bartos and Campbell 1998).

RESULTS

Four key environmental variables are used to compare

the results of each scenario over the next 50 years.

1) Acres dominated by pinyon–juniper within both 

the grass/shrub and woodland models. Pinyon-juniper

acreage has increased at least three fold from historic

levels. The treatment objective is to reduce this acreage.

2) The acres of open versus closed shrub in grass/shrub

areas. A relatively even proportion of open and closed

shrub provides near optimal sage grouse habitat and is

more reflective of historic conditions. Presently nearly

90 percent of sage stands are closed.

3) Acres dominated by aspen in forested areas. The intent

of the treatments is to increase the aspen acreage at the

expense of montane fir and spruce/fir. 

4) Acres of low-risk stands in dry forest areas. These are

the poletimber and larger stands historically maintained

by frequent surface fires. Most dry forest stands are

presently at high-risk for a stand replacement wildfire

due to missed fire cycles and the associated fuel buildup. 

Figure 2 shows the projected acreages dominated by

pinyon-juniper for the three scenarios. In the full-management scenario pinyon-juniper dominated acres decrease 41

percent from 115 to 68 thousand acres over the next 50

years. With no treatment this acreage would increase 35

percent to 155 thousand acres over this same period. At

current treatment levels pinyon-juniper acres are relatively

constant.

The ratio of open to closed shrub area is graphed in 

figure 3. This ratio increases substantially with the fullmanagement scenario and slightly when no treatments are

applied. Changes in the no-treatment option are in response

to projected wildfire acreage. If wildfires continue at the

160

Figure 2—Scenario comparison: Pinyon-Juniper dominant acres
rate of the past 10 years, there may be an increase in open

shrub areas, regardless of whether management treatments

are applied.

Acres dominated by aspen over the 50-year period are

graphed in figure 4. Dominant aspen acreage increases

from 9 to 24 thousand acres with full management and

from 9 to 17 thousand acres at current treatment levels. 

The increase in dominant aspen acreage associated with

the no-treatment scenario (from 9 to 11 thousand acres)

results from projecting the wildfire rates from the past

decade for the next 50 years. Presently there are approximately 17 thousand acres of spruce/fir forests that contain

aspen. Aspen is presently dominant on only 3 percent of

this area (500 acres). Since wildfires in sprucefir/aspen

stands normally result in an aspen dominated stand, the

acres of dominant aspen acreage increases over the next 

50 years even at a burn rate of 0.33 percent per year. Prior

to European settlement, aspen was dominant on about 35

thousand acres in the Beaver River drainage (Campbell 

and others, 2003).

Aspen reproduces primarily by suckering from the parent root system. This suckering is normally in response to a

disturbance such as wildfire. When aspen are lost from the

landscape they will not reseed an area, as do conifers (Bartos

and Campbell 1998). Overtopping by conifers, and reduced

vigor caused by grazing, is gradually causing aspen to die

out of many stands. Figure 5 shows the projected acreage

in which aspen is present, but not necessarily dominant.

Because all at-risk stands containing aspen cannot be immediately treated, the acres containing aspen declines with all

three scenarios.

Figure 6 compares the projected acres of dry forest that

are at a low risk for uncharacteristic stand-replacing fires.

These are the open ponderosa pine and mixed ponderosa

pine and Douglas fir stands. The acreage of low-risk stands

steadily increases with the full-management scenario but

ultimately decreases with the no-treatment and currentmanagement scenarios.

Management treatments in this study area are generally

designed to restore ecological conditions and reduce the

161

Figure 3—Scenario comparison: Ratio of open/closed shrub

Figure 4—Scenario comparison: Dominant aspen acres

Figure 5—Scenario comparison: Acres with aspen present

Figure 6—Scenario comparison: Dry Forest, low risk acres
intensity of wildfire. In range areas within the Beaver River

drainage these treatments result in a slight increase in

potential wildfires. This is due to the increased acreage of

annual and perennial grasslands, which are more prone to

fire. In forested areas management results in a small decline

in potential wildfire acres. The net result of these two factors

is that the total wildfire acreage is projected to be nearly the

same for the three scenarios.

DISCUSSION

VDDT is a non-spatial model intended mainly for broad

scale analysis. VDDT projects changes in vegetative conditions in response to succession, disturbances, and management treatments. In the Beaver River drainage key indicators

are used to estimate environmental conditions under different management scenarios, as projected by VDDT.

Major shifts in vegetative composition have occurred 

as a result of fire suppression and grazing. Historically

fires burned vegetative types within the Beaver River area

at frequencies ranging from every 5 to 80 years. In the past

decade only 0.57 percent of this area burned per year for an

average fire frequency of 175 years. In the previous two

decades (1973 – 1992) virtually no acres burned. Current

conditions in this drainage are reflective of conditions on

the Fishlake National Forest as a whole. On the Fishlake

National Forest spruce and fir acreage are estimated to have

increased 238 percent from historic levels, while aspen has

declined 259 percent. Pinyon-juniper acreage is believed to

have increased 357 percent, while sage/grass/forb acreage

has declined 295 percent (Campbell and others, 2003).

These changes result in reduced streamflows caused by

increased coverage of conifers including pinyon-juniper,

reduced habitat for sage and aspen dependent wildlife species,

and increased risk of uncharacteristically severe wildfires

in dry forest areas. 

This analysis suggests that if management treatments

are maintained at current levels, vegetative conditions will

most likely stabilize or decline further over the next fifty

years. Active management, with treatment levels designed

to mimic natural processes, would be necessary to restore

this area and reverse current environmental trends. 

The VDDT software, manual, and tutorial exercises can

be downloaded from the website http://www.essa.com.

LITERATURE CITED

Arbaugh, M.J.; Schilling, S.J.; Merzenich, J.; Van 

Wagtendonk, J.W. 1999. A Test of the Strategic Fuels

Management Model VDDT using Historical Data from

Yosemite National Park. In Proceedings of The Joint

Fire Sciences Conference and Workshop. Boise, ID

June 15-17, 1999. Vol II 85-89. Univ of Idaho press.

Bartos, D.L.; Campbell, R.B. 1998. Decline of Quaking 

Aspen in the interior West – Examples from Utah.

Rangelands. Vol 20, No. 1, February 1998. San

Antonio, Texas: Society for Range Management.

Beukema, S.J.; Kurz, W.A.; Pinkham, C.B.; Milosheva, 

K.; Frid, L. 2003. Vegetation Dynamics Development

Tool users’ guide, Version 4.4, ESSA Technologies Ltd.,

Vancouver, BC. 234 pp.

Bradley, A.F., Noste, N.; Fischer, W.C. 1992. Fire 

Ecology of Forests and Woodlands in Utah. Gen. Tech.

Rep. INT-287. Ogden, UT: U. S. Department of

Agriculture, Forest Service, Intermountain Research

Station. 128 pp. 

Campbell, R.B., Chappell, L., Greenhalgh, K. 2003. 

Personal communication, data from unpublished

Fishlake National Forest Wildland Fire Use Plan. 

Table 16: Analysis by Vegetation Cover Type for 

the Combination of all Subsections for the Fishlake

National Forest. Richfield, UT 84701.

Long, D. and Merzenich, J. (in prep). Chapter 19: 

Development and application of the LANDSUM and

VDDT models for the Southern Utah Fuel Management Demonstration Project. In: Long, Donald; Ryan, 

Kevin, editors. The Southern Utah Fuel Management

Demonstration Project: Development andApplication of

Models. General Technical Report RMRS-GTR vol 4.

FortCollins, CO: U.S. Department of Agriculture, Forest

Service, Rocky Mountain Research Station.

Merzenich, J.; Kurz, W.A.; Beukema, S.J.; Arbaugh, 

M.; Schilling, S. 1999. Long Range Modeling of

Stochastic Disturbances and management treatments

using VDDT and TELSA. In Proceedings of the SAF

1999 National Convention, 349-355, Bethesda, MD:

Society of American Foresters. 

162
O’Brien, R.A. 1999. Comprehensive inventory of Utah’s 

forest resources, 1993. Resource Bulletin RMRS-RB-1.

Ogden, UT: U.S. Department of Agriculture, Forest

Service, Rocky Mountain Research Station.105 p.

Payson, T.E., Ansley R.J., Brown, J.K., and others. 

2000. Fire in western shrubland, woodland, and grassland ecosystems. In: Brown, J.K., Smith, J.K., editors.

Wildland fire in ecosystems: Effects of fire on flora.

Gen. Tech. Rep. RMRS-GTR-42-volume 2. Ogden, UT:

U.S. Department of Agriculture, Forest Service, Rocky

Mountain Research Station, pp 121-159. 

Weise, D.R., R. Kimberlin, M. Arbaugh, J. Chew, J. 

Merzenich, J.W. Van Wagtendonk, and M. Wiitala

1999. A Risk-Based Comparison of Potential Fuel

Treatment Trade-Off Models. Proceedings: The Joint

Fire Sciences Conference and Workshop. Boise, ID

June 15-17, 1999. Vol II pp 96-102. Univ of Idaho 

press.

163
This page is intentionally left blank.
FIRE SUPPRESSION, 

FIRE PLANNING, AND

FUEL MANAGEMENT
This page is intentionally left blank.
In: Bevers, Michael; Barrett, Tara M., comps. 2005. Systems Analysis in Forest Resources: Proceedings of the 2003 Symposium. General Technical Report

PNW-GTR-656. Portland, OR: U.S. Department of Agriculture, Forest Service, Pacific Northwest Research Station.

1 M. Acuña is Graduate Research Assistant, Department of Forest Engineering, Oregon State University, Corvallis, OR 97331. 2 C. Palma is Professor 

of Forest Management, Departmento de Ciencias Forestales, Pontificia Universidad Catolica de Chile, Santiago, Chile. 3 A. Weintraub is Professor of

Logistics and Operations Management, Department of Industrial Engineering, University of Chile, Santiago, Chile. 4 D. Martell is Professor of Fire

Management, Faculty of Forestry, University of Toronto, Toronto, Canada. 5 W. Cui is a Postdoctoral Assistant, Faculty of Forestry, University of Toronto,

Toronto, Canada.

Systems Analysis in Forest Resources: Proceedings of the 2003 Symposium

INTEGRATED TIMBER HARVEST

AND FIRE MANAGEMENT PLANNING

Mauricio A. Acuña1, Cristian D. Palma2, Andrés Weintraub3, David L. Martell4, and Wenbin Cui5

INTRODUCTION

Fire, a natural component of most forest ecosystems,

has a significant impact on peoples’ safety, public and 

private property, and economic activities based on forestry

resources, such as wood production and recreation (Martell

and Boychuk, 1997). Likewise, fire plays a significant role

in innumerable functions of the ecosystem, such as biodiversity, the creation of habitat for animal life, etc. 

Different studies have examined the effect of fire on wood

supply (Reed and Errico, 1986; Van Wagner, 1979, 1983;

Martell, 1994), applying linear programming approaches

and stochastic losses. These studies all indicate that very

small burn fractions translate into substantial reductions to

long-term wood supply. Moreover, all these approximations

have considered fire losses as variables exogenous to the

problem of forest planning, without considering that forest

ABSTRACT

Harvest planners often consider potential fire losses and timber production plans can influence fire management, but

most timber harvest planning and fire management planning activities are carried out largely independently of each other.

But road construction, timber harvesting and silvicultural activities can influence the flammability of a forest, so the procedure applied in this work considers an integrated fire and forest management planning methodology that accounts for and

exploits such interactions using a fire spread model, a network model that identifies crucial stands that can influence the

spread of fires across a landscape, and a spatially explicit timber harvest scheduling model. The model and how it is being

applied to a forest management unit in the boreal forest region of Canada are described.

KEYWORDS: Mixed integer programming, heuristics, forest planning.

167

management decisions can have positive effects in terms of

reducing fire spread. Thus, for example, Hirsch and Pengelly

(1997) state that including certain forest management criteria, such as harvest planning and the design of cutting units,

could reduce the probability of fires by considerable magnitudes. In other words, the strategic fragmentation of the

forest that results from harvesting could have a firebreak

effect, which would translate into a reduction of the burn

area.

The forest management planning systems use the outputs

from the fire management system (annual burn fraction) as

an exogenous parameter. Forest management decisions,

however, particularly the areas and timing of harvesting

directly affect fuel continuity and therefore fire spread.

This integration requires the design of decision-making

support systems that can be used jointly by persons associated with both forest and fire management to develop,
evaluate and implement strategies and policies that contribute to sustainable forest management.

Our objective was to apply an integrated methodology that

includes the use of forest planning spatial models, heuristic

procedures and fire simulation models to evaluate the effect

of fire spread on harvesting in a real case. 

METHODOLOGY

Area of Study

The area of study selected to evaluate the integration 

of forest and fire management through spatial models is

located in the centre-west region of the province of Alberta,

Canada, and covers an area of more than 20,790 hectares

(rectangular area of 16.5 km x 12.6 km) managed by the

company Millar Western Industries. Of this total area, about

17,270 hectares correspond to forested areas used in this

evaluation, with the remainder occupied by bodies of

water, roads and towns. 

General Approach

Given the close interaction between forest and fire 

management decisions, their integration is proposed using

an iterative approach to planning as figure 1 indicates. 

The proposed methodology includes the interaction of

three systems: (1) a fire spread model, (2) a heuristic procedure for assigning a value to interrupting fire spread, and

(3) a spatial model for forest planning. Its main functions

are described below and a detailed description of (2) and

(3) is presented further along.

(1) Fire Spread Model: This model permits the simulation

of burn occurrences over time, based on spatial information for the study area (type of vegetation, topographical conditions, etc.) and weather information

such as wind direction and speed. The information 

provided by this model and relevant to our study is 

as follows:

- Fire spread times from each point (cell within the 

regular grid) to its neighboring points and the probability of a fire starting in each, to be used in (2).

- Statistics on burnt species to be able to estimate 

annual burn fraction, to be used in (3).

(2) Heuristic for assigning a value to interrupting fire spread:

this heuristic procedure assigns a relative “value” to

harvesting an area of the forest according to its positive

impact on interrupting the fire’s advance. It is based on

calculating the shortest routes for fire spread between

two points, with the result used in spatial planning of

forest harvesting. For implementation, the study area was

divided into 231,000 cells, each measuring 30 x 30 m.

(3) Spatial forest planning model: This corresponds to a

mixed integer programming model that maximizes the

net present value (NPV) and whose decision variables

include when and where to harvest the forest. The model

168

Figure 1—Iterative Planning Approach
includes losses of wood that result from fires, expressed

as the annual burn fraction, and considers the positive

impact that harvesting can have on interrupting fire

spread when it is associated with management decisions

to achieve minimum “protection” levels. 

This system proceeds as follows. Initial information

about the state of the forest feeds the fire module. This

generates (a) the information necessary (fire spread times,

probable fire starts, etc.) for calculating a “protection value”

using a heuristic procedure and (b) information on the fraction of the forest that burns given current conditions for

determining management decisions within the spatial model.

Meanwhile, protection values generated by heuristic procedures are included in the spatial model. Because the spatial

harvesting decisions obtained modify the original conditions

for fire spread and, as a result, its behavior and the fraction

of forest burned, spatial model decisions are again fed into

the fire module, generating new information for calculating

protection values and for the spatial model. The process is

stopped upon reaching a specific convergence criterion.

Procedure for Assigning a Value to Interrupting 

Fire Spread 

This heuristic procedure assigns a value to each cell 

for interrupting the fire spread. As a fire starting in a cell

should reach a destination cell along the route of least

resistance to the fire it can be represented by the different

shortest routes found for each pair of cells (Finney, 2002).

The importance of each route, R(i,j), will vary according to the ignition probability in the cell of origin I(i), the

value to be protected in the destination cell V(j), and the

fire spread time between both cells i and j T(i,j), and it is

calculated as R(i,j) = I(i) * V(j) / T(i,j). The higher the

probability of a fire occurring at a point or the higher the

value to be protected in the destination cell (which will

burn if the fire reaches it), the more important it is to block

this route. Similarly, if the fire takes too much time to reach

the destination cell, blocking its spread by harvesting a cell

will be less important, since it is more likely the fire will

be suppressed through active measures. In our case, the

value to protect in the destination cell corresponds to the

volume of wood available, but it could also consider the

presence of towns or critical infrastructure. The value of

the different routes, R(i,j), is added to cells that belong to

the corresponding route, so the protection value assigned to

each cell of forest will vary according to the quantity and

importance of the fires that it can interrupt. 

Finally, the procedure permits the creation of an index

or ranking that represents the cells’ belonging to different

routes for possible fires. Thus, regions prone to fire spread

will form part of many fire routes and their protection values will be higher. 

Spatial Forest Planning Model 

The model is for a mixed integer programming problem,

implemented by using the modeling language GAMS/IDE

v.19.6® and the solver CPLEX v.7.1®. The model establishes a spatial harvesting schedule that integrates fire protection. To do so, the model uses as its basis Model III by

Reed and Errico (1986) and the variation proposed by

Martell (1994), in order to maximize the NPV over an 80-

year planning horizon, divided into eight 10-year periods.

Model assumptions included:

- Completely accessed forest, so no road-building 

decisions were included. 

- Harvest flows (volumes of wood) were accepted 

varying ± 30% between periods.

- Harvested and burned areas regenerate naturally at no

cost. The model does not include reforesting decisions,

which means that once harvested the forest will always

reforest with the same species. 

- Transportation and harvesting (logging) costs per cutting unit were included, along with income per m3 of

wood harvested.

- There is an inventory of every age class when the 

planning horizon ends.

- Harvesting occurs at the mid-point of the period.

- If the decision is made to cut a stand during the period,

it is completely harvested.

- No adjacency restrictions are included.

- It is assumed that the fire consumes some fixed proportion of the area not harvested in each age class and that

this proportion does not vary over time. As indicated

previously, this is referred to as the burn fraction.

- Fire protection values are included for each cutting unit.

To do so, protection values for all the cells belonging to

a cutting unit are added together to reach its final protection value.

As output, the spatial model determines a harvesting

schedule for all cutting units over time and the supply of

wood arising from the allowable annual cut.

Information the Model Requires—As input, this

model receives the initial state of the forest, which consists

mainly of cutting units, species, the commercial volume

per species and age classes, all of which form part of Millar

Western Industry’s inventory information for the study area.

Cutting unit: To define the cutting units we used the geographic information system ArcInfo®, as our objective was

169
to generate the largest number of units of a single species

and age class, or a single species and different age classes.

Moreover, unit area was restricted to the 15-45 hectare

range, for a total of 628 cutting units averaging 27.5

hectares per unit. 

Wood volumes, ages and species: Using the above definition of cutting units, a species and age class representing

each cutting unit in the initial forest can be associated to

run the spatial harvesting model. Likewise, to calculate

commercial wood volumes, yield tables per species and

age class were used, as taken from inventory reports by the

company (Millar Western Forest Products Ltd., 2000).

Model Representation Using Networks—Our modeling focus at the forest level can be represented as networks,

one for each cutting unit (figure 2), using the outline presented in Martell and Boychuk (1996). 

Figure 2 depicts the area outcome over time. To do so,

the area in each age class, that is, a strata-based approach,

characterizes the forest at the beginning of each period.

The nodes represent the existing surface area of each age

class at the beginning of each period. The arcs represent

the surface flow among nodes. Thus, if the decision is

taken to harvest a cutting unit in a given period (figure 2a)

the entire area harvested of each age class is brought

together in a harvesting node, which will belong to the

regeneration age class (age class 1) in the next period. 

Similarly, if a cutting unit is not harvested (figure 2b),

the whole area in each age class is transferred to the next

age class for the following period. Logically, the area of

some nodes and arcs in each period of the network may 

be zero. 

For both situations within the same network a burn node

is used, which brings together the whole area for each age

class and each period that each cutting unit is affected by

fire. As indicated previously, the burn is a fixed fraction

that does not change over time or with each age class, so

independently of whether a cutting unit is harvested or not,

the area affected by fire is always discounted. This surface

area, as with harvesting, moves to the next regeneration

age class in the following period. With this formulation the

spatial detail is maintained during the periods as the area

harvested and burned as well as the area transferred to the

next period are tracked and known. 

RESULTS AND DISCUSSION

Information used by the model and dimensions of the

problem to be optimized

For the base problem a total of 464 cutting units were

considered. The burn fraction in each period was 5 percent,

the interest rate was set to 5 percent and the maximum

variation of volume between periods was set to 30 percent.

170

Figure 2—The network modeling approach. Cases in which the cutting unit is harvested (a) and not harvested (b).

a) b)
As for the problem to be optimized, this one had the

following features:

- Total number of variables: 345,260

- Number of binary variables: 3,712

- Number of equations: 349,443

- Number of iterations: 53,557

- Resolution time: 801 seconds

Analysis

A number of scenarios were evaluated in order to investigate the variation of timber volume in each period and the

economic effects of different strategies of fire protection.

Effect of different burn fractions on the harvest

schedule—The model was run varying the annual burn

fraction starting from a 0 percent level to a burn fraction 

of 10 percent. Five different scenarios were evaluated

whose burn fractions per period were 0, 1, 3, 5, and 10

percent, respectively.

The results obtained from these scenarios, their effects

on the objective function as well as the area and the total

volume harvested are shown in Table 1.

With an annual burn fraction of 10 percent the volume

and the objective function decrease about 24 percent in a

period of 80 years compared to the situation with no fire.

This situation can be more negative if the results are projected over a longer-term horizon, in other words, 200

years or more, where the losses can be dramatically higher

than the ones presented here (Martell, 1994). 

171

Table 1—Impact of different levels of burn fraction on volume harvested and NPV

Total area Total volume Objective

Scenario Burn fraction burned (ha) harvested (MM m3) function (MM $)

0 0% 0.0 1.758 21.438

1 1% 129.6 1.714 20.907

2 3% 388.9 1.635 19.899

3 5% 648.2 1.554 19.068

4 10% 1296.4 1.341 16.312

Table 2—Effect of fire protection on the total timber volume harvested

Burn fraction Burn fraction Volume before Volume after

before protection after protection protection protection Difference 

Protection (%) (%) (%) (MM m3) (MM m3) (%)

10 5 4.85 1.554 1.554 0

10 5 4 1.554 1.590 2.2

10 5 3 1.554 1.625 4.6

10 5 2 1.554 1.668 7.3

Table 3—Effect of fire protection on the Net Present Value (NPV)

Burn fraction Burn fraction NPV before NPV after

before protection after protection protection protection Difference 

Protection (%) (%) (%) (MM $) (MM $) (%)

10 5 4.85 19.068 19.068 0

10 5 4 19.068 19.376 1.6

10 5 3 19.068 19.697 3.3

10 5 2 19.068 19.978 4.8
In this work a burn fraction of 5 percent was used to run

the models. With this relatively small amount of fire on the

landscape, it is possible to observe losses of about 11 percent in volume and in NPV, corresponding to 204 million

cubic meters and 2.37 million dollars respectively. Logically,

these losses are considerably higher if the burn fraction is

increased (e.g. 30 or 40 percent). As the losses associated

to the presence of fire appear to increase exponentially,

they can become an important risk for the company’s wood

supply.

Effect of different fire protection levels on the 

harvest schedule and their impact on the annual burn

fractions—As described in the methodology, a number of

different fire protection levels were used for running the

model and evaluating its effects. 

Being an ongoing research the results are preliminary

and some of them are generated under assumptions (tables

2 and 3). Firstly the model was run with a 5 percent burn

fraction which represented the base level and whose value

was of 3.302E+7 units. Then we used two different protection levels, 10 and 20 percent, which were included into

the model along with the benefit of a reduced area burn

achieved with this higher protection. In this paper we only

present the results for a 10 percent of protection.

We analyzed the case of timber volume and NPV. For

the case with 10 percent of protection it is necessary to

reduce the actual burn fraction 0.15 percent, to a level of

4.85 percent in order to have no additional benefits in volume or in NPV. With a higher fire protection the model

becomes more constrained which means a decrease both in

the objective function and the volume harvested. However,

at the same time this negative effect in the objective function is neutralized with a reduced burn fraction as the

protection increases. 

If it is possible to reduce the burn fraction only from 5

percent to 4 percent, the total volume harvested increases 

a little more than 35,000 cubic meters, or a 2.2 percent.

Likewise, reducing the burn fraction 1 percent with a protection of 10 percent increases the NPV 1.6 percent, or

$307,991. 

Conversely, in an optimistic scenario with a protection

level of 10 percent over the base level and a reduction of

the burn fraction to 2 percent the volume increases in 7.3

percent or more than 114,166 cubic meters compared to the

scenario with no protection. The positive effect in NPV is

slightly lower than the benefit in volume. NPV increases

$910,000 for the entire horizon or 4.8%. As we used logging costs generated by a uniform distribution and a fixed

income for each cubic meter harvested, different results

would be obtained by using real information. 

CONCLUSIONS

From these results further research might be focused on

studying the effect of having higher initial burn fractions in

the landscape with no protection and the effects of greater

reductions in burn fraction as more protection is applied. It

is clear that the benefits are higher in landscapes strongly

affected by fire and when the spatial harvest planning generated by this methodology is able to considerably reduce

the annual burn fraction. 

It is also necessary to study the effects of different forest

structures and age class distributions and their impact on

the timber volume to be harvested in each period. The forest information and data used for carrying out this study

present much heterogeneity characterized by a considerable

number of different species and ages. In a commercially

managed forest the positive effect of applying this methodology both in volume and NPV can be higher than the

results presented here. Therefore, this methodology should

be tested in different kinds of forests with different structures and management in order to have more realistic

results and conclusions.

Finally it is important to continue studying and evaluating the methodology presented in this paper, in particular

the procedure for calculating the protection values as a

measure to stop fire propagation. We have just considered

the shortest paths between pairs of cells in order to know

the different routes for the fires. It might be useful to evaluate the effect of other paths, perhaps the K-shortest paths

that fires can take as alternatives routes. As this can be

computationally challenging and may not provide better

answers, other solution schemes such as using simulation

in a pro-babilistic way, for example Monte Carlo simulation, might be studied in a further research. 

LITERATURE CITED

Finney, M.A. 2002. “Fire Growth Using Minimum 

Travel Time Methods.” Can. J. For. Res. 32(8): 

1420-1424.

Hirsch, K. and I. Pengelly. 1997. “Forest Fuels 

Management in Theory and Practice.” Pp. 112-116 in

Proc. Symp. on Stand Density Management: Planning

and Implementation, C.R. Bamsey (ed.). Edmonton,

Alberta, Canada. 

172
Martell, D.L. 1994. “The Impact of Fire on Timber 

Supply in Ontario.” For. Chron. 70(2): 164-173. 

Martell, D.L. and D. Boychuk. 1996. “A Multistage 

Stochastic Programming Model For Level Timber

Supply Under Risk of Fire.” For. Sci. 42(1): 10-26. 

Martell, D.L. and D. Boychuk. 1997. “Level of Fire 

Protection For Sustainable Forestry In Ontario: Final

Report.” Report prepared for the Canada-Ontario Northern

Ontario Development Agreement, Sustainable Forestry

Development/Decision Support Project: Developing

Analytical Procedures for Establishing the Level of

Protection for Forest Fire Management to Support

Sustainable Forestry in Ontario.

Millar Western Forest Products Ltd. 2000. Detailed 

Forest Management Plan 1997-2000. Vol. 1, Ch. 1.

Edmonton, Alberta, Canada. 

Reed, W.J. and D. Errico. 1986. “Optimal Harvest 

Scheduling at the Forest Level in the Presence of the

Risk of Fire.” Can. J. For. Res. 16: 266-278.

Van Wagner, C.E. 1979). “The Economic Impact of 

Individual Fires on the Whole Forest.” For. Chron.

55(2): 47-50. 

Van Wagner, C.E. 1983. “Simulating the Effect of Forest 

Fires on Long Term Timber Supply.” Can. J. For. Res.

13: 451-457.

173
This page is intentionally left blank.
In: Bevers, Michael; Barrett, Tara M., comps. 2005. Systems Analysis in Forest Resources: Proceedings of the 2003 Symposium. General Technical Report

PNW-GTR-656. Portland, OR: U.S. Department of Agriculture, Forest Service, Pacific Northwest Research Station.

1 M.R. Wiitala is an Operations Research Analyst and 2 A.E. Wilson is a Resource Information Specialist at the Pacific Southwest Research Station, USDA

Forest Service, 620 SW Main Street, Suite 400, Portland Oregon, 97205; e-mail: mrwiitala@fs.fed.us.

Systems Analysis in Forest Resources: Proceedings of the 2003 Symposium

A SYSTEMS ANALYSIS MODEL FOR WILDLAND 

FIRE PREPAREDNESS PLANNING 

Marc R. Wiitala1 and Andrew E. Wilson 2

INTRODUCTION

The initial attack system remains a critical component

in the protection of the nation’s resources at risk to wildland fire. Planning the appropriate amount and use of firefighting personnel and equipment to effectively deal with

the fire fighting workload presented by impending fire 

seasons remains problematic for federal land management

agencies (GAO 2002).

Operations research technologies are available to improve

the quality of analytical tools for wildland fire planning.

The method of system simulation is particularly suitable

for application to the fire planning problem. Industry has

often used simulation to model complex systems such as

road transportation systems, emergency services delivery,

mining operations, and supply systems. In Canada, the

province of Ontario is successfully using system simulation

for fire preparedness planning (Martell and others 1994;

McAlpine and Hirsch 1999). System simulation offers a

175

great opportunity to capture the salient operating characteristics of the initial attack resource delivery system. An

appropriately designed computer simulation model would

provide fire managers a powerful analytical tool to evaluate

the relative performance of initial attack program alternatives and resource deployment policies. Investments in other

wildfire protection programs, like prevention, detection,

and fuels treatment, would also lend themselves to examination by the model.

The need for a modern wildland fire planning and policy

analysis tool is the target of a research and development

project funded by the National Fire Plan (USDA and USDI

2001). The project has designed and built a new initial

attack simulation model based on currently available simulation technologies. The model is currently undergoing testing and validation on the Umatilla National Forest in eastern

Oregon and Washington. This paper describes the anatomy

of this simulation model—the Wildfire Initial Response

Assessment System (WIRAS). Before beginning, we briefly

ABSTRACT

Planning an effective and efficient initial attack wildfire suppression organization in preparation for an uncertain fire season

is a challenging task for fire planners. This paper describes features of WIRAS, a discrete-event stochastic simulation model

that can assist fire planners in wildfire preparedness planning and policy analysis. WIRAS simulates, according to user-defined

rules, all phases of initial attack resource management and movement to meet the service needs in response to the ebb and flow

of fire workload over the landscape and through time. In an operational setting WIRAS is used to test alternative preparedness

programs against annual historical fire occurrence datasets to gauge performance with respect to spatially and temporally variable fire seasons. Planning capability at the national and regional scales is complemented by a unique ability of WIRAS to

evaluate local preparedness programs within a larger geographic context to account for the impact of external competition for

shared suppression resources, like airtankers, helitack, and smokejumpers.
look at the nature of fire planning problem and the initial

attack environment.

PLANNING PROBLEM

In preparation for fire seasons, wildland fire managers

choose the size, composition, and location of the fire suppression organization that best meets land management and

protection objectives, usually with respect to alternative

levels of fiscal constraint. Managers make this decision in

the face of uncertainty, where fire seasons in many locations of North America vary greatly in intensity over time

and across the landscape.

Choosing an appropriate presuppression organization

more specifically involves selecting, for example, the number, size, location, and term of hand and engine crews to

hire and long-term investments in capital equipment and

facilities like fire engines and helicopter bases. Fire preparedness planning also develops dispatch policies for using

not only locally controlled resources but also for those shared

among jurisdictions, for example, airtankers and smokejumpers. 

Initial Attack Environment

The initial attack setting has been aptly described by

Martell and others (1984) and Wiitala (1998). We summarize the account given by the latter, but with reference to

initial attack activity on federal lands in a region containing

northern California, the Pacific Northwest, Idaho, and

western Montana. 

The initial attack process in this broad area consists of

the collective action of many administrative units deploying wildland fire suppression forces in response to a seasonal flow of fires distributed temporally and geographically.

The number of fires in this region can exceed 10,000 per

summer season. Several hundred fires can occur on a given

day and approach 1,000 on several consecutive days 

(figure 1). 

Because thunderstorms are responsible for these peak

periods of fire activity, fires tend to cluster spatially and

temporally. This phenomenon creates localized heavy fire

workloads that can tax the capabilities of local initial attack

organizations and result in larger project fires. In a large

geographic area, multiple thunderstorm cells often give rise

to spatially dispersed fire clusters, creating competition for

aerial initial attack resources shared over larger geographic

scales.

176

Figure 1—Comparing the temporal distribution of single and consecutive three-day fire starts in the northwestern United States in 1994.
As a general characterization of the local initial attack

process, when fire activity commences, the resource dispatching process begins. Rules for resource dispatching

vary between administrative units. The type and effort of

response vary with uncertain conditions of the fire environment, land management objectives, and available resources.

In general, dispatchers favor use of local ground resources

with the shortest response times, provided the response

times are reasonable for the conditions of the fire and protection objectives. Otherwise, dispatchers will deploy aerially delivered fire fighters, if available. Dispatchers may

also request the support of airtankers to deliver retardant 

or water to the fire. Or, after delivery by helicopter, a crew

may request a helicopter to support fireline construction

with water drops.

Initial attack resources work a fire until it is contained

or becomes unmanageable. In the latter instance, a large

fire suppression organization assumes management of the

fire. In either case, how quickly initial attack resources

demobilize and return to a ready status depends on a variety of conditions including the time required for fire mopup, the level of fire activity, time of day, and accessibility

to transportation.

MODEL

WIRAS uses a process-oriented approach (Law and

Kelton 1991) in its design of a discrete-event, stochastic

computer simulation model representing the initial attack

system. An enhanced variant of the discrete-event simulation

language General Purpose Simulation System (Hendriksen

and Crain 1989) originally developed by Gordon (1962), is

the language used to build the simulation model. WIRAS

shares many attributes and capabilities with the clock-driven,

next-event simulation models of Fried and Gilless (1988),

Martell et al. (1994), and Wiitala (1998), which address how

the frequency and magnitude of multiple fire events and

spatial interactions of initial attack affect program performance. 

WIRAS distinguishes itself from its predecessors by

featuring the ability to plan both local and national fire preparedness programs. It also pioneers the ability to conduct

local planning within the context of the national fire environment and fire fighting activity. This embedded planning

approach allows fire planners to design a local program

while taking into account their success in competing for

nationally shared aerial resources. To describe the structure

of WIRAS, we present its major components.

Suppression Resources

WIRAS contains distinct subsystems to simulate the

activities of various types of aerial and ground resources in

the initial attack system. Each resource subsystem features

a unique design to capture disparate operational features.

Airtanker, Smokejumper, and Helitack resources—

The aerial subsystems deliver both water and retardant and

crews to fires. The airtanker subsystem consists of single

and multiple engine airtankers delivering retardant from a

system of geographically fixed airtanker bases. Additionally,

the single engine aircraft can operate from a mobile base.

Not all bases can accommodate all airtanker types. For this

reason, each airtanker and retardant base carries attributes

within the model to permit imposing base use or capacity

restrictions on some airtankers. Attributes affecting airtanker operations within the model include airspeed, retardant

carrying capacity, loaded weight, refueling and reloading

times, and a single engine airtanker identifier.

Major time-consuming operations modeled at an airtanker base include landing, taxiing, taking-off, refueling,

and loading. The simulation model monitors daily airtanker

flight hours to determine when to refuel an airtanker as

well as when to limit use due to daily flight-hour restrictions. The number of loading bays, refueling trucks, and

parking spaces specified for bases affects their ability to

service airtankers and ultimately affects retardant delivery

system performance. 

The model’s helitack subsystem is responsible for the

delivery of crews, water, and retardant. After delivering a

crew to a fire, the model may augment crew fireline building activity with helicopter water bucket support. In the

model, bucket support continues until there is a need for the

helicopter to refuel, deliver a crew, or return to base before

nightfall. Helitack is modeled as a locally controlled resource

operating from a system of bases with varying numbers of

deliverable crew members. This subsystem permits sharing

helicopters between bases but not crews.

The fixed-wing aircraft deliver smokejumpers primarily

to fires that are remote or when the supply of local resources

is exhausted. Smokejumpers are highly mobile at the national

scale. Wiitala and Dammann (2003) provide a detailed description of the smokejumper subsystem as used in WIRAS.

In summary, smokejumpers and their delivery aircraft

attack fires from a system of bases dispersed across the

western United States. The simulation model moves smokejumper assets between bases in response to demands created

by the ebb and flow of fire activity across the national

landscape. 

177
Ground resources—Because ground resources are very

numerous, they not only make the greatest contribution to

initial attack fire suppression, they also present the greatest

obstacle to designing a computationally efficient simulation

model. Initial attack ground resources in WIRAS currently

consist of different sizes and types of hand crews, engines,

water tenders, and dozers. The simulation model manages

these resources in two different ways. The first way, used

at the broad scale, gains computation speed at the expense

of modeling accuracy. Here, ground resources lose their

individual identity and specific locations as the model

aggregates line-producing capabilities of disparate ground

resources into a pool of generic capability called crew person equivalents. Aggregation takes place over user defined

geographic areas. As fires occur during a simulation, the

model dispatches crew person equivalents from the resource

pool for the geographic area in which the fire occurs. As a

pool of resources declines, the remaining resources are

assumed to be more widely dispersed on an increasingly

coarse grid and, thus, to have a higher mean response time

to each new fire. To implement this approach, the model

makes a random draw from a probability model describing

distance to nearest grid point (resource location) to determine a distance for calculating the travel time of the first

arriving resource. Not having to track within the model the

many individual ground resources populating the national

landscape substantially increases the speed of simulation

runs.

While the above scheme for capturing the suppression

contribution of ground resources is sufficient when planning nationally shared resource programs like airtankers

and smokejumpers, planning a local organization requires

more detail for modeling local ground resources. In this

second method, WIRAS resolves this problem by accepting

a list of local ground resources containing specific resource

attributes like location, line production capability, and travel

speed. This detailed information allows the model to simulate more closely the deployment policies and travel dynamics

for the local resources with only minor loss of computation

efficiency. By simultaneously simulating fire and suppression resource activity both locally and nationally, WIRAS

allows local planners to take into account the competition

for nationally shared resources, like smokejumpers, while

planning individual components of their local preparedness programs. 

Fire Workload Generation

The number of ignitions, fire behavior, and resistance 

to fire control define fire workload. The combined effect of

178

Figure 2—An example of the spatial distribution of fires on a day.
temporal variation in fire occurrence (figure 1) and spatial

variation in ignition pattern (figure 2) is an important determinant of initial attack program performance both locally

and nationally. To preserve the impact of spatial and temporal autocorrelation on program performance, WIRAS

simulates program performance for an entire set of fires

from a season. Because fire workload varies significantly

between years, WIRAS addresses the implications of fire

workload variation on initial attack performance by testing

a preparedness program against a range of fire seasons.

Historical fire records currently provide WIRAS with the

necessary time stream of geographically located fire ignitions. 

Strength of Response

WIRAS mimics an administrative unit’s dispatch rules

for determining the strength and priority of initial response

to a fire. Initial response policies revolve around dispatch

blocks. Each dispatch block is a geographic area that reflects

a difference in one or more of the following: management

and protection objectives, fire behavior potential, and level

of fire activity (figure 3). 

As is true operationally, when the first resource arrives

on a fire, the simulation model reassesses the needed

staffing level based on actual fire conditions (figure 4).

Depending on assessed needs and the strength of the initial

response, the model may look for additional resources or

halt the response of some resources in transit to the fire. In

accordance with local practice, the model may also initiate

airtanker support.

Resource Dispatching and Allocation

WIRAS controls resource dispatching through a system

of Boolean variables that mimics the rules and priorities

governing an administrative unit’s dispatch policies (figure

5). The model maintains and updates several state variables

used by the Boolean variables to reflect the types of information typically available to dispatchers. This system of

Boolean variables is sufficiently flexible to permit planners

to set the conditions for the dispatching ground resources

based on attack times, fire behavior potential, and protection objectives.

179

Figure 3—Factors determining the preplanned initial response to a fire.

Figure 4—Factors affecting reassessment of fire staffing needs upon first

resource arrival.

Figure 5—Resource dispatch decision process for staffing fires.
During a simulation, if the conditions for a ground

resource response are not met, the model attempts to dispatch helitack crews. Failing this, WIRAS searches the

smokejumper system for available suppression resources.

The search for resources to dispatch to a fire continues

until meeting fire staffing needs as determined by the initial response or by reassessment upon first resource arrival.

If no initial attack resources are available, the model queues

the staffing request. When a resource becomes available for

dispatching, WIRAS selects from the service queue the

highest priority fire serviceable by the resource. The model

user can set the priority scheme based on protection objectives and fire behavior potential. However, WIRAS always

attempts to complete unfulfilled fire staffing needs before

staffing new fires.

Response Time Calculations

Accurate modeling of resource response times is critical

to the successful modeling of the initial attack system. In

WIRAS, components of response time include getaway,

reloading, refueling, and travel time. Travel time is usually

the largest component. For aircraft, distance between travel

points and airspeed combine to calculate travel time. During

periods of high fire activity, the response time of airtankers

can lengthen if they have to wait for a retardant loading

bay at congested airtanker bases.

For ground resources, WIRAS combines vehicle travel

along a road network and off-road travel for ground

resource travel time (Wilson and Wiitala, this proceedings).

The simulation model integrates a computationally efficient

least-cost-path algorithm to estimate ground resource travel

times. To improve simulation time, many of the ground

resource response times are computed prior to simulation

based on known origins and destinations. During a simulation, a ground resource in transit to a fire may be diverted

to another fire. At this point in the simulation, its location

is calculated and the least-cost-path algorithm invoked to

calculate response times to all fires to which the resource

might respond. 

Fire Growth and Containment

WIRAS implements a time-dependent accounting

process operating on a five-minute time step to calculate

fire perimeter growth, track resource arrivals, and compute

fireline to determine when a fire is contained or escapes

initial attack. The forward rate of fire growth is estimated

at 2:00 PM from historical data for the fuel and environmental conditions at the location and on the day of fire

occurrence. Fire growth is diurnally adjusted to account for

wide differences in growth rates over a 24-hour period. A

double elliptical growth algorithm calculates perimeter

expansion over time (Anderson 1983). Fire line construction

occurs on the flanks of the fire from the point of ignition.

WIRAS contains fires according to a model developed by

Quintilio and Anderson (1976) that accounts for the effect

of fireline production in moderating the growth of a fire’s

perimeter. 

Fires are considered contained when fireline reaches a

user specified percentage of fire perimeter. WIRAS escapes

a fire (removes it from initial attack status) when fire

perimeter exceeds constructed fireline by a user specified

amount or fire area reaches a defined threshold.

Demobilization and Return to Service

Assumptions regarding how quickly initial attack resources

become available for dispatching after fire containment or

escape can greatly influence program performance. WIRAS

encodes these assumptions into user changeable rules that

provide the opportunity to evaluate how alternative demobilization and return to service policies affect system performance. 

Under the existing default rules, when finished fighting

a fire, WIRAS returns initial attack resources to service on

the next day following fire containment and on the second

following day if the fire escapes initial attack. The model

makes two exceptions to this rule. First, if ground resources

in roaded areas contain a fire sufficiently early in the day,

the model may dispatch them to other fires awaiting service,

response time permitting. Second, those crews simulated as

walking deep in wilderness or unroaded areas will remain

with the fire a full burning period after containment. This

policy assures full fire control and minimizes the risk of a

subsequent response to a reigniting fire in a remote area of

high cost for aerial delivery or long response times for

ground-based crews. 

MODEL OUTPUTS

During a simulation, WIRAS collects a variety of operational statistics. These statistics help planners judge the

physical and economic performance of alternative approaches

to organizing a fire preparedness program for their administrative units. For example, information collected on

resource utilization allows planners to identify underutilized

resources and bases.

Planners can also combine a program’s cost with suppression costs and projected natural resource losses estimated from model outputs to evaluate relative economic

180
efficiency of alternative preparedness programs. Program

budgeting issues are also addressable. Other statistics include

acres burned, escaped fire numbers, fire dispatches, aircraft

flight hours, gallons of delivered retardant, and frequency

distributions for contained fire sizes and resource attack

times. 

CONCLUSION

The purpose of this paper was to briefly describe the

anatomy of WIRAS—a stochastic simulation model for fire

preparedness planning. Stochastic simulation is proving to

be a powerful operations research method for planning

highly complex initial attack service delivery systems.

Taking a system analysis perspective in the design of the

WIRAS, the resulting initial attack simulation model can

address a wide range of program and policy issues not previously addressable by other fire planning technologies.

Not only can fire planners use WIRAS to determine the

best program composition, they can also use it to improve

their resource deployment and dispatch policies—locally,

regionally, and nationally.

Opportunities for the application of WIRAS are now

becoming apparent. In one recent application at the national level, WIRAS identified effective alternatives to the use

of heavy airtankers in support of initial attack activities.

Geographic area command centers are interested in the

potential of WIRAS to examine and implement aerial

resource prepositioning strategies during a fire season in

response to fire severity forecasts. Also under consideration

is the use of WIRAS to quantify improved initial attack

program performance and consequent benefits of reduced

suppression costs resulting from hazardous fuels treatment.

Each new application of WIRAS generally necessitates creating a new variant of the core simulation model to address

the unique characteristics of the new problem. 

LITERATURE CITED

Anderson, H.A. 1983. Predicting wind-driven wild land 

fire size and shape. Res. Paper INT-305. Ogden, UT:

U.S. Department of Agriculture, Forest Service,

Intermountain Forest and Range Experiment Station. 

26 p.

Fried, J.S., and J.K. Gilless. 1988. Modification of an 

initial attack simulation model to include stochastic

components. In: Proceedings of the 1988 system 

symposium on forest resources. USDA For. Serv. 

Gen. Tech. Rep RM-161. pp. 235-240.

U.S. Government Accounting Office. 2002. Improved 

planning will help agencies better identify fire-fighting

preparedness needs: Report to Congress. GAO-02-158

Wildland fire management. Washington, D.C. 32 p.

Gordon, G. 1962. A General Purpose System Simulator. 

IBM Systems Journal 1(1). pp. 18-32.

Henriksen, J.O., and R.C. Crain. 1989. GPSS/H 

Reference Manual. Third edition. Wolverine Software

Corporation. Annandale, Virginia. 

Law, A.M., and W.D. Kelton. 1991. Simulation modeling 

and analysis. McGraw-Hill,Inc., New York. 759 p.

Martell, D.L., R.J. Drysdale, G.E.Doan, and D. Boychuk. 

1984. An evaluation of forest fire initial attack

resources. Interfaces 14(5), pp. 20-32.

Martell, D.L., D. Boychuk, J.I. , MacLellan, B.M. 

Sakowicz, and R. Saporta. 1994. Decision analysis 

of the level of forest fire protection in Ontario. In:

Proceedings of the sixth symposium on systems analysis in forest resources, 1994 September 6-9, Pacific

Grove. Society of American Foresters, Bethesda, MD.

pp. 138-149.

McAlpine, R.S. and K.G. Hirsch. 1999. An overview of 

LEOPARDS: The level of protection analysis system.

The Forestry Chronicle 75(4): 615-621.

U.S. Department of Agriculture, Forest Service. 1985. 

National fire management analysis system user’s guide

of the initial action assessment (FPL-IAA2.2).

Washington, D.C.

U.S. Departments of Agriculture and Interior. 2001. A

collaborative approach for reducing wildland fire risks

to communities and the environment: 10-year comprehensive strategy. Washington, D.C.

Quintilio, D., and A.W. Anderson. 1976. Simulation study 

of initial attack fire operations in the Whitecourt Forest,

Alberta. Canadian Forestry Service Information Report

NOR-X-166. 43 p.

Wiitala, M.R. 1998. Modeling the spatial and temporal 

dynamics of initial attack fire suppression using GPSS.

In: Proceedings of the III international conference on

forest fire research and 14th conference on fire and 

forest meteorology, 1998 November 16-20, LuscoCoimbra, Portugal. pp. 2389-2403.

181
Wiitala, M.R., and C.L. Dammann. 2003. Simulating a 

smokejumper service delivery system. Systems analysis

in forest resources: proceedings of the eighth symposium held September 27-30, 2000, Snowmass Village,

Colorado, U.S.A., G.J. Arthaud and T.M. Barrett (eds.),

Kluwer Academic Publishers, Dordrecht. pp. 71-79.

182
In: Bevers, Michael; Barrett, Tara M., comps. 2005. Systems Analysis in Forest Resources: Proceedings of the 2003 Symposium. General Technical Report

PNW-GTR-656. Portland, OR: U.S. Department of Agriculture, Forest Service, Pacific Northwest Research Station.

1 Andrew G. Kirsch is a Program Analyst with the National Park Service, NIFC Boise, ID 83705. 2 Douglas B. Rideout is a Professor of Fire Economics 

at Colorado State University, Fort Collins, CO 80523.

Systems Analysis in Forest Resources: Proceedings of the 2003 Symposium

OPTIMIZING INITIAL ATTACK EFFECTIVENESS 

BY USING PERFORMANCE MEASURES

Andrew G. Kirsch1 and Douglas B. Rideout2

INTRODUCTION

Estimating the preparedness organization for an upcoming wildfire season on public lands is a complex process.

Federal land management agencies expend considerable

effort to estimate the budget and to decide upon the resources

that will be needed to protect public lands from wildfire

each year. In preparedness planning, a fire manager’s decisions often include the type and number of firefighting

resources to have available. The difficulty of determining

the number of resources is often compounded by the occurrence of simultaneous ignitions and having to plan for the

possibility of receiving alternative appropriation levels. 

A suitable analysis of the tradeoffs of all of the possible

combinations of appropriation levels, resources, and uses,

requires a sound analytic system and a definitive measure

of performance. 

Many approaches have been suggested to assist decision

makers in analyzing parts of the preparedness problem.

ABSTRACT

Increased scrutiny of federally funded programs combined with changes in fire management has created a demand for 

a new fire program analysis model. There is now a need for a model that displays tradeoffs between initial attack effectiveness and alternative funding levels. The model is formulated as an integer linear program that operates in a performance

based, cost-effectiveness analysis (CEA) environment. Using the performance measure of weighted area protected (WAP),

the model employs a non-monetized approach to interagency fire planning. The model optimizes the initial attack deployment for a user-defined set of fires that a manager would like to be prepared for across alternative budget levels. The paper

shows how an integer programming formulation can use optimal deployment to address the annual fire program. It is also

shown, for the first time how simultaneous ignitions can be incorporated into an optimization approach. This compact and

robust model can provide the basis for a wider scale formulation with the potential to measure an organization’s performance and promote a higher level of accountability and efficiency in fire programs. 

183

Approaches employing optimization have their roots in

World War II and the birth of operations research. Much 

of the optimization literature extends Parks (1964) who

designed a deterministic model to minimize the cost of

suppression plus damages to find an optimal constant

workforce. See Parlar and Vickson (1982), Parlar (1983),

Aneja and Parlar (1984) for models that extend the Parks

formulation. Boychuk and Martell (1988) also evaluated

seasonal forest fire fighter requirements utilizing the least

cost plus loss framework by using Markov chains. Recently,

Donovan and Rideout (2003) used an integer linear programming model to optimize a firefighting resource allocation to a single fire using a cost plus net value change

framework. 

Despite the difference in optimization approaches, most

previous models have relied heavily on derivatives of the

cost plus loss theory first applied to wildfire in the early

20th century (Sparhawk 1925). The cost plus loss framework requires monetization of the damages caused by wildfires. As a fire advances across a landscape it can affect a
myriad of natural resources, each posing its own challenge

for accurate valuation. The valuation of fire damage continues to be an extremely difficult and costly task (Cleaves

1985; Pyne and others 1996). 

Non-monetized approaches have also been used in a

few fire management applications. For example, Kourtz

and O’Regan (1968) used a cost effectiveness analysis

(CEA) to assess fire detection systems, Nautiyal and Doan

(1974) used iso-dissatisfaction curves to assess trading

planned cut for wildfire protection expenditures, Omi and

others (1981) used a damage reduction index to optimize

fuel treatments, and Mees and Strauss (1992) used a utility

measure of the relative importance of holding different

constructed fireline segments to evaluate strategies for the

tactical deployment of resources to a single large fire. Our

research has been enriched by these previous studies as we

have developed an optimization model in a non-monetized

system. 

Here we provide the first non-monetized application 

of an optimization model applied to preparedness planning.

This approach shows how a strategic fire management

planning effort can utilize integer linear programming to

organize and optimize the initial attack response, in a CEA

framework. The output of the model gives the planner an

optimal list of resources to have available for an entire 

season across a range of possible budget appropriations.

The model accounts for simultaneous ignitions, which are

important planning elements that have not previously been

considered in an optimization model. 

DEVELOPING THE MODEL

Our approach demonstrates how a performance based

optimization model can inform strategic planning decisions.

For example, the model can provide information to aid

with the development of a menu of initial attack firefighting resources that would be available for the fire season

while considering both single and simultaneous fires. This

strategic scope is in contrast to tactical or real-time decision-making such as the deployment of specific resources

to a specific fire event. Because we optimize initial attack

performance, as opposed to extended attack or large fires,

only local resources affecting the planning unit’s budget

allocation are included in the model. 

We use a deterministic integer linear programming (ILP)

optimization approach to model the initial attack planning

problem for several reasons. First, strategic planning involves

decisions that are often integer or binary. These decisions

include the number, location, and type of resources to have

available for deployment to a set of fires. Second, an ILP

allows for the evaluation of thousands of decisions in a

compact and flexible formulation. Our optimization model

employs a performance-based CEA to display the tradeoffs

between funding levels and initial attack performance. Cost

effectiveness analysis helps decision makers allocate limited resources efficiently when performance measures are

used in lieu of monetized benefit estimates (Robinson and

others 1995; Osborne and Plastrik 2000). The results of a

CEA are often displayed by forming a frontier as in figure

1, where costs are compared with effectiveness. While all

points on the interior of the frontier are possible, they are

technically inferior to points that comprise the frontier. For

any interior point, there is at least one point on the frontier

that can be shown to be preferable. Optimization enables

the analyst to focus on solutions or points that define the

frontier.

To compare the effectiveness of different initial attack

organizations, we provide a performance measure defined

as the weighted area protected (WAP). A WAP is a geographical area (for example, acre or hectare) assigned a

proportional numerical weight representing the importance

of protecting that area relative to another area from damaging wildfires. The WAP combines both qualitative and

quantitative assessments of damage caused unwanted wildfires. Not all fires are created equal and the WAP captures

the importance of protecting our natural treasures in a way

that dollars values could not. Tradeoffs are clear to see and

are directly related to performance on the ground. Defining

WAP as a performance measure will allow the model to

weigh the possibility of using scarce resources to contain

more important fires while letting less important fires escape.

Using a frontier estimated from the collection of a program’s

analyses, the overall effectiveness of the initial attack

organizations can be displayed and analyzed as in figure 1. 

The model requires an input set of fires that reflect the

expected workload of a future fire season. Initially, the

model needs predetermined initial response analysis period.

This time can take on any value and each time step does

not have to be constant allowing for a flexible analysis.

Each fire is defined by an initial reporting size at time zero

and its total cumulative perimeter and area burned for each

time step in the analysis period. The fire’s perimeter is

directly related to cost through resource production rates

and the fire’s burned area is directly related to performance

through WAP. Other fire behavior characteristics such as

flame length and fire intensity can be reflected in the firefighting resources’ ability to construct fireline. This allows

managers to incorporate tactical firefighting standards,

such as a fire with flame lengths of four to eight feet can

184
be too intense for a direct attack with hand tools, but 

bulldozers, engines, and aerial drops can be effective

(BLM Standards Chapter 9 2003). We use the free burning

fire containment rule from previous deployment models

(for example, USDA Forest Service 1991; Donovan and

Rideout 2003) that states a fire is contained when the total

fireline produced by firefighting resources is greater than

the free burning fire perimeter. A fire is defined as having

escaped initial suppression efforts if it cannot be contained

in the initial attack time period because of either a lack of

funds or inadequate fireline production. 

The attributes of the firefighting resources constitute the

main set of inputs to the optimization routine. The model

requires a list of resources to choose from in order to maximize the WAP. This list includes all of the resources that

are potentially affected by a planning unit’s budget. Each

resource is defined by a total fireline production and by its

fixed and variable cost. Fireline production is input to the

model as cumulative values for every time step of each

fire. An advantage of this integer time step format is that

the production rate does not have to be constant or linear

and can reflect fatigue and other disruptions in the production such as water refills and refueling. Arrival and other

travel times are also reflected in these production values.

Resources produce zero chains of fireline during travel

periods. The model uses the production information along

with costs to solve for the optimal deployment. 

The costs of initial response resources are important

considerations in estimating the optimal deployment. While

previous approaches included cost information, costs did

not directly affect resource deployment. Our model relies

upon fixed and variable costs that are directly input to aid

with the management of optimal deployment (Donovan and

Rideout 2003). The fixed cost is modeled as a one-time

charge that is incurred only if the resource is deployed to

any fire during the season. Each resource’s variable cost is

modeled as an hourly cost that reflects its operating costs

on each fire, including maintenance, fuel, regular hourly

wages, overtime and hazard pay. 

Including simultaneous ignitions in the optimization

model adds depth and advancement to the analysis. To

model simultaneous ignitions we force certain resources 

to choose to fight a maximum of one of the simultaneous

ignitions. We assume that once a resource is deployed to 

a fire, it is not released when containment is achieved for

another deployment. 

The following is the mathematical representation of the

model:

Maximize WAP

I De

WAPo – ∑ ∑ (Wid * fid * Aid ) (1)

i=1 d=0

Subject to:

D

∑ xird ≤ ur i,r (2)

d=1

De

∑ fid = 1 i (3)

d=0

R D D

∑ ∑ (xird * Lird) ≥ ∑ fid * Pid i (4) r=1 d=1 d=0

D D

∑ d * fid ≥ ∑ d * xird i,r (5)

d=0 d=1

I R D R

∑∑∑ (xird * Hrd)+∑ur * Fr ≤TC (6)

i=1 r=1 d=1 r=1

D

∑ ∑ xird ≤ ur n,r ∈ Rs (7)

r∈Sn d=1

185

Figure 1—Generic cost-effectiveness frontier

A

A

A

A

A
Decision Variables

xird = Binary(0,1) = 1 if resource (r) is used for (d) time

periods on fire (i). 

= 0 if resource (r) is not used on fire

(i) for (d) time periods.

fid = Binary (0,1) = 1 if fire (i) burns for (d) time periods.

= 0 if fire (i) does not burn for (d)

time periods. 

ur = Binary(0,1) = 1 if resource (r) is deployed to any

fire for any duration.

= 0 if resource (r) is never used.

Parameters

I = set of all fires indexed by i.

R = set of all firefighting resources indexed by r.

D = resource deployment and contained fire duration

indexed by d.

De = escaped fire duration. De is defined as D+1.

Sn= (n)th set of simultaneous ignitions. 

Sn I.

Rs = set of firefighting resources that are restricted to fight

at most one of the fires that ignite simultaneously. 

Rs R.

Fr = fixed cost for resource (r).

Hrd = total hourly cost accrued for resource (r) for (d) time

periods.

Lird = total (cumulative) fireline produced by resource (r)

for (d) time periods during fire (i).

Wid = relative weight for the area burned by fire (i) after

(d) time periods.

Pid = total fire perimeter of the burned area for fire (i) after

(d) time periods. 

Aid = total area burned by fire (i) after (d) time periods.

Calculated from Pid.

TC = total cost of initial attack input to the model.

WAPo = total weighted area of the fire planning unit.

The objective function (1) maximizes effectiveness

defined as the weighted area protected for a given budget.

Equation (2) limits a resource to at most one deployment

per fire. For example, a firefighting resource cannot be

deployed for two time periods and four time periods to 

the same fire. The equalities in (3) force each fire to have

exactly one burn duration. Equations (4) are the containment constraints. For each contained fire, the total amount

of fireline produced by all resources must be greater than

or equal to the fire’s burn perimeter. Equations (5) force

each fire to burn at least as long as the longest duration 

of any resource deployed to that fire. For example, if on a

given fire there were two resources used, one for three time

periods and the other for six time periods, the fire would

burn for six time periods. Inequality (6) is the budget constraint. The total cost of all resources deployed to all fires,

both hourly and fixed, must be less than or equal to the

total cost denoted as TC. Equations (7) are used for fires

that are modeled as simultaneous ignitions. Fires in these

groups compete for firefighting resources. We assume that

a resource can only be deployed to one fire in each group

of simultaneous ignitions. 

CONCLUSION

Our formulation provides an approach to performance

based initial attack planning that incorporates the tenets of

CEA in an optimization model. The performance measure

of WAP shows how a non-monetized performance based

system could be applied while addressing key elements of

preparedness planning. While this model expanded previous work to optimize simultaneous ignitions, future efforts

could continue to expand the scope of analysis. For example, one could change the deterministic model to incorporate stochastic elements, such as a range of likely fire

occurrences. Other possible scope extensions could include

linear or piecewise linear approaches to reduce the solution

times. 

LITERATURE CITED

Aneja, Y.P., and M. Parlar. 1984. Optimal staffing of a 

forest fire fighting organization. Can. J. For. Res. 14:

589-594.

Boychuk, D., and D.L. Martell. 1988. A Markov chain 

model for evaluating seasonal forest fire fighter requirements. For. Sci. 34(3): 647-661.

Bureau of Land Management. 2003. Interagency 

Standards for Fire and Aviation Operations. BLM

Handbook 9213-1.

Cleaves, D.A. 1985. Net Value Changes in Fire Management:

Issues and Application Problems. Fire management: the

challenge of protection and use: proceedings of a symposium. Logan, Utah.

Donovan, G.H., and D.B. Rideout. 2003. An integer 

programming model to optimize resource allocation for

wildfire containment. For. Sci. 49(2): 331-335.

Kourtz, P.H., and W.G. O’Regan. 1968. A cost-effectiveness analysis of simulated forest fire detection systems.

Hilgardia 39(12): 341-366. 

186

⊃

⊃
Mees, R., and D. Strauss. 1992. Allocating resources to 

large wildland fires: a model with stochastic production

rates. For. Sci. 38(4): 842-853.

Nautiyal, J.C., and G.E. Doan. 1974. Economics of 

forest fire control: trading planned cut for protection

expenditure. Can. J. For. Res. 4: 82-90.

Omi, P.N.; J.L. Murphy; L.C. Wensel. 1981. A linear 

programming model for wildland fuel management

planning. For. Sci. 27(1): 81-94.

Osborne, D.E., and P. Plastrik. 2000. The reinventor’s 

fieldbook: Tools for transforming your government.

Jossey-Bass, San Francisco. 

Parks, G.M. 1964. Development and application of a 

model for suppression of forest fires. Man. Sci. 10(4):

760-766.

Parlar, M. 1983. Optimal forest fire control with limited 

reinforcements. Opt. Cont. Appl. Meth. 4(2): 185-191.

Parlar, M., and R.G. Vickson. 1982. Optimal forest fire 

control: an extension of parks’ model. For. Sci. 28(2):

345-355.

Pyne, S.J.; P.L. Andrews; R.D. Laven, 1996. Introduction 

to wildland fire. Wiley, New York. 769 p.

Robinson, R.; W. Hansen; K. Orth. 1995. Evaluation of 

Environmental Investments Procedures Manual. IWR

Report 95-R-1. U.S. Army Corps of Engineers, Institute

for Water Resources. Alexandria, VA. 81 p.

Sparhawk, W.N. 1925. The use of liability ratings in 

planning forest fire protection. J. Ag. Res. 8: 693-762.

USDA Forest Service. 1991. National fire management 

analysis systems users’ guide of the initial action assessment. (FPL-IAA 2.3). USDA For. Serv., Washington,

D.C.

187
This page is intentionally left blank.
In: Bevers, Michael; Barrett, Tara M., comps. 2005. Systems Analysis in Forest Resources: Proceedings of the 2003 Symposium. General Technical Report

PNW-GTR-656. Portland, OR: U.S. Department of Agriculture, Forest Service, Pacific Northwest Research Station.

1 A.E. Wilson is a Resource Information Specialist and 2 M.R. Wiitala is an Operations Research Specialist at the Pacific Southwest Research Station, 

USDA Forest Service, 620 SW Main Street, Suite 400, Portland Oregon, 97205; e-mail: aewilson@fs.fed.us.

Systems Analysis in Forest Resources: Proceedings of the 2003 Symposium

AN EMPIRICALLY BASED MODEL

FOR ESTIMATING WILDFIRE 

SUPPRESSION RESOURCE RESPONSE TIMES

Andrew E. Wilson1 and Marc R. Wiitala2

INTRODUCTION

Planning for wildfire suppression is one of many ways

to help increase the efficiency of a forest fire-fighting

organization. Pre-suppression planning and analysis involves

assessing how many resources to use, where to locate them,

and how they will be dispatched. While there are several

approaches to planning a pre-suppression organization, most

account for and place great importance on the time it takes

the suppression forces to move from their location to a fire

incident. Mees (1986) found response time to be one of the

most influential factors in his study of suppression resource

base locations. 

In that same study, Mees also describes suppression

resource response time as having two major components:

getaway time and travel time. For the initial response to a

fire, the getaway time is the time between a dispatcher’s

request for a resource to respond and the beginning of

transport to the fire. Travel time begins with the start of

transport and ends with arrival at the fire. 

ABSTRACT

Suppression resource response times are an important component in planning models used to evaluate the effectiveness

of wildfire preparedness organizations. Most fire planning models rely on subjective estimates of travel times between predetermined points on the landscape. This paper describes a raster-based travel time model that provides greater accuracy,

objectivity, and flexibility in estimating travel times between points on the landscape. The methodology used to create the

travel time prediction model integrates existing GIS data with empirical travel data collected using GPS. Although the need

to develop this rose from development of WIRAS (Wildland Initial Response Assessment System), it will be useful to many

wildland fire pre-suppression planning efforts where there is a relationship between the time it takes fire suppression resources

to respond and the final size and resultant costs of suppression.

189

Different types of fire fighting resources travel differently over the landscape. Travel time for aerial resources is

easily determined by assuming flights are made in straight

lines between locations. Accounting for terrain, transportation infrastructure, equipment capabilities, and other uncertainties complicates the computation of ground resource

travel times. 

The need for improved techniques for predicting travel

times arose during the development of the Wildland Initial

Response Assessment System; also know as WIRAS (Wiitala

and others, this proceedings). WIRAS is a simulation system designed for use in planning and analyzing strategic

pre-suppression decisions, such as quantities and locations

of fire fighting resources and deployment rules. To properly assess the effectiveness of the system, the response times

for aerial and ground resources are evaluated on an eventby-event basis as the scenario unfolds. WIRAS uses these

response times in two primary ways: first, to select the

quickest arriving fire fighting resources from a pool of
available resources, and secondly to simulate the proper

delay before individual resources start fireline production. 

Within WIRAS, there are three possible location situations where a resource would be available to respond to a

fire (fig. 1): 

• at a base such as guard station or district office, 

• at a fire, or 

• at an intermediate location along a dispatch route. 

This last situation arises when a resource in transit to 

a fire becomes available to respond to another fire. This

could happen when the first resource arriving on the fire

scene determines some of the resources en route to the fire

are not needed for fire containment. The surplus resources

are then available to respond to other incidents or, if not

needed, return to their base. 

WIRAS also simulates the demobilization of resources

from fires back to bases. Travel times for these movements

between bases and fires as well as those from fire to fire

are computed prior to simulation. However, the travel times

of in-transit resources becoming available for reassignment

must be computed during the simulation from their calculated locations at the time of reassignment. 

WIRAS provides the capability to address the variable

nature of response times. These variations arise from the

imperfect information on fire location that results in additional search time and other random time-consuming

events related to getaway and travel times. 

To meet the unique needs of WIRAS, a development

effort was started to improve travel time prediction capability for the simulation model. This paper describes an 

on-going study to develop a model with greater accuracy,

objectivity, and flexibility in estimating ground resource

travel times between points on the landscape. The methodology used to create the travel time prediction model integrates existing GIS data with empirical travel data

collected using GPS.

RELEVANT TRAVEL TIME STUDIES

A prominent study of vehicle travel time on forest roads

is the Logging Road Handbook (Byrne and others 1960),

verified by Moll and Copstead in 1996, which focused on

log truck operations. To predict speeds, the study used

engineering characteristics of the road (such as slippage

due to surface type and curvature derived through surveying particular stretches of roads) and the physical

characteristics of the vehicles (including gross vehicle

weight, horsepower, tire pressure and weight of the load).

This methodology required intensive attribution of equipment and travel routes, and in the context of forest road

engineering, it was able to ignore the problem of off-road

travel by foot. 

Mees (1978) developed a model of firefighting resource

allocation as a network problem (nodes connected by lines

with travel times) using linear programming optimization

to determine minimum travel time from an origin (a base)

to destination (a fire). The network method works well when

travel is restricted to linear features, like roads and trails.

Mees resolved the problem of fires not falling on the transportation network by creating pseudo-nodes at the points

on the inter-node connections closest to the incident and

connecting these points to the incident with perpendicular

segments using a fixed ground speed. There was not a discussion on how this heuristic worked in the face of terrain

barriers and the methodology of determining inter-node

times was not disclosed.

190

Figure 1—Potential origins when responding to fire 2 include a base (A),

another fire (1) or an intermediate point (@) where the resource, initially

responding to fire 1 is reassigned to fire 2.
Travel speed by fire fighters on foot has been reported

in relation to escaping from fires (Dakin 2002; and Butler

and others 2000). While the former study focused on short

(250 m) courses, the latter used historical records of fire

fighter travel rates at two well-documented incidents to

make estimates of sustainable travel rates over broad slope

categories. Howley and Franks (1992), in the context of

sports medicine, used human physiology to determine

walking speeds over varying slopes based on a fixed rate 

of energy output. A comparison in figure 2 of the physiological walking model derived from this work and the sustained travel rates that Butler and others found shows the

results are similar.

The intensive development of physical characteristics

required for the engineering approach to this problem suggested by the Logging Road Handbook did provide insight

into important factors in travel time, the time and costs of

implementation were not compatible with the broad nature

of fire pre-suppression planning across forest or landscape

levels. The simple approach suggested by the linear programming approach forwarded by Mees did not adequately

address the problems of access across varied terrain where

the shortest distance is not necessarily the quickest distance.

STUDY AREA

The Umatilla National Forest was selected for the site

of the study. The forest occupies about 1.4 million acres in

an area 110 miles by 130 miles in the Blue Mountains of

Northeast Oregon and Southwest Washington. It has four

districts, two in the steeper, wetter north half of the forest

and two in the drier, more rolling south half. Outside wilderness and inventoried roadless areas, fires are moderately

accessible with 64 percent of the forest within 4 mile of an

open road, a distance within which fires can generally be

attacked by an engine hose lay. 

The Umatilla fire season has about 136 fires in a year

with the majority of them spread across July, August and

September. The forest directly employs about 20 engine

and hand crews. Cooperators and contractors provide additional sources of dozers, tenders, engines and hand crews.

A typical initial attack resource deployment at a district

includes three engines and two or three hand crews. The

Pendleton Interagency Communications Center, which dispatches fire suppression resources for the Umatilla National

Forest and cooperatively protected lands, provided much of

information on initial attack resource assignment required

by this study.

METHODOLOGY

The study’s methods consist of two parts: the development of a raster-based model to calculate travel times and

an empirical analysis of travel times. To predict travel times,

a raster-based rather than network approach was selected to

maximize computational speed and facilitate travel time

calculations anywhere across the landscape. This model

uses a travel resistance surface that incorporates both vehicle and foot travel. The travel time model incorporates a

minimum cost path algorithm (Hatfield and others, this

proceedings), which can calculate in a few seconds the

responses of twenty or more resources spread across the

forest to a single fire. Using the raster approach allows a

better representation of the walk-in portion of the response

time that can be affected significantly by terrain.

A major component of the travel time study is the

empirical collection and analysis of travel time data.

Empirical travel times provide a basis for testing accuracy

of the minimum cost-path algorithm time calculations and

adjusting the model for identified systematic bias. The

deterministic predictions of the raster-based travel time

model are expected to be different from actual travel times.

As noted earlier, these statistical differences arise from

variable search times due to imperfect initial fire reports

and other delays experienced during an initial attack

response. The empirical travel time data will be used to

191

Figure 2—A comparison of walking model estimates of a firefighter’s sustainable walking speed with sustainable walking speeds found in Butler

and others (2000).
quantify these underlying statistical processes for incorporation into the travel time prediction mode for use in

WIRAS. 

To acquire the necessary data to develop the travel time

surface and support our verification and statistical analysis,

we chose to collect empirical data on response speeds over

the fire season. Three different data sets were collected:

position-velocity logging, dispatch records, and radio transmission position logs. 

Position-velocity logging units were installed in engine

and crew vehicles and kept a record of the vehicles velocity based on time and position. The data (the time, the location of the vehicle and its speed) were determined within

the unit with GPS and recorded at 15-second intervals as

long as the vehicle was moving, regardless of whether the

vehicle was responding to an initial attack incident or not.

The data was stored on removable storage media and manually downloaded periodically.

To determine the data points where the vehicle was

responding to a dispatch, we gathered the dispatch record

information on resources ordered for each incident, including the unit identifier, when the dispatch order was sent and

where the unit was, and when the unit arrived at the incident or was reassigned/released. The unit identifier and the

dispatch start and stop times were used as key fields to filter the excess data from the logged velocity points. By

examining the time series of logged position-velocity data

points associated with each run, we could detect and quantify travel delays for resources responding to incidents. The

magnitude of a delay is measured by the difference in

expected response time and actual response time. This

information is used to determine the nature of the variability of these delays. 

There is always the possibility of confusion over which

resources were involved in a dispatch and where they were

located at the time of dispatch (especially in the reassignment case). We attempted to reduce the possibility of this

confusion by augmenting the resource position information

in the dispatch record with a resource position data set

based on a radio-based real-time position tracking system.

This position tracking system provided a GPS-derived

coordinate/time pair for each radio transmission from the

fire fighting resource that helped precisely determine unit

initial attack response start and stop times and positions.

Since the position information was displayed in real-time

in the coordination center, it improved dispatch knowledge

of resources and subsequent record keeping. 

Once the velocity data was associated with initial attack

incidents, the position information was used to determine

the attributes of the predictors. Based on forest engineering

studies like Byrne and others (1960) and discussions with

engineering staff, we selected a set of available and fairly

complete data elements from the Forest Service database

(INFRA) that stores travel route information linked to a GIS

route system. While the INFRA database had relatively fine

attribute detail, it was sparse beyond the boundaries of forest Service administrative boundaries. We augmented the

INFRA data these areas with information from the GTRN

transportation database (BLM 2002) developed by the

Bureau of Land Management. The other needed piece of

spatial information was the terrain database. We resampled

the 10-meter database from the national elevation dataset

(USGS 2003) to 30-meter resolution, as we believed that

this was the size of ground relevant to the movements at

the forest scale. The specific road attributes examined

included surface type (paved, gravel, improved native and

native), operational maintenance level (closed, high clearance, suitable for passenger cars, moderate degree of user

comfort and high degree of user comfort) and two derived

attributes, Grade and Alignment. Deriving these two attributes required segmenting the linear road dataset to a finer

resolution detail – rather than segments that connect intersections, the segments were split so that the maximum

length was 4-mile. This length was chosen to approximate

the average area of influence the attributes would have on 

a vehicle’s speed and, given the rate of sampling, to ensure

that pieces would not be missed in normal transit. To derive

Grade, the elevation difference between each segment’s

endpoints was divided by the connecting road length segment. The Alignment measure was determined by dividing

the straight-line distance between two points by the actual

road distance depicted in the data. This resulted in a continuous measure of alignment, where the straighter the road

was, the closer to unity the alignment value would be.

The other portion of the trip to the fire involves walking

to the fire from where the vehicles are parked. Ideally,

empirical data of the entire trip could be collected much

like the velocity-position logging that was gathered on

vehicles. Unfortunately, a cost-effective method to collect

this type of data has not been discovered. Instead, the physiological-based model introduced above was applied and

will be tested against the empirical data from the dispatch

record and radio transmission logs. This model uses the

assumption of constant energy expenditure based on energy

requirements for walking (Howley and Franks 1992). The

sum of vertical and horizontal net oxygen costs determines

the change in speed as slope varies. This simplifies to a

quadratic slope–speed relationship.

192
We assume that travel by vehicle will be at least as fast

as travel by foot, so the final model uses the vehicular travel

model where the roads exist, and the walking time model

elsewhere. The combined model, measured in minutes per

unit measure, can then be summed to determine travel time

across the landscape, the basis for determining the minimum

duration path. 

RESULTS

Position-velocity loggers were installed in initial attack

vehicles in the south half of the forest in late-summer 2002

and the dispatch records associated with these vehicles were

coded into the database. Difficulties in the engineering of

the radio-based tracking system delayed implementation of

those units until spring, 2003. All resources being examined

had both position-velocity logging units and the radio position tracking hardware installed in them. The 2003 season

will have the full complement of data for examination when

it is complete in mid-November, allowing a full analysis of

the data. 

Even though there was limited data from the 2002 season, the procedures were tested to ensure that the data could

be filtered and analyzed to allow preliminary examination

of relationships. Given this preliminary nature of the data,

we saw that a primary factor in vehicle speed was the surface type (fig. 3) – most of the paved roads were highways

with resultant high average observed speed of 44 mi/h with

subsequent drops in speed as the vehicle moved to lower

quality surfaces of gravel (24 mi/h) and improved and

unimproved native materials (12 and 10 mi/h). 

Examination of the data at the tails of the curves showed

that there were errors induced by lack of coincident positioning between the observed GPS data and the GIS data.

Many of these errors occurred at intersections between road

classes. In the GIS, roads are lines without widths. Even a

very accurate observation near an intersection would not

fall precisely on either road in the GIS and would have

some probability of being closer to the intersect road. As

the accuracy of the position decreases, the probability of

being closer to the intersecting road increases. While it is

tempting to remove all points near intersections to remove

this error, potentially there are factors affecting speed at

these same intersections as well. 

Other problems associated with the GPS measurements

included occasional incorrect positions (usually not more

than 300 yards from the actual location) when a poor set 

of satellites was chosen from the constellation by the GPS

unit, causing one or more locations to be recorded incorrectly. Resolving these errors required a thorough manual

examination of the time series of data points to catch these

errors. An even rarer error (about one in 10,000) was the

recording of an anomalous velocity attribute. These were

always single incidents in the data stream and were identified through looking at a moving window of statistics in

the time series. Another source of error was the dispatch

record. It had to be thoroughly examined and compared to

the observed data to ensure that it was correct. Many of the

problems occurred during days when there were multiple

fires and were caused by the combination of more resources

to track and additional dispatchers that were not familiar

with the area.

Although the data for the first year was from the south

half of the forest and only included the last part of the fire

season, we were able to observe some apparent relationships. Once the data was classed on surface type, it appears

that alignment was significant in determining speeds on the

paved roads, while it had little impact on the improved native

roads. There was a high degree of correlation between surface type and maintenance level when looking at all categories, however, the closed road class in maintenance level

was not as correlated and provided predictive power on the

non-paved roads.

193

Figure 3—Frequency distribution of speeds observed on different road

surface types.
It should be noted that using the raster approach as

described above might not be an appropriate tool for realtime routing of emergency response services. Individual

routes determined optimal with the model may, in reality,

be infeasible or may have longer durations than the model

specifies. As well, complex traffic controls, such as oneway roads and restricted intersections are more difficult to

address. However, the model works well within this strategic planning realm of wildfire suppression.

CONCLUSION

This paper describes development of a model to estimate the response time of initial attack resources for use in

strategic planning. It discusses data collected and used in

developing the model. Problems with the GPS-derived data

are pointed out and methods to remove errors are posited.

Finally, preliminary results are presented and indicate that

this will be an effective improvement in the determination

of response times in analysis of wildland fire pre-suppression planning.

LITERATURE CITED

Bureau of Land Management, Oregon State Office. 

2002. BLM ground transportation (roads & trails).

Retrieved April 2, 2002 from http://www.or.blm.gov/

gis/data/catalog/dataset.asp?cid=39

Butler, B.W.; Cohen, J.D.; Putnam, T.; Bartlette, R.A.; 

Bradshaw, L.S. 2000. A method for evaluating the

effectiveness of firefighter escape routes. Pages 42-45

in Butler, B.W. and K. Shannon (eds.) Proceedings from

the International Wildfire Safety Summit (Oct 10-12,

2000, Edmonton, Alberta). Int. Assoc. Wildland Fire,

Montana City, Mont.

Byrne, J.; Nelson, R.; Googins, P. 1960. Logging road 

handbook: the effect of road design on hauling costs.

Agriculture Handbook No. 183. U.S. Department of

Agriculture, Forest Service.

Dakin, G. 2002. Ground rates of travel by fire crews 

using escape routes: an interim report. Advantage 3(15).

Forest Engineering Research Institute of Canada. 

Howley, E.T.; Franks, B.D. 1992. Energy costs of activity. 

In Health fitness: instructor’s handbook. 2nd edition.

Human Kinetics, Champaign, Illinois, pp. 131-139.

Mees, R.M. 1978. Computing arrival times of firefighting 

resources for initial attack. Gen. Tech. Rep. PSW-27.

U.S. Department of Agriculture, Forest Service, Pacific

Southwest Forest and Range Experiment Station.

Mees, R.M. 1986. Locating suppression resources by 

travel time to wildfires. Res. Note PSW-387. U.S.

Department of Agriculture, Forest Service, Pacific

Southwest Forest and Range Experiment Station.

Moll, J.; Copstead, R. 1996. Travel time models for 

forest roads: a verification of the Forest Service 

logging road handbook. Pub. 677 1202 – STDC. 

U.S. Department of Agriculture, Forest Service, 

San Dimas Technology Development Center. San

Dimas, California.

U.S. Geological Survey, EROS Data Center. 2003.

National elevation dataset. Retrieved February 2003

from http://gisdata.usgs.net/ned/default.asp.

194
In: Bevers, Michael; Barrett, Tara M., comps. 2005. Systems Analysis in Forest Resources: Proceedings of the 2003 Symposium. General Technical Report

PNW-GTR-656. Portland, OR: U.S. Department of Agriculture, Forest Service, Pacific Northwest Research Station.

1 Fried, Christensen, Weyermann, and Hiserote are, respectively, team leader, forester, computer specialist, and computer specialist in Forest Inventory

and Analysis, 2 Barbour is program manager in Focused Science Delivery, and 3 Fight and Pinjuv are, respectively team leader and forester in Human and

Natural Resources Interactions, USDA Forest Service, Pacific Northwest Research Station, Portland Forestry Sciences Laboratory, Portland, OR 97205.

Systems Analysis in Forest Resources: Proceedings of the 2003 Symposium

MODELING OPPORTUNITIES AND FEASIBILITY

OF SITING WOOD-FIRED ELECTRICAL GENERATING

FACILITIES TO FACILITATE LANDSCAPE-SCALE FUEL

TREATMENT WITH FIA BIOSUM

Jeremy S. Fried1, Glenn Christensen1, Dale Weyermann1, R. Jamie Barbour2,

Roger Fight3, Bruce Hiserote1, and Guy Pinjuv3

INTRODUCTION

As large, stand-replacing fires have become common

over the past few years in the western United States, interest in assessing the potential for landscape-scale fuel treatment has intensified. Doubts have grown about the feasibility

and/or wisdom of widespread use of prescribed fire due to

concerns about air quality, liability, narrow windows of

opportunity to implement treatments, and potentially undesirable fire effects. Fuel treatment has for most forestry

professionals and most of the public come to be virtually

synonymous with mechanical removal or thinning of the

ABSTRACT

Utilization of small diameter trees is viewed by many as the key to making landscape-scale fuel treatment financially

feasible. But little capacity currently exists for utilizing such material and capacity of sufficient scale to have a significant

impact on the economics of small diameter removals will only be added if predictable feedstocks can be assured. The FIA

BioSum modeling framework that incorporates Forest Inventory and Analysis (FIA) plot data, a transportation cost model, a

treatment cost accounting module, a log valuation model, and a crown fire hazard evaluator was applied to a 28 million acre

study area containing 6200 FIA plots spanning the Eastern Cascades, Southern Cascades, Klamath Mountains and Modoc

Plateau ecosections of western Oregon and northern California. Up to nine fuel treatment prescriptions with a high likelihood of producing a substantial reduction in crown fire hazard were simulated for each plot, and 221 potential biomass processing sites were considered. With four 50 MW biomass-fueled power plants strategically distributed over the study area,

up to 5.3 million acres could be effectively treated with net revenue of 2.6 billion dollars, a merchantable yield of 9.5 billion cubic feet, and a biomass yield of 79 million green tons, if net-revenue maximizing fuel treatments are selected. If merchantable volume minimizing treatments are selected instead for these 5.3 million acres, net revenue would be negative 2.6

billion dollars, merchantable yield would be 3.6 billion cubic feet and biomass yield would be 75 million green tons. With

the constraint that every acre generate positive net revenue, only 2.6 million acres would be treated, even if the net revenue

maximizing treatment is selected.

195

forests to reduce fire severity and the likelihood of standreplacing fire, especially since the advent and passage of

the Healthy Forest Restoration Act of 2003. 

The conventional wisdom has been that in order to be

effective, such treatments would require the removal of

large numbers of small stems, at considerable cost, but that

this harvested material would have little or no value. In part

to address this concern, but also with an eye towards promoting renewable energy options and increasing employment opportunities in rural forests, proposals have been

and continue to be floated in many parts of the U.S. to
develop markets for submerchantable-sized trees by promoting the construction of biomass processing facilities of

various capacities that would convert biomass to electrical

energy via direct combustion or through various intermediate pathways (for example, gasification). However, such

facilities require a sizable, up-front investment—one that is

unlikely to be made by the private sector without some

confidence that there will be an adequate source of supply over a long enough period for the investment to generate a positive financial return. 

Fire and forest planners and managers, rural community

economic development staff, and potential investors in biomass processing capacity believe they could benefit from

knowledge of where the “hot spots” of potential biomass

supply are to be found and the kinds of materials, both submerchantable and merchantable, that could be reasonably

expected to flow from landscape-scale fuel treatments. This

analysis was undertaken to identify such hot spots, but has

evolved, thanks to wide-ranging discussions with diverse

clients, to encompass a great deal more—including comparison of a number of fuel treatment prescriptions, assessment of the economic feasibility of fuel treatments with a

complete accounting of harvest and haul costs, and a modelbased characterization the fire risk reduction accomplishment of such a fuel treatment program. Questions we sought

to address progressed from “how much biomass is out

there” and “where is the best place to site a biomass processing facility” to “how much biomass is legally and economically accessible”, “how much will treatments cost”,

and “will a subsidy help”. 

We constructed a geographically explicit analytic framework for assessing and summarizing biomass production

opportunities, dubbed FIA BioSum (fig. 1). We used publicly available, field-collected forest inventory plot data,

publicly available road and ownership GIS layers, a suite

of publicly available models, and sets of assumptions,

parameters and decision rules developed in consultation

with local fire, fuels, silviculture, and logging experts as

well as biomass plant operators to address the feasibility

issue for a 28 million acre, 4 ecosection study area in

Oregon and California. These ecosections (Klamath, Modoc

Plateau, southern Cascades, and eastern Cascades) were

selected because current fire regime condition class maps

of the U.S. show much of the area in these ecosections as

in condition class 3—lands that have significantly altered

vegetation composition, diversity, and structure due to

altered fire return intervals so that they verge on the greatest risk of ecological collapse due to loss of key ecosystem

components from fire and are thus likely to receive high

priority for fuel treatment (USDA, USDI 2002). 

METHODS

Forest inventory data for this analysis were drawn from

six inventories undertaken at various times during the

1990s by the USDA Forest Service PNW Research

Station’s Forest Inventory and Analysis Program, USDA

Forest Service Region 5, USDA Forest Service Region 6,

and the Bureau of Land Management. While there were

numerous design differences among these inventories, each

was a statistically representative sample of a portion of the

total landscape and included measurements of tree attributes such as diameter, height and species with this information compiled to produce plot level estimates of volume,

biomass, basal area and density, for example. A total of

6200 field plots representing 22.2 million acres of forested

land fell within the study area boundary. This plot set was

culled to remove from further analysis plots that were evaluated in the field as non-forest or located in designated

wilderness, natural areas, parks, preserves, monuments,

national recreation areas, national wildlife refuges, and

inventoried roadless areas that had been under consideration for possible protection in the 1990s. We also omitted

plots on steep (> 40 percent) slopes that were more than

2000 feet from the nearest mapped road on the grounds

that current skyline harvesting systems cannot economically reach beyond this limit, and the cost of building new

roads to bring such equipment closer to the plot would

likely be prohibitive and is in any case beyond the scope 

of this analysis. Finally, we omitted plots containing no

trees over 5 inches d.b.h.. Due to unforeseen data management problems, the 20 percent of the plots that contained

196

Figure 1—Flow diagram of the FIA BioSum modeling framework.
more than one condition (in other words, the plot contained

road or non-forest land in addition to forest or contained

more than one type, size or density of forest) could not be

processed in time for this analysis, so the results presented

here cannot be considered comprehensive over all lands

and surely underestimate costs, revenues, area treated and

product flows by an unknown amount if considering the

entire landscape; still, these results do apply to the 80 percent of the landscape represented by single condition plots.

The net result of all these reductions is a set of plots that

represents about 10.4 million acres of forest land—8.2 million acres of this federally administered and 2.2 million

acres privately owned. We did not account for special use

areas with harvesting restrictions such as riparian buffers or

late succesional reserves designated under the Northwest

Forest Plan. Nor did we include restrictions such as those

found in the Healthy Forest Restoration Act of 2003 that

attempt to focus treatments within a fixed distance of communities. Either of these types of management restrictions

could be accounted for in future analyses if they are accurately mapped.

The plots used in the analysis were loaded into the

Forest Vegetation Simulator (FVS) (Stage 1973; Wykoff

and others 1982) to simulate fuel treatment prescriptions

and, via its Fire and Fuels Extension (FFE) (Reinhardt and

Crookston 2003), the likely change in fire hazard these

treatments would produce. The two FVS variants applicable

to this study area, South-central Oregon and Northeastern

California (SORNEC) and East Cascades (EC), were used

for modeling and their FFEs were used to derive several

fire-related stand attributes and ultimately compute indices

of crown fire potential for each plot, specifically, torching

index (TI) and crowning index (CI). TI represents the wind

speed at which fire could be expected to move from surface fuels into crown fuels and is highly influenced by vertical stand structure (ladder fuels) and height to crown base

(derived from crown ratio); CI is the wind speed at which a

crown fire could be expected to be sustained and is heavily

influenced by crown bulk density. For the purposes of this

study, increases in CI or TI brought about by fuels treatment were assumed to reduce fire hazard.

Nine fuel treatment prescriptions representing two treatment approaches were developed in consultation with fire

and fuels specialists and applied to all plots for which they

were valid (table 1). For example, a prescription calling for

a post-treatment residual basal area of 125 ft2/ac could not

be applied to a plot containing only 80 ft2/ac. Five of the

treatments were designed for density reduction (personal

communication with John Gerritsma on 3 May 2002), with

a primary focus of reducing the propagation of a crown

fire, and involved thinning proportionately across all diameter classes to a target residual basal area, with the proviso

that 70 percent of the basal area cut be from trees that are

5.5-14.5” d.b.h. with the remainder from trees >14.5 inches. The other three were designed as ladder fuels reduction

treatments (personal communication with John Szymoniak,

6 May 2002), with a primary focus of reducing the initiation of crown fire, that thinned from below (>5.5” d.b.h.)

to a residual basal area target. The treatments had a range

of residual basal area targets (60-125 ft2/ac) and maximum

acceptable diameter for cut trees (10” to no limit). For both

kinds of treatments, but mainly for ladder fuels reduction,

if the maximum diameter limit was reached before the

residual basal area target, then the latter was not achieved.

For all prescriptions and on all plots, trees less than 3.5”

d.b.h. were cut into pieces and scattered, and trees 3.5-5.5”

were cut and scattered on steep slopes (>40 percent) and

harvested and collected to the landing as biomass on gentle

slopes. Because the most aggressive treatment had a residual basal area of 60 ft2/ac, plots with less basal area were

excluded, leaving 6.9 million acres represented in the

analysis.

197

Table 1—Fuel treatment prescriptions applied in a test of the FIA BioSum analytic framework.

Density reduction Fuel reduction

Prescription A B C D E F G H J

Target residual 

basal area ft2 125 125 90 90 90 80 60 60 60

Maximum 

d.b.h. (in.)

allowed for 

harvested

trees 21 none 21 16 none 21 21 none 10
These prescriptions were applied in FVS-generated “cut

lists” for each plot and, with post-processing via various

program scripts, numerous plot and tree attributes used at

later stages in the analysis (for example, slope, number of

stems harvested, average stem size harvested, and volume

and biomass by species and size class) were computed and

stored. Post-treatment estimates of TI and CI were computed

with FFE for each prescription for each plot. 

Data extracted from the FVS output was processed with

STHARVEST (Fight and others 2003), a spreadsheet model

composed of regressions and look-up tables for logging

cost components derived from empirical data on small timber sales. STHARVEST requires specification of assumptions regarding logging system (for example, whole tree,

cut-to-length, cable), range of tree diameters to be included,

volume per acre, and disposition of residue. Miscellaneous

costs (for example, “brush-cutting”, waterbars) were also

accounted for. Whole tree logging systems were assumed

for slopes ≤40 percent (except for trees too large to be handled without bucking); cable systems were used on slopes

>40 percent, with manual felling, bucking and limbing in

the woods (in other words, limbs and tops not recovered

for biomass utilization). For each plot and prescription,

STHARVEST provided an estimated on-site cost of implementing the prescription (in other words, from stump to

truck). Harvested material was allocated to one of two categories: Merchantable, which consisted of boles of trees

>7.0” d.b.h. to a 5” diameter top, and Biomass, which consisted of trees brought to the landing that were 3.5 – 7.0”

d.b.h., the limbs and tops of merchantable trees, and all

harvested hardwoods. Delivered biomass was valued at $18

per green ton, the going rate circa 2003 in places like

northern California where multiple biomass facilities compete for material, and merchantable material values were

assigned from a look-up table based on species group and

tree diameter. Plot-treatment combinations that generated

less than 300 ft3/ac of total volume (biomass and merchantable combined) were deemed unrealistic and discarded; for some plots, no treatment cleared this hurdle and

these plots were excluded, leaving a set of plots representing 5.7 million acres amenable to one or more of the fuel

treatment prescriptions. 

Because this analysis targets fuel treatments that reduce

the risk of stand-replacement fire, only treatments that are

effective in achieving this goal are included. Reaching a

consensus among fire and fuel managers on which crown

fire potential index is most critical and on how much an

index must change for the effect to be significant proved

elusive. Consequently, we developed, as a working criterion to demonstrate FIA BioSum’s analytic potential, an

arbitrary, minimum threshold of fuel treatment effectiveness: a 20 mph improvement (increase) in either TI or CI

with no reduction in the other. Treatments that fell short of

this threshold were deemed ineffective and discarded; this

resulted in some plots being untreatable and reduced the

area in the analysis to 5.4 million acres. 

Because the problems we tackled are strongly associated with location (where can treatment be feasibly applied,

where to site biomass processing facilities) and our plans to

eventually integrate this analysis with maps of the wildland

urban interface, we chose to pursue a spatially explicit analytical framework that would account for differences in

haul costs and identify locations with a sufficient potential

accumulation of biomass to justify investment in a processing facility. An approximately 20x20 km grid of 221 potential processing sites (psites) was established over the study

area, approximate in the sense that we relocated psites that

fell on public lands. GIS road layers obtained from BLM,

USDA Forest Service Region 5, USGS and others were

combined, edge-matched, cleaned of gross anomalies and

massaged to produce a study area-wide, topologically

imperfect vector GIS road coverage with each road segment

attributed as to likely rated speed. This road coverage was

then tessellated into 250 m grid cells, with each cell’s value

set to the cost per ton-mile of traversing the cell on the

fastest (lowest cost of transit) road in that cell, resulting in

an impedance surface that formed the basis of haul cost

calculations for FIA BioSum. Cells containing no roads

were assigned infinite cost of transit, and plots occurring in

such cells were “moved” to the nearest road for purposes

of haul cost calculation, retaining the distance moved as 

an input to the yarding cost calculations used by STHARVEST. For each psite, a cost accumulation grid was generated in Arc/Info, and overlaid on the plot grid to provide

haul cost in dollars per ton to that psite from every plot in

the study area. The haul costs were combined with the outputs from FVS, FFE, and STHARVEST to create a database

that stored the biomass and merchantable yields, harvest

and haul costs, gross and net revenues, and change in TI

and CI associated with every combination of plot, psite 

and prescription. 

To test the simulation framework and develop preliminary estimates of biomass processing opportunities for this

broad region, we evaluated seven fuel treatment policy scenarios under the assumption of the construction of 4 biomass processing facilities capable of generating 50MW each.

The locations were selected from among the sites with the

best accumulation potential – of biomass, merchantable

volume, net revenue, and acres treated—with the additional

requirement that they be geographically separated (in other

198
words, one per ecosection) so as to have little influence on

one anothers’ markets. Each of the 1556 plots (representing

the 5.4 million acres that could be treated effectively) was

allocated to the psite with the lowest haul cost: Grants Pass

OR, Burney CA, Klamath Falls OR, or Bend OR, and the

rest of the analysis proceeds with those allocations (fig. 2).

Four of the treatment scenarios were: treat all treatable plots

and select the prescription for each plot that 1) maximizes

net revenue, or 2) maximizes improvement in TI, or 3)

minimizes the merchantable material removed, or 4) maximizes improvement in CI. The other three scenarios 1A,

2A and 3A, corresponded to scenarios 1, 2 and 3, except

that only plots with non-negative net revenue would be

treated. Scenario 1 was designed to make as much profit as

possible subject to the non-trivial constraint that all treatments were designed with an eye towards reducing crown

fire potential, not maximizing economic return or even

facilitating a positive cash flow. Scenarios 2 and 4 represent

a no-holds-barred policy of reducing fire hazard as the

paramount objective, while scenario 3 might be more

appealing to groups or individuals who favor the achievement of some measure of fuel reduction but do not support

the removal of sawtimber-sized trees and or the generation

of profits from fuel treatment activity.

RESULTS

The fuel treatment policy scenarios 1-3 lead to quite

different distributions of treatments. The most frequently

selected treatments for scenario 1 were H, G, and J, for

scenario 2 were G, H, and F, and for scenario 3 were J, G

and F (fig. 3). Given the substantial differences in the treatments selected, we expected to see substantial differences

in the aggregate results; however, the most substantial differences among scenarios were in net revenue, while the

biomass yield was rather robust under this range of scenarios (table 2). Under scenario 1, net revenue was positive

for every psite when all effectively treatable plots were

treated. Even when selecting prescriptions that minimize

merchantable yield (scenario 3), there would still be 3.6

billion ft3 of merchantable volume produced, though the

aggregate net revenue under this scenario is extremely negative. Under every scenario, the acres feeding material to

the psite at Grants Pass generated the greatest merchantable

volume and the greatest biomass. 

With the added constraint that, as proposed in the

Healthy Forests Initiative, “every acre pays its own way”,

the area that could be treated would be reduced to no more

than 12 percent of the forested landscape (table 3) in scenarios 1A, 2A and 3A. Of the four sites, Grants Pass again

results in the greatest amount of net revenue, acres treated,

199

Figure 2—The Oregon-California study area with ecosections, accessible

plots and the four potential processing sites considered in this analysis.

Figure 3—Number of acres for which each fuel treatment prescription

was selected; numbers below prescription letters are target residual basal

area in ft2/upper diameter limit in inches.
biomass yield, and merchantable yield under scenario 1A.

Scenario 3A treats about the same area as scenarios 1A and

2A, but produces less than half the net revenue, largely

because there is much less revenue from merchantable 

volume in the prescriptions selected for this scenario that

could be used to offset treatment costs. 

All scenarios produced quite different opportunities in

terms of the potential longevity of a biomass plant (table

4). Under scenario 2, sufficient biomass could be collected

to operate 50 MW generating stations at Grants Pass and

Burney for about 50 years, and in all likelihood, much

longer given that this was a static analysis (that did not

account for regeneration and continued growth of the residual trees and involved no harvest scheduling component),

and preliminary analysis of post-treatment plots modeled

forward with FVS indicates a need to re-treat as soon as 20

years after the initial treatment for fire hazard reduction to

be maintained. Requiring every acre to pay its own way

dramatically cuts the projected longevity of biomass supply, and scenario 3A (minimum merchantable yield and

nonnegative net revenue) would appear to all but eliminate

the opportunity to attract capital investment for the construction of new biomass generating capacity except at

Grants Pass. 

200

Table 2—Aggregate net revenue, area treated, and yields of merchantable volume and submerchantable biomass by potential processing site for three, fuel treatment policy scenarios 

in which all effectively treatable plots are treated.

Scenario 1: Net Acres Merchantable

Maximize net revenue revenue treated volume Biomass

(Millions of (Millions (Millions

Processing site dollars) (Thousands) of cubic feet) of tons)

Burney 542 1,546 2,635 24

Klamath Falls 593 1,101 1,741 14

Bend 343 741 998 12

Grants Pass 1,162 1,962 4,153 30

Total 2,640 5,351 9,527 79

Scenario 2:

Maximize torching index Net Acres Merchantable

improvement revenue treated volume Biomass

(Millions of (Millions (Millions

Processing site dollars) (Thousands) of cubic feet) of tons)

Burney -12 1,546 2,256 30

Klamath Falls 415 1,101 1,602 15

Bend 144 741 822 11

Grants Pass 376 1,962 3,940 37

Total 923 5,351 8,620 94

Scenario 3:

Minimize merchantable Net Acres Merchantable

yield revenue treated volume Biomass

(Millions of (Millions (Millions

Processing site dollars) (Thousands) of cubic feet) of tons)

Burney -795 1,546 974 21

Klamath Falls -157 1,101 773 12

Bend -148 741 384 9

Grants Pass -1,465 1,962 1,444 34

Total -2,565 5,351 3,574 75
Not surprisingly, average effectiveness of fuel treatments

in terms of both TI and CI was usually greatest under scenario 2, as shown in figure 4 for CI. What was surprising

was that the average improvement in both indices was

comparable between with scenario 1, in which the fuel

treatment with the greatest net revenue was selected.

Scenario 3 was substantially less effective on the plots

feeding every psite, with an average CI improvement of 20

to 40 mph, most likely because prescriptions that remove

less merchantable material also leave stands with higher

canopy bulk-density.

Table 5 summarizes the aggregate results for scenarios

1, 2 and 3 plus scenario 4 (maximize improvement in

crowning index). Scenario 4 has greater yields of both 

biomass and merchantable material, most likely because 

it favors treatments that reduce canopy bulk density (by

removing more of the medium-sized trees greater than 7”

d.b.h. and generally less than 21” d.b.h.). It is the additional merchantable component that accounts for the boost in

net revenue as compared with scenario 2. The pre-treatment TI and CI averaged 7 and 23, respectively, so on

average, all of the scenarios lead to substantial improvements in these indices. 

201

Table 3—Aggregate net revenue, area treated, and yields of merchantable volume and submerchantable biomass by potential processing site for three, alternate fuel treatment policy

scenarios in which only the effectively treatable plots that generate positive net revenue are

treated.

Scenario 1A: NR+ Net Acres Merchantable

Maximize net revenue revenue treated volume Biomass

(Millions of (Millions (Millions

Processing site dollars) (Thousands) of cubic feet) of tons)

Burney 1,112 735 2,017 11

Klamath Falls 891 524 1,482 7

Bend 536 335 833 8

Grants Pass 2,382 1,023 3,478 16

Total 4,921 2,616 7,809 42

Scenario 2A:

Maximize torching index Net Acres Merchantable

change, NR+ revenue treated volume Biomass

(Millions of (Millions (Millions

Processing site dollars) (Thousands) of cubic feet) of tons)

Burney 944 711 1,873 12

Klamath Falls 798 516 1,407 7

Bend 486 329 780 8

Grants Pass 2,184 1,014 3,295 17

Total 4,413 2,569 7,355 44

Scenario 3A:

Minimize merchantable Net Acres Merchantable

yield revenue treated volume Biomass

(Millions of (Millions (Millions

Processing site dollars) (Thousands) of cubic feet) of tons)

Burney 442 735 1,320 12

Klamath Falls 357 524 898 7

Bend 238 335 539 8

Grants Pass 942 1,023 2,182 21

Total 1,978 2,616 4,939 47
202

Table 4—Years of feedstock for a 50 MW biomass-based electrical generating plant 

consuming 1750 green tons per day of biomass in an operation 24/7.

Scenario Burney Klamath Falls Bend Grants Pass

1. Maximize net revenue 38 22 18 47

1A. Maximize net revenue, treat 

only stands with nonnegative NR 18 11 12 25

2. Maximize improvement in TI 47 24 18 58

2A. Maximize improvement in TI, 

treat only stands with nonnegative NR 19 11 12 27

3. Minimize merchantable yield 33 18 14 53

3A. Minimize merchantable yield, treat 

only stands with nonnegative NR 19 10 12 33

Table 5—Aggregate yield, net revenue, and improvement in crown fire potential (from a pre-treatment 

mean CI of 22 and TI of 7) for the four ecoregion study area for four fuel treatment policy scenarios, when 

all effectively treatable acres are treated.

Biomass Merchantable Net revenue Mean change Mean change 

(millions volume (millions of in CI in TI 

Scenario of tons) (millions of ft3) dollars) (miles per hour) (miles per hour)

1. Max. net rev. 79 9,527 2,639 63 96

2. Max ∆ TI 94 8,620 922 83 98

3. Min Merch. Vol. 75 3,574 -2,566 23 35

4. Max ∆ CI 93 9,765 1,817 80 114

Figure 4—Average improvement in FFE-predicted

crowing index, by proposed processing site to which

recovered material would be hauled, by scenario for

scenarios 1, 2 and 3.
DISCUSSION

There is plenty of biomass to supply four 50 MW power

plants for decades under the most aggressive scenarios (for

example, scenario 2), but supply under the most conservative scenarios (for example, 3A) would be far more limited.

Another alternative that we have not yet analyzed would be

to operate several smaller or less heavily capitalized power

plants for a shorter period of time in order to reduce haul

costs and make treatment feasible over a greater area.

By any basis of calculation, most of the material removed

in these treatments is merchantable and would fetch higher

prices than the $18 per green ton currently paid for biomass-sized material by power plant operators in northern

California. For scenario 1, 75 percent of the total tonnage

of wood removed is merchantable, with the remaining 25

percent in biomass-sized material composed of trees less

than 7 inches and the tops and limbs of larger trees. And of

the biomass-sized material, only 20 percent is in trees less

than 7 inches (i.e., only 5 percent of the total wood tonnage

removed). On a value basis, merchantable wood recovered

accounts for 90 percent of the total value recovered. This

helps explain why prescriptions that do not maximize net

revenue or that reduce the quantity of merchantable-sized

material recovered fare so poorly on a net revenue basis. 

Both harvest costs and haul costs are considerable,

averaging $1634 and $560 per acre, respectively, under

scenario 1, with the haul costs for biomass alone averaging

$127 per acre. And with the configuration of four processing sites we specified, haul costs alone averaged $8.53 per

green ton, nearly half the $18 per green ton that biomass

was assumed to fetch upon delivery, suggesting strongly

that biomass rarely pays its own way out of the woods,

given that harvest costs averaged $22 per ton over all size

classes of material, and are surely greater for biomass-sized

material. 

The average net revenue per acre under scenario 1 

was $493, though on some acres, the net revenue is much

greater, and on others, it is negative, sometimes considerably so. For some plots, if none of the biomass-sized material was hauled to processing sites but was instead left at

the landing, the net revenue per acre would be greater, but

in all likelihood, in situ disposal costs for this material, for

example by burning, would likely drive net revenue to zero

or below. Leaving the biomass-sized material in the woods

would reduce harvest costs and increase net revenue, but

make the estimated fire risk reduction benefits more

dubious.

While it is possible to focus on acres where net revenue

is positive, for example, under scenario 1A, this results in

a 51 percent reduction in the area treated relative to the area

that could be treated effectively, and calls into question

whether the fuels reduction achieved would be effective at

a landscape scale. Furthermore, some caveats must be taken

into account in interpreting this analysis that could even

reduce the estimate of 2.2 million acres that could be treated effectively and generate positive net revenue. At almost

every decision point requiring an assumption, we assumed

maximum land availability, minimum cost and maximum

revenue. For example, we assumed that all non-reserved

acres would be allocated to a fuels-reduction treatment. In

fact, some private landowners might well choose other prescriptions that generate more net revenue and achieve less

or no fuel reduction benefit, while others would be unwilling to consider any treatment. And, much of the National

Forest land in the study area may be in late successional

reserves, riparian reserves, or other designated uses that

could be considered incompatible with the fuel reduction

treatments we’ve considered. Nor did we include the costs

of planning, administration, site clean-up, environmental

assessments, or litigation. In estimating accessibility, we

assumed that any road on the map actually exists on the

ground and would be passable by the equipment needed to

transport equipment to the site and material from the site to

a processing facility. In haul cost calculations, the nature 

of the grid-based cost accumulation carries the implicit

assumption that any 250 m cell containing a road contains

a road heading in the direction of the processing site. Furthermore, the price used for biomass-sized material derives from

a competitive market price paid by a large-capacity biomass-fired power plant; yet in this analysis, we dispersed

biomass conversion facilities to assure a source of supply

sufficient to provide feedstocks for a number of years, but

this geographic separation also makes them unlikely to

compete with one another. So the cost and value calculations and the area availability calculations embedded in this

analysis represent best-case conditions and are optimistic.

We also assumed that all treatable lands were treated immediately; however, this would be unrealistic and even if

spread out over a few years, would still generate wood at 

a rate so much greater than current harvest patterns, that

wood product prices would likely be substantially depressed,

reducing net revenue even below the estimates reported

here.

Other caveats are that we did not consider existing processing capacity, either for biomass or merchantable material, so our results are not a justification for adding new

capacity, particularly in places like Shasta County that

203
already have substantial capacity. Nor did we factor in the

capital costs of building any processing capacity. Adding

biomass processing facilities at additional sites (beyond the

four modeled in this analysis) would reduce the haul costs

for many acres, thereby increasing net revenue. However,

the plants would need to have smaller capacity, and economies of scale, and the assumed biomass price, could be

compromised. 

Key strengths of the BioSum analytic framework are: 

1) it is based on a statistically representative sample of the

entire landscape, 2) it is connected to detailed ground observations that permit an extremely high level of resolution in

characterizing outcomes (for example, knowing the species

and size distribution of the material recovered, and being

able to model the improvement in crown fire potential), 3)

the raster-based haul cost module enables a level of spatial

detail not previously attempted without imposing enormous

costs for first building an error-free vector based road network, and 4) the approach is adaptable to any place in the

United States, by virtue of its reliance on FIA data (available nationwide as part of a National Program), publicly

available models (for example, FVS and STHARVEST),

and even topologically inferior GIS road coverages. 

ACKNOWLEDGEMENTS

We gratefully acknowledge helpful reviews of this paper

provided by Kent Connaughton, Ted Bilek, and Dave

Azuma; the data collection efforts of the PNW-FIA data

collection team, John Szymoniak and John Gerritsma for

guidance on prescriptions, Bruce Hartsough and Glenn

Murphy for assistance with STHARVEST and the harvest

cost assumptions, Nick Crookston and Don Vandendriesche

for assistance with FVS, and those who provided critical

funding support: Roger Condit of the Western Forest Leadership Coalition, David Cleaves of the National Fire Plan,

and the USDA Forest Service PNW Research Station’s

FIA, Human and Natural Resource Interactions, and

Focused Science Delivery Programs.

LITERATURE CITED

Fight, Roger D.; Zhang, Xiaoshan; Hartsough, Bruce R. 

2003. Users guide for STHARVEST: software to estimate the cost of harvesting small timber. Gen. Tech.

Rep. PNW-GTR-582. Portland, OR: U.S. Department 

of Agriculture, Forest Service, Pacific Northwest

Research Station. 12 p.

Reinhardt, Elizabeth; Crookston, Nicholas L. (Technical

Editors). 2003. The Fire and Fuels Extension to the

Forest Vegetation Simulator. Gen. Tech. Rep. RMRSGTR-116. Ogden, UT: U.S. Department of Agriculture,

Forest Service, Rocky Mountain Research Station. 209 p.

Stage, A.R. 1973. Prognosis Model for stand development. 

Res. Pap. INT-137. Ogden, UT: U.S. Department of

Agriculture, Forest Service, Intermountain Forest and

Range Experiment Station. 14 p. 

USDA, USDI (U.S. Department of Agriculture Forest 

Service, U.S. Department of the Interior). 2002. A

collaborative approach for reducing wildland fire risks

to communities and the environment, 10-year comprehensive strategy. US Department of Agriculture Forest

Service, US Department of Interior, Washington DC. 

21 p. Online: http://www.fireplan.gov/reports/11-23-

en.pdf.

Wykoff, W.R.; Crookston, N.L.; Stage, A.R. 1982. User’s 

guide to the stand Prognosis Model. Gen. Tech. Rep.

INT-133. Ogden, UT: U.S. Department of Agriculture,

Forest Service, Intermountain Forest and Range

Experiment Station. 112 p.

204 Continue