Wilderness shapes contemporary fire size distributions across landscapes of the western United States

In many U.S. federally designated wilderness areas, wildfires are likely to burn of their own accord due to favorable management policies and remote location. Previous research suggested that limitations on fire size can result from the evolution of natural fire regimes, specifically in places where fuels were recently reduced by previous burning. To explore the broader-scale importance of fire management on wilderness landscapes, we selected three study regions representing diverse ecosystems in the western U.S. and modeled the change in fire size distributions across a gradient defined by wilderness/ non-wilderness boundaries. For randomly selected locations across the gradient, we derived a scaling parameter (a) using fire size-frequency data for public lands (1984–2007); the parameter reflected the magnitude of change in the right tail of the fire size distribution where the largest fires reside. We then used quantile regression to model changes in a across the wilderness gradient, interpreting the results in terms of constraints on the relative role of large fires in structuring the fire size distribution. In the Southwest study region, the influence of large fires on size distributions decreased across the gradient toward wilderness at some places, suggesting that increased occurrence of natural burning, favored by wilderness management, led to limitations on fire sizes within recent timeframes. In contrast, we were unable to support the expectation that wilderness fire management limits the role of large fires in the Sierra Nevada and Northern Rockies study regions. Rather, the predominance of large fires increased toward wilderness interiors. Among spatial climate and topographic roughness variables included in our study, only winter and fire season precipitation limited fire size in the Northern Rockies, whereas several constraints on large fire occurrence operated in other regions. In southwestern ecosystems, evidence is needed to document stability in fire size distributions through time. In ecosystems of the Sierra Nevada and Northern Rockies, a longer temporal extent of observations may better match scales of disturbance and recovery. Our findings reflect the role of wilderness in addressing a fire deficit which has resulted from strong human influences on forests and fires over the past 150 years.


INTRODUCTION
It is likely that no place on earth is entirely free of human influence, but places set aside as wilderness have the potential to represent the far end of the continuum toward naturalness (Wiens 2007). From this vantage, the dynamics of ecosystems can evolve relatively unhindered by human interference (Wright 1974, Baker 1993. Through time, disturbances shape wilderness landscapes as part of the natural state of continual change (Pickett andOstfeld 1995, Gillson andWillis 2004).
Many federally designated wilderness areas, especially in the western U.S., include ecosystems where recurring disturbance is considered critical to defining and maintaining wilderness qualities (Cole and Landres 1996). To that end, management policies within many wilderness areas have encouraged so-called ''natural burning'' whereby lightning-ignited fires are allowed to burn with little human interference (Parsons et al. 1986). The likelihood that fires will burn of their own accord in wilderness is favored by policy; in addition, prospects for natural burning increase at times when suppression in non-wilderness areas takes higher priority, or when accessing remote areas prevents timely and effective suppression. However, natural burning policies are not applied uniformly across wilderness areas. Rather, the need to protect property and respect management policies on adjacent lands plays an important role in deciding when and where fires are allowed to burn (Kelson andLilieholm 1999, Karieva et al. 2007).
Taking advantage of opportunities to use fire in remote tracts of wilderness has been encouraged to restore and maintain ecosystems (Allen et al. 2002, Miller 2007, and a growing body of work is focused on understanding the dynamics of wilderness ecosystems under these policies (e.g., Rollins et al. 2001, Collins et al. 2009, Collins et al. 2011, Nesmith et al. 2011. Recent studies have focused on the idea that landscape self-organization evolves as feedbacks occur between burned areas and future fire events (Peterson 2002, Moritz et al. 2011. Self-organization of landscapes occurs when an ecological process, such as fire, becomes shaped by its past modifications of a landscape and thus influences ecosystem dynamics (Peterson 2002).
A case in point, Collins et al. (2009) examined fire extent and effects in mixed conifer forests of the designated wilderness of Yosemite National Park, where natural burning had been practiced for decades. Their findings suggested that fire size was limited in places where fuels were recently reduced by previous burning. However, the applicability of these findings across western ecosystems is unknown, especially given the interaction of multiple factors that can influence fire size.
Although fire sizes can be limited by previous fires' effects on fuel amount and continuity, other factors can influence the spread of individual fires and the aggregate distributions of fire sizes that result (Johnson and Miyanishi 2001, Rollins et al. 2002, Stephens 2005, Littell et al. 2009). Climate trends, especially drought and high temperatures, were key drivers of area burned in mountain ecosystems during the 20th century (Littell et al. 2009). Weather conditions, including winds and changes in temperature also affect fire spread. Furthermore, mountainous topography forms a relatively static influence on fire spread as it presents varied conditions that retard or facilitate combustion (Moritz et al. 2011). For example, linear ridges or canyons can provide a barrier to the spread of fire while steep slopes contribute to pre-heating of fuels which facilitates burning. In both cases, terrain complexity interacts with wind and fuel availability to determine fire spread behaviors. These interactions shape the dynamics of disturbance and recovery and thus affect fire size.
Here, we explored the broader-scale importance of fire management in structuring fire size distributions on wilderness landscapes through feedbacks on the burn mosaic. We asked: Have natural burning policies limited the occurrence of unexpectedly large fires by promoting fires of smaller size? We looked for evidence of this limitation in recent fire history data for three regions of the western U.S.: the Southwest, Sierra Nevada, and Northern Rockies (Fig. 1). The strength of our design lies in the potential to interpret our findings relative to the unique climatic, topographic, and ecosystem characteristics that constitute key influences on fire size in each region.
We posited that fire size distributions in wilderness interiors, or core areas, had the greatest likelihood of exhibiting changes due to management and based our study design on a continuous surface, or wilderness gradient, which we scaled to reflect that expectation. For each region, we examined the influence of topography and fuels on shaping fire size distributions across the wilderness gradient. However, we thought it unlikely that topogra-phy, fuels, or wilderness management practices would consistently influence fire size distributions at all locations. Instead, we expected one or more factors would be active constraints in certain situations represented in our sample. For example, variations in fuels and topography can affect fire spread when temperature and wind are less influential, but during times when burning conditions are extreme, landscape con- v www.esajournals.org nectivity increases and patterns in topography and fuels become less effective barriers to fire spread Romme 1994, Miller andUrban 2000).
We also presupposed that the effects of management history on fire sizes along the wilderness gradient would be evident in certain places where multiple fires were limited in size due to reduced fuels in older, adjacent burned areas (Finney et al. 2007, Rhodes andBaker 2008). These situations are more likely to occur in systems where fuel treatments have been shown to be most effective in limiting fire size and/or severe effects (Schoennagel and Nelson 2010), for example, mixed-conifer ecosystems of the Sierra Nevada (Keifer et al. 2006). Other systems where fire spread is limited more by fuel amount than moisture (e.g., seasonally dry mixed conifer and some ponderosa pine (Pinus ponderosa) forests in the Southwest) or systems where both fuel amount and fuel moisture limit the spread of fire (e.g., dry ponderosa pine forests in the Southwest; Finney et al. 2005, Meyn et al. 2007 would also be likely to exhibit self-organization on the burn mosaic. Higher elevation ecosystems which comprise a large proportion of Northern Rockies landscapes are prone to large, infrequent stand-replacing fire driven by extreme drought (Rollins et al. 2002, Schoennagel et al. 2004). However, generalizations regarding area burned in fuel versus moisture limited systems have been questioned by recent studies (Littell et al. 2009). In addition, gradients in elevation and associated variability in ecosystems allow the possibility of observing landscape self-organization in any of the regions.

Study regions
The three study regions are located in the western United States: Southwest, Sierra Nevada, and Northern Rockies (Fig. 1). The three regions are notable because each contains large wilderness areas that have famously implemented policies to encourage burning of naturally ignited fires over the past several decades. Specifically, the Gila and Aldo Leopold (Southwest), Yosemite and Sequoia-Kings Canyon (Sierra Nevada), and the Selway-Bitterroot and Frank Church-River of No Return (Northern Rockies) Wilderness Areas have followed practices described as prescribed natural fire or wildland fire use since the 1970s (van Wagtendonk 2007).
Climate.-Fire regimes in the three regions fall under different climatic influences. In the Southwest, years with more area burned occur during dry years associated with La Niña events which followed wet years associated with El Niño events (Swetnam and Baisan 2003). The Sierra Nevada has a Mediterranean climate with welldefined dry and wet seasons; drought years are coincident with more area burned, with fö hn winds promoting extreme fire behavior in the late summer and autumn (Minnich 2006). Regional-fire years in the Northern Rockies are associated with broad scale synoptic climate (Baker 2003. Infrequent high-pressure blocking systems promote extremely dry regional climate patterns Despain 1989, Bessie andJohnson 1995) and may respond to global climate anomalies (Baker 2003). Climate regimes in the three regions interact through teleconnections (Baker 2003, Macias Fauria et al. 2011; for example, years with more area burned in the Southwest tend to be quieter years in the Northern Rockies, and fire synchrony has been documented in fire scars from the Sierra Nevada and the Southwest (Swetnam and Baisan 2003).
Topography.-The diversity of landforms that occurs across the regions is attributable to their specific origins, geology, and modifications by on-going physical processes (U.S. Geological Survey 2004). In the Southwest, linear valleys, deep canyons, high mesas, mountain ranges, and rims reflect extreme stretching, uplifting, and faulting through time. The Sierra Nevada is a comprised of mostly granitic rocks which form rugged peaks; this fairly continuous North-South chain is tilted steeply on the eastern escarpment, with gradual western slopes. The Rockies are a discontinuous series of broad, high mountain ranges with distinct origins. In some places, the mountains exhibit no major crest, in others, ridges are broken with rows of granite peaks adjacent to steep canyons. Mountain ranges run primarily North-South in all three regions; the Continental Divide runs through the Northern Rockies and the Southwest, and the Great Basin Divide is located on the eastern edge of the Sierra Nevada.
v www.esajournals.org Ecosystems.-Ponderosa pine forests are widespread on the broad, flat mesa tops of the Southwest, however, wilderness areas have diverse composition (Rollins et al. 2002). The deep canyons and steep forested ridges of mountain ranges in the Aldo Leopold and Gila Wilderness Areas are inhabited by piñon-juniper (P. edulis-Juniperus spp.) and oak (Quercus spp.) at lower elevations, ponderosa pine at middle elevations transitioning to mixed forests of Douglas-fir (Pseudotsuga menziesii ), subalpine fir (Abies lasiocarpa), white fire (A. concolor), Englemann spruce (Picea engelmannii ) and aspen (Populus tremuloides) at higher elevations. The region also includes the mesas and hills of the Mazatzal Wilderness, composed of ponderosa pine-oak forest, piñon-juniper woodland and upland scrublands (Wilderness.net 2012). The Southwest region's western edge abuts the Sonoran Desert.
A gradient in forest composition also occurs from the gradual western slopes across the steep eastern escarpment in the Sierra Nevada. Foothill grassland and chaparral occur at lower elevations, while ponderosa pine woodlands and mixed-conifer forests are found at progressively higher elevations (Miller andUrban 1999, Collins et al. 2009). For example, in Sequoia-Kings Canyon Wilderness, Jeffrey pine (P. jeffreyi ), red fir (A. magnifica) and lodgepole pine (P. contorta) occur in upper-montane forests; lodgepole and western white pine (P. monticola) are more common in subalpine forests. The Mojave Desert forms Sierra Nevada region's southeastern edge.
In the Northern Rockies, wilderness areas including Frank Church-River of No Return and Selway-Bitterroot exhibit a gradient in forest composition that changes with elevation. Higher elevation forest species include whitebark pine (P. albicaulis), lodgepole pine, subalpine fir, and Englemann spruce; ponderosa pine, Douglas-fir, and western larch (Larix occidentalis) are present at middle elevations (Rollins et al. 2002. In wetter climates, including areas of the Cabinet Mountains Wilderness, western redcedar (Thuja plicata) and western hemlock (Tsuga heterophylla) are found in valleys with heathshrublands inhabiting open ridges (Wilderness. net 2012). Upmost elevations in all three regions are known for their sparsely vegetated rocky and often icy conditions. Wilderness size, configuration, and context.-Southwestern wilderness areas are generally smaller and more scattered than in the other regions (Wilderness.net 2012). Pine Mountain (8,118 ha) represents the small end of the spectrum, and the Gila is the largest (225,820 ha); it is separated from the Aldo Leopold by only a single Forest Service road. The Mazatzal and Aldo Leopold rank intermediate in size relative to other southwestern wilderness areas (102,138 ha and 81,752 ha respectively). Most wilderness areas in the Sierra Nevada are fairly large, including Sequoia-Kings Canyon (310,888 ha) and John Muir (263,851 ha) and form a fairly continuous North-South chain. Some Northern Rockies wilderness areas are even larger than those in the Sierra, for example, Frank Church-River of No Return (957,820 ha) and Selway-Bitterroot (542,516 ha); these two areas are adjacent to each other, separated by the Magruder Corridor and Highway 12. Some relatively smaller areas are scattered on the surrounding landscape, for example, Welcome Creek (11,385 ha).
Most wilderness areas are surrounded by publicly owned lands managed for multiple uses including livestock grazing, timber harvest, and recreation (U.S. Geological Survey 2012). The Sierra Nevada is distinguished by its proximity to major agricultural areas in California's Central Valley, but agricultural land use is widespread in all three regions. None of the regions contain major population centers, but some larger towns are present, e.g., Mesa, Arizona (2010 population 439,041: U.S. Census Bureau 2012), which is near the Mazatzal. Otherwise, scattered and sparsely populated development occurs in all regions.

Fire size data
We developed a fire size database that included point locations and area burned of small fires (!0.40 and 400 ha) from federal agency reports (U.S. Department of Interior 2012) and polygon centroid coordinates and polygon area for all fires in the Monitoring Trends in Burn Severity (MTBS) fire perimeter database (Eidenshink et al. 2007, Monitoring Trends in Burn Severity 2012). The MTBS maps include fires of approximately 400 ha or greater. In the process of combining the two datasets, we eliminated duplicate records in the federal report database. The location of the v www.esajournals.org fires in federal agency reports can be either polygon centroids or ignition location, and size information may be estimated in various ways. Using the MTBS data, we improved consistency in the location and size data for larger fires because these quantities were always based on digitized perimeters.
The time span of the final fire size dataset was 1984 through 2007. Although fire atlas data for earlier time periods are available for some regions, older data vary in accuracy (e.g., Rollins et al. 2002). In addition, longer records tend to underestimate the occurrence of small fires which can hamper the ability to interpret analysis of fire size distributions. We deemed the federal agency datasets and MTBS data more suitable for our purpose due to consistent standards for accuracy, despite being somewhat limited in time span. We used ArcGIS and ArcInfo software (ESRI 2010) for spatial analyses of fire size data and derivation of wilderness gradient and climate season grids described, below. All other analyses were done using the statistical software R (R Development Core Team 2011).

Sampling design
We sampled fire size data and other variables across a wilderness gradient (Fig. 2) which was constructed using the following steps. First, we derived a gridded map of wilderness/non-wilderness (assigned values of 10/0 respectively) for each study region, based on wilderness boundaries (U.S. Geological Survey 2012). Next, we calculated a uniformly weighted sum of values within ;130,000 ha moving windows across each region. The analysis window was small enough to fit entirely within many of the wilderness areas and big enough to include an adequate sample of fires for distribution modeling in most places. Thus, some of the samples represented size distributions characteristic of wilderness regimes and others characterized fire sizes across landscapes that varied in wilderness context. Variations in the shape and size of the wilderness areas resulted in a continuous gradient of values (0 , x , 10; Fig. 2). The wilderness gradient achieved its maximum value (10) in interior portions of large wilderness areas, and the minimum value (0) in areas farthest outside wilderness boundaries.
A sample of locations was selected by sorting all the wilderness gradient values for each region, randomizing the data within 20 equal-interval bins, and keeping the first 15 records from each bin (except for the Southwest, where the first 6 records were retained). In this way, samples were uniformly distributed across the gradient. The wilderness gradient surface matched the resolution and projection of the climate and elevation data described below (Albers Equal Area Projection, North American Datum of 1983, 830-m projected cell size).

Fuel condition/climate models
We modeled seasonal climatic variation as a proxy for effects of climate on fuel condition, and for effects on fuel production at the extreme ends of the spectrum in temperature and precipitation (Hannah et al. 2002, Krawchuk et al. 2009, Ackerly et al. 2010. Representing the moisture and flammability of fuels using spatial climate models provides gradients in temperature and precipitation that preserve continuous information. This is an advantage over currently available spatial fuel models that are categorical (e.g., Fuel Characteristic Classification System: U.S. Forest Service 2012). Another advantage is that it enables modeling the influence of climate on fuel condition during regionally specific time periods (i.e., before and during fire seasons).
We used the Parameter-Elevation Regressions on Independent Slopes Models (PRISM) climate data (PRISM Group at Oregon State University 2011) to represent variability in precipitation and temperature across the study regions. The PRISM models have made significant progress in explaining the extreme and complex variations in climate that occur in mountainous regions where western wilderness areas are often located (Daly et al. 2002). We used the publicly available monthly climate normals  for precipitation and maximum and minimum temperature which was the closest match to the temporal extent of the fire data.
We developed climate variables with the objective of capturing temperature extremes and amount of precipitation that tend to occur within regionally defined seasons of winter, spring, fire season, fall, and monsoon (Table 1). For precipitation, we summed the mean monthly precipitation for months of interest. For v www.esajournals.org minimum and maximum temperature, we averaged the means of those months that defined a season. These calculations resulted in multiple seasonal climate grids for each region. Each grid was summarized within the random analysis windows by applying a uniformly distributed weights matrix and then summing the values within a sample window. To define seasons that could affect fuel condi-tion, we relied on a synthesis of literature, a review of timing of fire seasons in recent years (e.g., National Interagency Fire Center 2012), and our own familiarity with the study regions. We included several combinations of the seasonal climate variables in quantile regression models described, below. Random locations were chosen to achieve a uniformly distributed sample across the wilderness gradient. At each location, fire sizes, climate seasons, and topographic roughness were summarized within ;130,000 ha windows (black box illustrates scale). Values for the wilderness gradient ranged from 10 (analysis window entirely within wilderness) to 0 (analysis window entirely outside of wilderness).

Topographic roughness models
A digital elevation model for the conterminous U.S. was obtained from the PRISM Climate Group (PRISM Group at Oregon State University 2011) and subset to the boundaries of our three study regions. Topographic roughness was quantified using fractal dimension, which provides a useful measure of roughness for many continuous surface phenomena (Burrough 1981). In a recent study, Gneiting et al. (2010) investigated several approaches for calculating fractals to assess the roughness of surface data. Based on their recommendation, we used the transectincrement method in which a variance estimate is calculated for each row and column in the grid (i.e., matrix of elevation values); the fractal dimension is the median over the set of estimates based on second differences (Gneiting et al. 2010). The algorithm is implemented in the R package fractaldim (Sevcikova et al. 2011).

Data analysis
Modeling fire size distributions.-We modeled the distribution of the fire size data to identify a parameter of the distribution that would be influenced by factors which varied across the wilderness gradient. Theoretically, disturbance phenomena including fire size-frequency data are likely to follow a distribution with a heavy right tail because they are characterized by many common and few rare events (Newman 2005). In particular, the power-law distribution has been employed because the scaling factor (a) increases where small fires are predominant in structuring fire size distributions (Malamud et al. 2005). Conversely, lower values of a indicate greater influence of larger fires on the size distribution.
We examined the linear correlation between a and univariate quantiles of the fire size sample data to provide an interpretation of the response variable, a, in terms of a simpler statistic (Fig. 3). Across regions, a was most strongly correlated with fire sizes at middle to upper quantiles; specifically, a weak negative correlation between a and fire size was observed at lower quantiles, but the relationship became increasingly negative and reached an asymptote roughly between the 60th and 80th quantiles. At the highest quantiles, a was less negatively correlated with fire size; we reasoned that variability in the far right tail of the fire size distribution, as well as the downweighting of the right tail in the exponential cutoff model lessened correlation strength.
Modeling changes in fire size distributions across landscape gradients.-Quantile regression provided an apposite method to test our research hypothesis because of its ability to identify upper and lower limits imposed by variables of interest Noon 2003, Slocum et al. 2010). The method is especially useful when measured and unmeasured factors interact to produce heterogeneous responses. In these cases, estimates based on upper and/or lower intervals of quantiles tend to be more informative than estimates based on the mean (Cade et al. 2005). Importantly, quantile regression allows a more complete picture of relationships between variables when a single rate of change based on central tendency may not adequately characterize patterns in the data Noon 2003, Cade et al. 2005).
Quantile regression models have the same interpretation as ordinary least squares models, but the models can be developed for specific portions of the data distribution. For any given quantile (s), slope coefficients represent rates of change after adjusting for the effects of the other variables in the model. At upper and lower quantiles, results can be interpreted as constraints because a majority of the sample data occurs below or above the modeled interval (s), in contrast with other regression techniques where interpretation focuses on the mean response (Cade and Noon 2003).
We estimated rates of change for five quantiles of the response distribution conditional on the matrix of predictor variables, specifically, s ¼ 0.10 and s ¼ 0.20 for lower quantile models, s ¼ 0.80 and s ¼ 0.90 for upper quantile models, and s ¼ 0.50, to compare the median response.
We used the linear quantile regression method to estimate parameters for each s and the rank score-based method to compute confidence intervals assuming independent but not identically distributed errors; these options are available in the quantreg package written in R (Koenker 2012). The scaling parameter, a was most strongly correlated with fire sizes from the 60th and 80th univariate quantile of the data distribution; uppermost quantiles were somewhat less so. The negative correlation indicates that the role of large fires in structuring size distributions decreased as a increased.
v www.esajournals.org Examining models at lower quantiles provided the advantage of modeling situations where large fires are proportionately more important (lower a) than at upper quantiles (higher a; Fig. 4). Factors which exerted lower limits on a were of particular interest because a positive constraint imposed by the wilderness gradient at the lower edge of the data distribution (s ¼ 0.10 or 0.20) would be consistent with our research hypothesis. That is, when a increases in response to a predictor variable at lower s, the relative role of large fires is lessened. Thus, a positive relationship between a and the wilderness gradient at lower s means that large fires were less important in structuring fire size distributions in wilderness, consistent with our original hypothesis. Conversely, a negative relationship at lower s means that large fires increase in relative importance as the intensity of wilderness increases, counter to our original hypothesis. On the other hand, upper limits on a operate where large fires are less predominant relative to small fires. Therefore, a positive, upper constraint is consistent with relatively more small fires.
The goal of our modeling process was to evaluate the effects of the wilderness gradient on the fire size distribution at lower, upper, and middle s, after accounting for the effects of topographic roughness and seasonal climate.
First, we formulated eight climate models for each region which represented either seasonal variability in multiple attributes or variability in one attribute (temperature or precipitation) within a given season (Table 1). A variable was considered an active constraint if the upper and lower bounds of 90% confidence intervals were either both greater than zero or both less than zero for a specified s.
Then, we used AIC c (s) to compare models with topographic roughness and climate (temperature, precipitation, or seasonal) only, with models containing the wilderness gradient as an additional predictor. The full model was: where TR ¼ topographic roughness, Climate ¼ climate variables (Table 1) and WG ¼ wilderness gradient. The AIC c (s) is a modification of Aikake's Information Criterion for small sample size (Hurvich and Tsai 1990) that incorporates differences in coefficients of determination, that is, the weighted sum of absolute deviations minimized in estimating the sth quantile regression with p parameters (Cade et al. 2005). We compared differences between AIC c (s) for more complex models and the simplest model with a constant only (b 0 ) to represent model improvement with addition of parameters (DAIC c (s) as Fig. 4. Plots of regression quantile model fits for the Southwest (left), Sierra Nevada (center), and Northern Rockies (right) study regions, where y ¼ a and x ¼ Wilderness Gradient. The black lines represent regression quantile fits for s ¼ 0.10, 0.20, 0.80, 0.90; the blue line is the median fit (s ¼ 0.50), and data points are shown as black points. In the Southwest, patterns suggest that interacting factors across the wilderness gradient resulted in opposite trends at upper (À) and lower (þ) quantiles. Patterns in the Sierra Nevada data indicate additive effects (i.e., roughly parallel lines) but rates of change are less dramatic than the strong, negative trends across quantiles for the Northern Rockies. Interpretations of these patterns were based on simulations done by Cade et al. (2005).
v www.esajournals.org defined in Cade et al. 2005). The resulting scale of DAIC c (s) provides a positive measure that is proportional to model fit, as well as a measure of relative improvement as variables are added to simpler models. We calculated DAIC c (s) across all s (0.05 to 0.95 in increments of 0.05) to observe changes in models representing lower, middle, and upper constraints on a.

Southwest
Our findings suggest that increased occurrence of unhindered burning, favored by wilderness management practices, led to interactions among fires that limited fire sizes on the Southwest burn mosaic within recent timeframes. In support of our hypothesis, we identified locations where the wilderness gradient actively constrained the role of large fires in structuring fire size distributions. Interacting factors across the wilderness gradient resulted in positive trends at lower quantiles (Fig.  4). Specifically, the wilderness gradient played a role in limiting the lower range of a after considering topographic roughness and climate in two models: Minimum Temperature and Spring ( Fig. 5; Southwest: Positive Limits). In these two models, the wilderness gradient imposed a significant positive constraint on a at lower s, indicating a reduced role of large fires in structuring the distribution of fire sizes in locations where larger fires are more prevalent (i.e., places with the lowest a).
Notably, where large fires played a greater role in structuring size distributions (i.e., lower s), all significant variables worked to limit fire size in the Southwest (i.e., increased a; Fig. 5; Southwest: Positive Limits). All models except Maximum Temperature (Table 1) included variables that imposed positive limits on a at lower s ( Fig.  5; Southwest: Positive Limits). Significant predictors included topographic roughness, minimum temperature and precipitation during the Winter and Fire Season and minimum temperature during Monsoon.
Where large fires played a lesser role (i.e., upper s), our findings suggested that wilderness management was one of several factors which contributed to a relative decrease in small fires (i.e., proportionately more large fires). At middle and upper s, the wilderness gradient imposed a negative constraint on a in all models, indicating places where a was restricted to smaller values (i.e., the role of small fires decreased; Fig. 5; Southwest: Negative Limits). Upper and lower constraints resulted in wilderness landscapes exhibiting a restricted range of a in which large fires have an intermediate influence on size distributions.
Although the positive, lower limit imposed by the wilderness gradient was statistically significant in Minimum Temperature and Spring Models (Fig. 5), it added relatively little information to models representing lower constraints (Fig. 6). Greater improvement resulted from addition of the wilderness gradient at middle and upper s. Models tended to be better overall, in terms of fit and explanatory power, at the lowest values of s, then DAIC c (s) decreased to a minimum where s ¼ 0.2 to 0.3 before steadily and, in some cases, rapidly increasing to the best fit at the highest s (Fig. 6).

Sierra Nevada
Our analysis of the Sierra Nevada study region produced no evidence to support the expectation that wilderness management limits the overall proportion of large fires. Visual trends in univariate models ( Fig. 4; gradual or no change in slope) and significance tests in multiple regression models ( Fig. 5; Sierra Nevada: Positive and Negative Limits without ''WG'') indicated that the influence of large fires on size distributions did not depend on the wilderness gradient in many cases. In places where large fires had greatest influence (smaller a and lower s), several variables worked to favor smaller fires in the Sierra Nevada ( Fig. 5; Sierra Nevada: Positive Limits); the wilderness gradient and other variables had the opposite effect by increasing the predominance of large fires ( Fig. 5; Sierra Nevada: Negative Limits).
A diverse set of variables were counted as positive lower limits in the Sierra Nevada ( Fig. 5; Sierra Nevada: Positive Limits); specifically, topographic roughness in the Spring model, and seasonal factors including precipitation in the Fire Season Model, and Winter in the Precipitation Model. In addition, a was less influenced by large fires in places with higher maximum temperature in the Winter Model, and in the Maximum Temperature Model. Fig. 5. Stacked bar charts summarizing the quantile multiple regression results. The height of each bar represents the number of significant variables for a given model at lower (gray), middle (light blue), and upper (dark blue) s. Variables that imposed positive limits on a (e.g., less large/more small fires at lower s) are shown in the left column; negative limits (e.g., less small/more large fires at lower s) are counted in the right column. Model names (x-axis labels) indicate the combination of climate variables tested in multiple regression models (Table 1). Bar segments for models in which the wilderness gradient was one of the counted (significant) variables are marked ''WG.'' Interpretations and discussion of specific variables can be found in the text.
v www.esajournals.org The wilderness gradient was generally a worthwhile addition to simpler models based on climate and topography across all s (Fig. 7). The best models usually occurred at uppermost s, but in some cases (e.g., Monsoon and Fire Season) models at lowest s were comparable. The Precipitation model had the highest DAIC c (s) and the wilderness gradient added greatest benefit to model fit as an upper constraint on a.  (Table 1), before and after adding the wilderness gradient (WG). The wilderness gradient improved models at middle and upper s, but added little to models representing lower constraints.

Northern Rockies
In the Northern Rockies, conditions favorable to greater influence of large fire events increased dramatically toward wilderness interiors (Fig.4); these patterns persisted after taking into account the role of topographic roughness and fuel condition (i.e., spatial climate). Several variables, including topographic roughness, worked to-gether with the wilderness gradient to increase the role of large fires (Fig. 5;Northern Rockies: Negative Limits). Constraint of a to smaller values (i.e., more large/less small fires) occurred at lower, upper, and middle data intervals (Fig. 5;Northern Rockies: Negative Limits).
We identified only two factors that limited the role of large fires at lower s in the Northern Rockies. In places with more precipitation (Precipitation model: W 1-3 and FS 6-9; Table  1), large fires were less predominant (i.e., a was restricted to higher values; Fig. 5; Northern Rockies: Positive Limits). These factors continued to impose limits on a at middle and upper s, indicating a consistent relationship between higher precipitation during winter and fire season and proportionately more small fires. The additive effects suggested by parallel quantile regression slopes (Fig. 4), and the scarcity of measured variables limiting the role of large fires indicated that environmental conditions act in concert to structure fire regimes across the Northern Rockies wilderness gradient.
Of the three regions, the Northern Rockies models showed the most dramatic improvement with addition of the wilderness gradient to the simpler climate and topography models, based on DAIC c (s) (Fig. 8). Patterns in the middle range of s were pronounced in this regard, reflecting a strong constraint on a by the gradient in wilderness at the median of the data distribution. Models at middle s were better fit than those at upper or lower parts of the data distribution. A significant, negative limit was imposed by the wilderness gradient in all models (Fig. 5;Northern Rockies: Negative Limits). Across s, the role of large fires increased and the role of small fires decreased (lower a) along the nonwilderness-wilderness continuum.

DISCUSSION
The current belief that natural burning fostered in wilderness promotes self-organization on burn mosaics has been subjected to limited tests, but still guides management practice today. In our empirical study, we were able to support this assumption only in the Southwest, with contradictory results in the other two regions. Ultimately, detecting effects of natural burning practices requires sufficient instances where fires result in landscape heterogeneity that affects the spread of subsequent fires (Turner et al. 1989). In southwestern ecosystems, the short-term synchronicity of fire with (2-6 yr) fluctuations in La Niña/El Niño (Swetnam and Baisan 2003) likely provides an environment conducive to rapid evolution of organization in landscape structure. In places where wilderness management was implicated in reducing the role of large fires, our findings suggest that the cumulative effect of interactions between disturbance events has resulted in a balance between disturbance and recovery, as patchy disturbance corresponds to predictable regrowth (Sprugel 1991).
Due to the composition and arrangement of ecosystems in the Sierra Nevada, longer timeframes may be required for the maintenance of structure in burn mosaics to become apparent at a regional scale. The region includes a diverse range of vegetation spanning gradients in elevation and moisture along which fire regimes vary (Miller and Urban 2000). Indeed, observations that suggested self-organization in Yosemite included a longer time span (1971Collins et al. 2009) than the data used in our study . We cannot rule out the possibility that with a natural fire program operating in wilderness over a longer time period, trends observed in the Yosemite study will become apparent at broader scales.
Alternatively, limitations of our study may explain the lack of evidence in support of our original expectations in the Sierra Nevada. When posing our hypothesis, we implicitly used the wilderness gradient as a proxy for a fire management gradient, with more fires burning unimpeded in the interior of wilderness. Our difficulty associating observed constraints on fire size with wilderness management may be due to the inconsistent implementation and duration of natural fire programs. Because management of wildfires, even in wilderness, is fraught with uncertainties, fire suppression is deeply seated in institutional principles (Stephens andRuth 2005, Miller et al. 2011). Decisions to suppress fires in wilderness areas can have dramatic effects on landscape conditions (Davis and Miller 2010). Furthermore, management and land use of adjacent areas often hinders achievement of wilderness objectives (Kelson and Lilieholm 1999). Consequently, it is possible that wilderness is a poor proxy for natural burning over a sufficient period to affect self-organization, and it may be poorer in some regions (e.g., Sierra Nevada) than others (e.g., Southwest).
In addition, the wilderness gradient is potentially confounded with other landscape attributes that can influence the fire regime. In particular, outside of wilderness, anthropogenic landscape fragmentation caused by roads and other human structures, as well as land use practices such as agriculture and livestock grazing, can disrupt the continuity of fuels and impede the spread of fires Moritz 2009, Shinneman et al. 2010). In contrast, fragmentation in wilderness areas is more often defined by natural boundaries that modulate fire spread including rock, ice, and topographic variability. The structure of landscape ''fences and corridors'' (sensu Moritz et al. 2011) differs across the wilderness gradient, and potentially contributes to quantifiably different fire regimes. These landscape attributes may be counteracting the propensity for natural burning to decrease the relative role of large fires, and more so in some regions than others. The limitations of our study design mentioned in relation to the Sierra Nevada may also apply to the Northern Rockies, but the dramatic influence of the wilderness gradient on predominance of large fires compelled us to an alternative interpretation of our findings. First, fire-adapted landscapes in the Northern Rockies study region may be less likely to exhibit properties suggesting self-organization within the timeframe of our study due to the configuration of ecosystems and fire-climate regimes. Opportunities for feedbacks among fire events would be limited to lower-and mid-elevation forests where fire frequency and severity can be more variable (Fulé et al. 2003, Sherriff andVeblen 2006); the predominance of subalpine ecosystems (Rollins et al. 2001) makes it more likely that fires are infrequent, and typically large . Second, there is no consistent relationship between time elapsed since the last fire and fuel abundance in subalpine forests (Brown and Bevins 1986), making the results of interactions between fires more difficult to predict. Finally, the timeframe of our study coincides with an increased occurrence of large fires in the Rocky Mountains associated with broad-scale climate change (Westerling et al. 2006); our findings suggest that the manifestation of recent regional fire years is more apparent in wilderness core areas, where management practices and environmental conditions are also conducive to extensive burning.
From a long-term perspective, recent large fires in the western states occurred during a fire deficit which has resulted from strong human influences on forests and fires over the past 150 years (Marlon et al. 2012). If wilderness areas were less affected by human activities, both historically due to remoteness and more recently with the addition of management policy, wilderness fire regimes may exhibit closer correspondence to climate dynamics. Given conditions that are warmer and drier than those that occurred previously (Westerling et al. 2006, Marlon et al. 2012, the predominance of large fires in wilderness suggests that allowance of natural burning may effectively promote a coupling of fire and climate that would likely be unwisely obstructed by fire exclusion and suppression in the long run.

CONCLUSION
Through the framework provided by a gradient across wilderness landscapes, we gained insights into regional differences in fire regimes and how they are structured by spatial patterns of fuels and topography. In the Southwest, evidence of a stable fire size distribution over longer timeframes would give stronger credence to the idea that wilderness management can lead to systems in which natural disturbance processes maintain structure and function through internal feedback mechanisms. If the self-organization we observed in southwestern wilderness ecosystems leads to resiliency, attributes such as community composition and biodiversity may fluctuate within a relatively constant range of variability (Holling 1973).
The management message in regions where fire size distributions reflect a greater role of large fires is unclear. Large disturbances can result in wide swings in communities and ecosystems over long timeframes (Sprugel 1991), but these changes must be considered relative to spatial and temporal scales of disturbance and recovery (Turner et al. 1993). Extending the temporal scale of observation and taking into account exogenous factors that influence resiliency to disturbance would be required to fairly interpret contemporary fire regimes. Opportunities to observe the dynamics of natural processes in wilderness are essential to these endeavors.

ACKNOWLEDGMENTS
Critical reviews of an earlier draft by Brian Cade, William Romme, and four anonymous reviewers were helpful in improving the manuscript. Thanks to Brian Cade, Ethan Plunkett, and Nick Povak for assistance with analysis and R code. Discussions with Paul Hessburg and Don McKenzie on fitting power-law distributions improved the rigor of our approach to the analysis. Stephen Howard and Vicky Johnson were extremely helpful with access to and interpretation of burn perimeter and fire history datasets. The project was funded by a Joint Venture agreement between the USDA Forest Service Rocky Mountain Research Station (Research Work Unit Number 4901) and the University of Massachusetts.