Assessing the Minimum Number of Time Since Last Fire Sample-Points Required to Estimate the Fire Cycle: Influences of Fire Rotation Length and Study Area Scale

Boreal forest fire history is typically reconstructed using tree-ring based time since last fire (TSLF) frequency distributions from across the landscape. We employed stochastic landscape fire simulations to assess how large a study area and how many TSLF sample-points are required to estimate the fire cycle (FC) within a given accuracy, and if those requirements change with length of the simulated fire rotation (FRS). FRS is calculated from simulated fire-year maps used to create the TSLF map, and is the “true” measure of fire history that FC estimates should equal. Fire-year maps were created by (i) using a spatially homogenous landscape, (ii) imposing large variations in annual area burned, and (iii) having no age-related change in the hazard of burning. We found that study areas should be ≥3× the size of largest total annual area burned, with smaller-scale areas having a bias that cannot be fixed by employing more samples. For a study area scale of 3×, a FC estimate with an error <10% was obtained with 187 TSLF samples at 0.81 samples per 100 km2. FC estimates were not biased in study area scales that were ≥3×, but smaller-scale areas with a short FRS had an overestimated FC and smaller-scale areas with a long FRS had an underestimated FC. Site specific variations in environmentaland age-related variations in the hazard of burning may require more sample-points; site specific simulations should thus be conducted to determine sample numbers before conducting a TSLF field study.


Introduction
There is concern that fire frequency is increasing globally [1] and that the resulting carbon emissions will enhance climate change [2].Fire frequency is here defined as the probability that a forest stand burns [3,4].An understanding of how fire frequency has changed, and the climatic and anthropogenic factors that have driven those changes, can help anticipate future changes in fire frequency.In areas of high fire frequency, the use of decades-long historical and remotely sensed data can indicate temporal changes in wildfire frequency [5], while areas of lower fire frequency require longer-term biophysical records such as provided by fossil charcoal and tree-rings [6].
Fire frequency is equal to the inverse of each of the fire cycle (FC), fire rotation (FR), and mean fire interval (MFI), if the hazard of burning is constant over historical time and forest age [4,7].FC and FR are both equal to the number of years required to burn an area equal in size to the study area, but they differ in the data employed for their estimation.FR is estimated from a series of annual fire-year maps created using either historical records or tree-ring reconstructions [8].FC is estimated from a map that displays the times since last fire (TSLF) across the landscape; a frequency distribution of TSLF ages typically displays a negative exponential form and its FC can be obtained using a maximum likelihood estimate (MLE) [3].MFI is the mean number of years between consecutive fires, as estimated from peaks in fossil charcoal records, fire scars on a tree, or polygons on annual fire maps created using historical records or stand origin maps [7].This paper will focus on FC estimation, but will use the simulated value of FR (FR S , see below) to evaluate the accuracy of FC estimates.
FC estimation using TSLF data typically has three steps [7].First, a study area is identified that is ≥3× larger than the largest fire year, a size believed to ensure that individual years do not over-influence the shape of the TSLF frequency distribution.Second, the TSLF at each sample-point is dated using tree-ring pith origins or fire scar data and, if available, historical fire records.While the TSLF usually represents the last stand replacing fire as indicated by tree ages, it could also represent the last sub-lethal fire as indicated by fire scars [7].Third, the TSLF data from across the landscape is statistically analyzed to estimate the FC.The third step can be conducted using point-sample data, or using a complete TSLF map.However, it has been suggested that TSLF studies should employ point-sample data as the creation of TSLF maps through interpretation of aerial photographs can result in pseudo-replication [9].Additional methods may be employed, depending on the study area.For example, for study areas with spatial and temporal heterogeneity, methods have been developed to estimate variations in FC over time [10] and space [11].In other areas, the potential problem of censoring of old TSLF dates due to tree senescence [12] has been addressed [13].
In this study, we focus on three issues related to FC estimation using TSLF data: sample number, study area size, and FR length.First, no research has suggested how many TSLF points should be employed to obtain an accurate estimate of FC.Landscape fire simulation analyses have, however, shown that increased numbers of TSLF samples result in decreased variance in FC estimates and FC estimates closer to the "true," simulated FR [4,13].Since the FR S is calculated with all simulated fire-year maps used to create the TSLF map, FR S represents the "true" value that FC estimates should equal.Second, some results [13] appear to show that as FR S lengthens, the variance in FC estimates also lengthens, which may indicate the need for more TSLF samples.Third, the suggestion that the study area should be ≥3× the largest total annual area burned [7] has not been evaluated, though studies have exhibited influences of scaling on related measures.Smaller areas tend to have a greater variance in stand ages [14], and greater variation around the expected negative exponential TSLF distribution [15,16].Related to this, the capacity to detect relations between historical fire records and environmental drivers was stronger for larger study areas [17].Two empirical studies whose studied areas were spatially nested indicate that the FC estimate from the smaller area [3] was shorter than the FR estimate from the larger area [8].
Eighteen representative field studies that estimated FC using TSLF sample-points varied greatly in study parameters: study area ranged from 78 to 44,870 km 2 (median 2250 km 2 ); ratio of study area to largest total annual area burned that could be calculated for five studies ranged from 1.7 to 16.2 (median 4.1); number of TSLF sample-points ranged from 75 to 3168 (median 197); and sample density ranged from 0.4 to 185.6 per 100 km 2 (median 10.3 per 100 km 2 ) [8,9,[18][19][20][21][22][23][24][25][26][27][28][29][30][31][32][33].Variations in study area size were often due to the size of the conservation area that was studied, while variations in sample number were likely due to operational constraints such as time and money.It is possible that if guidelines on TSLF sample number requirements had been available then those studies might have taken more TSLF samples to increase FC accuracy, or perhaps fewer TSLF samples so as to be more cost-effective.
We will use a simulation approach to assess the number of TSLF sample-points required to estimate the FC with a given accuracy in response to study area size and the FR S .Since the FR S is calculated using all simulated fire-year maps used to create the simulated TSLF map, a good FC estimate from the TSLF map should be equal to the value of FR S .Over a dozen computer models simulate landscape fires [34], many of them have been used to simulate TSLF maps, e.g., [35,36], and simulated TSLF maps have been used to assess the influence of number of TSLF samples on FC estimates [4,13].First, we develop a statistical approach to indicate the minimum number of TSLF samples required for a fire history study to have a FC estimate within a predefined error.
Second, we evaluate how the error in the FC estimate and how the number of required TSLF samples varies with the ratio of study area size to largest total annual area burned.Third, we evaluate if error in the FC estimates varies with FR S .To address these issues we simulate fire in a homogenous landscape with no long-term climatic variation, as methods to address spatial and temporal homogeneity have already been developed.

Landis-II Model
To simulate fires and develop TSLF maps, we used LANDIS-II (Landscape disturbance and succession model, version 2) [37] and its Dynamic Fire and Fuel System (DFFS) extension [36].LANDIS-II was designed to simulate broad-scale (>100 km 2 ) forest landscape dynamics [38], including succession, insect disturbance, carbon fluxes, and impacts of climate change, e.g., [39][40][41].Landscapes contain interacting cells with user defined resolutions, and individual cells have homogenous forest cover and soil conditions [40].It allows landscape conditions and dynamics to be parameterized using empirical data that reflect historical conditions [42,43].The DFFS extension allows simulated fire behavior to respond to temporal variations in fire weather [39] thus increasing realism in fire shape.

Historical Forest Fire Data
The historical fire and weather data used to parameterize LANDIS-II and DFFS were from Wood Buffalo National Park (WBNP), a 44,807 km 2 area of boreal forest in Alberta and adjacent Northwest Territories.These data were collected by the WBNP warden service from 1969 to 2012.During that 44-year period there were 1556 forest fires that ranged in size from 0.001 to 1812.5 km 2 .The largest 10% of the fires caused 99% of the total area burned, exhibiting a fire size distribution similar to other areas [44].The annual number of fires ranged from 5 to 65, with a mean of 36.The annual burned area ranged from 0.133 to 6519 km 2 , with a mean of 560 km 2 .The historical fire rotation (FR H ; see Section 2.8) for this 44-year period was thus 80.0 years (i.e., 44,807 km 2 divided by 560 km 2 per year).
The following fire weather data associated with each historical fire, were obtained for the model fire parameterization (Section 2.3): fire season (spring: March to May; summer: June to August; fall: September to November), leaf status (leaf-on: summer; leaf-off: spring and fall), wind speed velocity (km/h), wind direction (degrees), fine fuel moisture code (unitless), buildup index (unitless), and weather class (moderate to extreme, as related to fire size by LANDIS-II) [36].As 253 fire records lacked some of that information, only 1303 fire weather records in the 44 years could be used in the weather information table.

Model Fire Parameterization
In DFFS, fire sizes are assumed to follow a lognormal distribution that is described using the mean and standard deviation of the natural logarithm of fire sizes [36].Our preliminary simulations were parameterized using the mean and standard deviation of the complete 44-year historical record of fire sizes and annual numbers of fires.However, we found that those simulations predicted much less inter-annual variation in fire numbers and area burned than occurred in the historical record.To allow simulations to predict the full observed range of fire activity, we thus divided the historical record into six groups of years based on natural breaks in a rank-ordering of their total annual area burned (Figure 1).In classes with progressively smaller total annual area burned, the historical record of fire had progressive decreases in the size of the largest individual fire, mean percent of the landscape burned annually, and percent of total annual area burned (Table 1).
Figure 1.The 44-year historical record of total annual area burned, rank-ordered from most-to least-area burned in a year.Classes of total annual area burned: Class 6: <1 km 2 ; Class 5: 1-9.9 km 2 ; Class 4: 10-99.9km 2 ; Class 3: 100-399.9km 2 ; Class 2: 400-999.9km 2 ; Class 1: ≥1000 km 2. For each of those six classes of total annual area burned we calculated the mean and standard deviation of the log-transformed fire sizes and the mean number of fires per year (Table 2).The values in Table 2 differ from the WBNP records in two ways.First, the mean number of fires per year in WBNP was multiplied by 1.88 as the 2900 × 2900 grid simulated area (see Section 2.4) was 1.88 the size of WBNP.Second, the many fires <0.01 km 2 (<1 ha) were discarded so as to not over-influence the parameter values obtained for mean and standard deviation of fire sizes.Each of the six sets of fire parameters (Table 2) were used to simulate 200 individual years of fires, for a total of 1200 years of simulated fires.A size-based, rather than duration-based, fire simulation approach in DFFS was employed using four steps to increase naturalness of the spatial pattern of fires in annual fire-year maps.First, a fire was initiated and its size was randomly selected from our predefined lognormal distribution for one of the six classes of total annual area burned.Second, the fire spread to cells selected using spread equations adapted from the Canadian Fire Behavior Prediction System [45].The weather parameters for these spread equations were randomly selected from the weather table for the weather class associated with that fire size; those weather conditions were used to simulate the entire fire event.Third, the fire simulation was terminated when either the number of cells selected multiplied by the cell area equaled or exceeded the predetermined fire size, or it reached the edge of the simulated landscape.Fourth, additional fires were simulated until the total number of fires that year was equal to the number determined by a random draw from the Poisson distribution [36] for that class of total annual area burned (Table 2).

Landscape Conditions Input
LANDIS-II and DFFS require information on landscape size, topography, waterbodies, tree species coverage, and forest age.Our simulated landscape contained 2900 × 2900 cell that were each 100 m × 100 m (0.01 km 2 ).Following the simulation of a fire year, the outer 400 cells (40 km) were clipped from all four edges of the landscape to account for reduced fire occurrence due to the "probability shadow" [46].By removing that buffer area, all areas have an equal chance to be burned by peripheral fires.The remaining 2100 × 2100 cell landscape was 44,100 km 2 , similar in size to the 44,807 km 2 WBNP from where the fire parameterization data were obtained.
Each of the 1200 years of fire simulations was initiated on a homogenous 60-year-old jack pine (Pinus banksiana Lamb.) forest with no spatial variation in environmental conditions or in the probability of fire ignitions.This was done as all landscapes differ, and thus this simple landscape would serve as a baseline condition.In each of our simulations the same forest age (i.e., 60 years) was used so that multiple simulations could be run at once, and so that each year's fire map could be used in the creation of multiple TSLF maps (see Section 2.5 ).
While each of the 1200 years of fire simulations was initiated on a 60-year-old forest, the information that we extracted from that simulation was the spatial pattern of fires created in that year.From each of the 1200 simulations we thus obtained one fire-year map of the areas burned in that year, providing a total of 1200 fire-year maps.

Time Since Last Fire (TSLF) Map Creation
We simulated 200 years of fires in each of the six classes of total annual area burned, for a total of 1200 years.Since most of the individual fires in Class 6 years were smaller than the cell size of 0.01 km 2 , we simply used a "no fire" map to represent total annual area burned in those years.Thus, there were actually 1000 years of simulated fires.
TSLF maps were created by overlaying the 1200 maps of total annual area burned in the following sequence of classes of total annual area burned: 6, 5, 4, 3, 2, 1, 6, 5, 4 . . . 1 until all 1200 maps had been stacked (Figure 2).The stacking operation was undertaken using a program we wrote in Python.Each annual fire-year map was randomly selected (without replacement) from the appropriate class of total annual area burned, until all 1200 fire-year maps had been selected.The first-used fire-year map represented 1200 years ago, the second-used fire-year map represented 1199 years ago, and the last-used fire-year map represented one year ago.When cells burned in earlier years were overlaid (i.e., over-burned) by more recent fires, their age changed to that of the newer burn.

Creating TSLF Maps with Different Fire Rotations (FRs)
TSLF maps with a range of FRs were created by leveraging the independence of the 1200 fire-year maps with which they were built, and the influence that the largest fire-years have on the FR [44].We created a total of 201 TSLF maps that contained progressively fewer of the large, Class 1 fire-year maps: 200, 198,197 . . .2, 1, 0. To insure that the 201 TSLF maps were unique in their spatial pattern of fires, they were each created using a new random sequence of fire-year maps (described in Section 2.5).Once stacked, j of the Class 1 fire-year maps were removed by first dividing all fire-year maps into j intervals, and then removing one randomly chosen Class 1 fire-year map from each of the j intervals.This was done for j = 1, 2, 3 . . .200.By removing the 200 Class 1 fire-year maps in this strategic manner, we ensured that we would obtain TSLF maps that would have a wide variety of FR S .
Uniqueness of the spatial pattern of forest ages in the TSLF maps was evaluated by comparing TSLF ages in the 11 TSLF maps that contained 200, 180 . . .20, 0 of the Class 1 fire-year maps.The 11 TSLF maps had their 2100 × 2100 grid overlain with a lattice which contained 900 nodes that were each separated by 68 grid-cells.The 900 nodes were at the following XY grid-pairs: 68-68, 68-136, 68-204 . . .136-68, 136-136 . . .2040-2040.The TSLFs at those 900 nodes were extracted from each of the 11 TSLF maps.Spearman rank correlation coefficients were calculated between all pairs of the 11 node-constrained arrays of TSLFs, and the frequency distribution of those 55 paired correlations was examined.Spatial autocorrelation in the 11 TSLF maps was evaluated using the global Moran's I [47].

Creating TSLF Maps at Different Study Area Scales
To assess the influence of study area scale, the 201 TSLF maps with different FRs (see Section 2.6) were evaluated at nine different study area scales.The scales were of the study area size relative to the largest historical total annual area burned.For example, the 2100 × 2100 grid simulated landscape is 6× larger than the largest total annual area burned in the 1200 years of fire simulations (see Section 3.1).

Fire Rotation (FR) and Fire Cycle (FC) Estimation Methods
Fire rotation (FR) is calculated using annual fire maps: study area size divided by the mean total annual area burned.The historical fire rotation (FR H ) for WBNP was calculated using its 44-year historic record.The simulated fire rotation (FR S ) for the 2100 × 2100 grid simulated landscape was calculated using its 1200 simulated annual fire-year maps.A Python program was written to conduct a Monte Carlo simulation to estimate the mean FR S and its standard deviation.The Python program randomly selected without replacement 100 of the 200 annual fire maps in each of the six classes, and then used those 600 maps to calculate one FR S ; this operation was repeated 1000 times providing 1000 estimates of the FR S from which the mean and standard deviation were calculated.
The FC of the TSLF map created using all 1200 annual fire-maps was calculated using two methods: the mean age (FC MA ), and the maximum likelihood estimate (FC MLE ).The FC MA was estimated as the mean age of 1000 TSLF sample-points; the FC MA was introduced as a simple MLE of a TSLF distribution with a negative exponential form [3]. The FC MLE of 1000 TSLF points was estimated using a more elaborate maximum likelihood function [48] that accounted for two sources of randomness: the random occurrence and spread of fires, and the random selection of sample-points [49].The 1000 points used to calculate FC MA and FC MLE were distributed across the landscape of the TSLF map in a structured random approach: discretizing the landscape into 100 equal-sized non-overlapping areas of 210 × 210 grid (441 km 2 ), and then obtaining the TSLF from 10 different randomly selected, without replacement, cells from each of the 100 areas.For both methods of estimating the FC, the historical hazard of burning in the simulation would be constant, because the six classes of total annual area burned that we employed would simply cause a six-year cycle around a multi-century trend.The two methods were employed 1000 times.Differences between estimates of FR S , FC MA and FC MLE were assessed using a two-sample t-test.

Estimating the Minimum Required Number of TSLF Sample-points
To assess the minimum number of sample-points required to estimate the FC MA at a given level of accuracy, we wrote a Python program to conduct a five-step procedure.First, the TSLF map was discretized to obtain 100 square grids numbered 1 to 100.Second, 1000 sample-points were randomly selected in a repeated sequence from the grids 1 to 100; reaching 1000 points required 10 cycles through the 100 grids.Third, each time a new sample-point was added, the FC MA was recalculated, until FC MA was estimated with 1000 sample-points.FC MA was estimated as the sum of TSLF ages for the n points divided by n.Fourth, the number of sample-points required to estimate mean FC MA at a given level of accuracy was determined by the absolute difference in the measured mean (ADMM): where abs is the absolute value, Mean Age(1000) is the mean of the 1000 randomly selected sample-points, and S order is the order of the 1000 samples.Accuracy in this study was assessed using an ADMM threshold of 10%, but the program can be modified to employ any value.Fifth, the minimum sample size was identified as the lowest sample size at which the ADMM stayed lower than the threshold value of 10% for the remainder of the 1000 sample-points.The 1000 different estimates of the minimum sample size were then expressed as their 95th percentile.

Influence of Study Area Scale on Estimated FC and TSLF Age Variability
The influence of study area scale (study area size divided by largest total annual area burned) on FC MA estimation was evaluated in four steps.First, the FRs was calculated for all 201 TSLF maps for the whole 2100 × 2100 grid landscape.Since FRs was calculated using all of the individual fire-year maps that were used to create the TSLF map, FRs provides the "true" value that a good estimate of FC MA should also provide.Second, FC MA was estimated for all 201 TSLF maps at all nine scales, for a total of 1809 values.Third, at each of the nine scales, the root mean square error (RMSE) and R 2 adj were calculated between the 201 FR S values and FC MA estimates.Fourth, the RMSE and R 2 adj values were graphed across the nine scales to assess if there was a scale at which the FR S and FC MA no longer agreed well.
The influence of study area scale on age variability across TSLF maps was evaluated in four steps.First, the FRs was calculated for all 201 TSLF maps for the whole 2100 × 2100 grid landscape.Second, 1000 TSLF points were selected from each of the 1809 TSLF maps and their standard deviation was calculated; this was done 1000 times for each map and their mean was calculated.Note: in a negative exponential distribution, the standard deviation and the mean should have equal values [3].Third, at each of the nine scales, the RMSE and R 2 adj were calculated between the 201 FR S values and the standard deviation of the TSLF ages.Fourth, the RMSE and R 2 adj values were graphed across the nine scales to assess if there was a scale at which the two measures indicate poor agreement between FR S and the standard deviation.

Influence of FR S and Study Area Scale on Minimum Required Sample-Points
The influence of the FR S on the accuracy of the FC MA estimate was evaluated in three steps.First, the value of FR S minus FC MA (i.e., FR S -FC MA ) was calculated for all 201 TSLF maps at all nine scales, for a total of 1809 values.FR S was always for the whole 2100 × 2100 grid landscape, as it represents the basic conditions behind the creation of the TSLF map.Second, linear regressions were conducted between the 201 values of FRs against FRs minus FC MA , with this done separately for each of the nine scales.Third, the intercept and slope parameters and p-values were tabulated; study-area scales at which the intercept and slope were not significantly different than zero were inferred to have had no influence on the length or the FR S on the length of the FC MA estimate.
The minimum number of TSLF samples required to estimate FC MA at a given level of accuracy was assessed using four steps.First, for each of the 201 TSLF maps at each of the nine scales we determined the 95th percentile of the minimum number of TSLF sample-points required to obtain an ADMM of 10% (see Section 2.9).Second, those values were expressed in two ways: the number of TSLF points per map and, since the nine study area scales differed greatly in total area, the number of TSLF points per 100 km 2 .Third, linear regressions were conducted between the 201 values of FR S against the 95th percentile of the minimum number of sample sizes for each of the 201 maps; this was done separately at each of the nine scales.The FR S was for the whole 2100 × 2100 grid landscape.Fourth, the 95th percentile of the 201 minimum sample values was calculated for each of the nine scales and graphed.

Simulated Fires and Fire-Year Maps
The largest simulated fire in the 1200 fire-year maps was 1562 km 2 (Table 1; 3.5% of the simulated landscape).The largest of the 1200 simulated fire years burned 7385 km 2 (16.7% of the landscape), and the mean total annual area burned for the 1200 fire maps was 717.5 km 2 .The fire characteristics in the six classes of total annual area burned was similar for the 1200 years of simulated fires as it was for the 44 years of historical fires (Table 1) in terms of both the absolute size of the largest individual fire, as well as the mean percent of the landscape burned annually and the percent of the total area burned.
Rank-ordered distributions of the total annual area burned in the 1200 years of simulated fires and the 44 years of historical fire records both exhibit a negative exponential distribution (Figure 4).The distributions are more similar for large fire years than for small fire years; for example, the largest 30% of the historical fire years and the largest 34% of the simulated fire years each burned >1% of the study area per year.The historical data exhibit a steeper slope in the frequency of total annual area burned than do the simulated data.The Monte Carlo estimate of FR S from the 1200 fire-year maps was 61.5 years (range: 57.6 to 66.3 years).

Time Since Last Fire Map and FC Estimates Using Initial Fire Parameters
The TSLF map for the 2100 × 2100 grid landscape, created using 1200 fire-year maps based on the fire parameters in Table 2, contains much young forest and progressively less old forest (Figure 3a).The cumulative frequency distribution of TSLFs for all 4,410,000 cells from that TSLF map (Figure 3a) displays a negative exponential distribution (Figure 5), the slope of which indicates a FC of 61.7 years.When 1000 TSLF points were sampled from that TSLF map 1000 times, the mean FC was 61.6 years for FC MA (range: 57.9 to 66.9 years), and 61.7 years for FC MLE (range: 57.0 to 66.3 years).The FC MA was not significantly different from FC MLE (p = 0.35) or FR S (p = 0.18; from Section 3.1), as indicated by two-tailed t-tests.

Examples of Calculating the Minimum Number of Sample-Points
With an increase in the number of sample-points from the TSLF map (Figure 3a, created using all 1200 fire maps and the 2100 × 2100 grid landscape) the FC MA estimate initially varied greatly and then stabilized (Figure 6).When an ADMM of 10% was employed 1000 times from this TSLF map, the minimum number of samples required ranged from 65 to 571, with a median of 147 and a 95th percentile of 239 sample-points (Figure 7).The frequency distribution for 1000 simulations of the minimum number of TSLF sample-points required to estimate FC MA with an ADMM of 10% when using the TSLF map for a 2100 × 2100 grid simulated landscape that included all 1200 fire-year maps.The minimum number of TSLF sample-points is 65, median is 147, the 95th percentile is 239, and maximum is 571.

Creation of TSLF Maps with Different FR S
The 201 TSLF maps had an FR S that ranged from a low of 61.1 years when all 200 Class 1 fire years were employed, to a high of 200.5 years when no Class 1 fire years were employed.TSLF maps with fewer Class 1 fire years had progressively less young forest and more old forest (Figure 3).Spearman rank-correlations between 55 permutations of 11 TSLF maps ranged from −0.164 to 0.343 with a median of 0.026; when classified into bins that were 0.05 correlation-units wide, the mode was in the zero-bin of −0.024 to 0.025.The global Moran's I ranged from −0.127 to 0.120, and was positively correlated with FR S for those maps (r s = 0.964, p = 0.0000, n = 11), but not correlated with the 95th percentile of TSLF sample-points for those maps (r s = −0.092,p = 0.789, n = 11).

Influence of Study Area Scale on Estimated Fire Cycle (FC MA ) and TSLF Age Variability
At a study area scale of 3× of the largest total annual area burned, TSLF maps based on a higher FR S also had a higher FC MA estimated from the TSLF map (Figure 8).In contrast, at a study area scale of 0.5×, TSLF maps based on a higher FR S did not yield higher FC MA estimates.When those calculations were done for all nine scales (Figure 9), TSLF maps with a study area scale of ≥2× had a high R 2 adj and low RMSE between FR S and FC MA , while scales of ≤1× had a low R 2 adj and high RMSE.When similar analyses were conducted between the FR S and the standard deviation of the 201 TSLF maps (Figure 10), TSLF maps with a study area scale of ≥2× had a high R 2 adj and low RMSE between FR S and standard deviation, and study area scales of ≤1× had a low R 2 adj and high RMSE.

Influence of FR and Study Area Scale on Minimum Required Sample-Points
At a study area scale of 3× of the largest total annual area burned, TLSF maps based on a longer FR S did not exhibit any difference between FR S for the whole study area and FC MA estimates from the 3× scale TSLF map (Figure 11).In contrast, at a study area scale of 0.5×, TSLF maps based on a longer FR S did exhibit differences between FR S for the whole study area and FC MA estimates from the 0.5× scale TSLF map; the differences were typically negative for low values of FR S and positive for high values of FR S .Results from linear regressions between FR S and the FC MA from the nine different scales (Table 3) indicate that scales ≥4× had intercept values that were not different than zero, while progressively smaller scales had progressively more negative intercepts.The regression results also indicate that scales ≥3X had slopes that were not significantly different than zero, while progressively smaller scales had steeper slopes that were significantly different than zero.
At a study area scale of 3× of the largest total annual area burned, TSLF maps based on a longer FR S did not require more TSLF sample-points to obtain an ADMM of 10% than did TSLF maps based on a shorter FR S (Figure 12).There was no significant trend with FR S in the number of sample-points for the 0.5× study area, nor for the density of sample-points for the 3× or 0.5× study areas.The 3X study area scale did require more TSLF total sample-points, but at a lower density of sample-points per 100 km 2 , than did the 0.5× study area scale.In general, an increase in study area scale resulted in a linear increase in the minimum total number of required TSLF sample-points, and a negative exponential decrease in the number of sample-points per 100 km 2 (Figure 13).Each sample-point pairs the FR S from the whole 2100 × 2100 grid simulated landscape with the 95th percentile estimate of the 1000 trials of how many TSLF points are required; this was done for all 201 TSLF maps.The minimum number of TSLF sample-points is expressed in two ways: the total number for that study area scale, and the density of sample-points (points/100 km 2 ).Of the nine scales at which these examples were conducted, two are depicted here: 3× and 0.5×.

Discussion
Our key findings suggest that: TSLF fire histories should employ study areas ≥3× the largest total annual area burned [7]; the minimum number of TSLF sample numbers and sample density that should be employed differs with study area scale (study area divided by largest total annual area burned); and an influence of FR S length on the FC estimate only occurs if the study area scale is <3×.Although the 201 TSLF maps were created by randomized resorting of the same 1200 simulated fire-year maps, Spearman correlations suggest that the TSLF maps contain unique spatial patterns.Our results most directly apply to areas of boreal forest that are predominantly burned by large stand-replacing fires, but scaling relations could be used as a first approximation for study areas with small stand-replacing fires.
Three aspects of our results support the suggestion that study areas should be ≥3× the size of the largest total annual area burned [7].First, relations between FR S and the FC MA are best for scales of ≥2×.That is, for study area scales ≥2×, the estimated FC MA was strongly similar to the FR S for the complete 2100 × 2100 grid simulated landscapes.Remember that since the TSLF map is based on the FR S , a good FC MA estimate should equal the FR S .Second, the individual regressions between FR S against FR S minus FC MA at each study area scale indicate that only scales of ≥4× had an intercept not different than zero, and only scales ≥3× had a slope not different than zero.Third, standard deviation in TSLF ages is most similar to FR S for scales of ≥2×.Since a fundamental property of the negative exponential distribution that many TSLF distributions display, is that the mean and standard deviation should have equal values, it is not surprising that our scaling results were similar for the two parameters.Thus, while our results suggest that a study area scale as small as 2× could provide accurate estimates of the FC, since some of our relations only become stable at 4×, our findings agree with others that progressively larger areas provide more accurate results [14][15][16].Also, since a small study area scale will likely provide a TSLF distribution that is not representative of the FR, increasing the number of TSLF sample-points cannot fix that misrepresentation.
Given that real landscapes contain spatial heterogeneity in environmental conditions, unlike our spatially homogenous simulated area, study area scales >3× are likely required to obtain accurate FC estimates [50].Of the 18 TSLF fire history studies whose size and sample numbers were given in the Introduction, five provided information on the largest total annual area burned; they had a median study area scale of 4.1 and a range from 1.7 to 16.2 [8,22,27,28,32].This suggests that TSLF fire history studies should, to reduce possible errors in their results and to maximize the interpretive potential of their data, strive to study larger areas.In study areas restricted by property boundaries such as the edge of a park that make larger study area scales infeasible, recognition should be made that the FC MA could misrepresent the real FC even if the TSLF sample size is high.
The number and density of TSLF samples required to estimate the FC with an error (ADMM) of less than 10% for study area scales of 3× were, respectively, 187 samples and 0.81 per 100 km 2 .While this coarse-scale sampling will provide information on the regional FC, it will not provide enough information to understand small-scale drivers.An increase in sample number or density would help with detecting environmental small-scale drivers of spatial variations in FC, e.g., [18,27].If we assume that the 18 TSLF-based fire history studies mentioned in the Introduction had a study area scale of 3×, then two of them did not meet the sample density requirement, and nine did not meet the sample number requirement.This indicates that increased TSLF sample sizes should be employed to reduce error in FC estimates.
Our results suggest that only study areas with a scale of <3× will have FC estimates influenced by the length of the FR.In general, small-scaled areas with a short FR have an overestimated FC, and small-scaled areas with a long FR have an underestimated FC.This pattern is supported by the only studies we are aware of that allow a multi-scale comparison of FC and FR: the FC of 47 years from a 1.1× scaled study area by Van Wagner [3] underestimates the FR of 204 years calculated from a 2.3× scaled study area by Heinselman [8] within which it was nested.
The reduced standard deviation of TSLF ages that we observed in smaller areas could result in uncommon old stands over-influencing the FC estimate in areas with a short FR, and uncommon young stands over-influencing the FC estimate in areas with a long FR.It is possible, however, that field-based TSLF studies could avoid those biases as the exact area of recent burns can be known from historical records and/or remotely sensed imagery, and very old burns would be lost due to mortality of the original tree cohort and be replaced by younger successional forests, e.g., [13].
We recognize that our simulation methods and study design have some limitations, and could be improved in four key regards.First, our uniform landscape did not consider spatial variation in factors such as soil, topography, hydrography and climate, known to influence FC [27,51].Future studies could assess how requirements of study area scale and TSLF number change with increased spatial heterogeneity.Our preliminary research using neutral landscape models showed that increases in positive spatial autocorrelation did require increased number of TSLF sample points (Supplementary Materials), but variations in spatial autocorrelation in our 11 TSLF maps did not appear to influence the required number of TSLF sample points.A second improvement would involve modifying LANDIS-II and DFSS (Dynamic fuels and fire system) so that each simulated year could sample from different probability distributions of fire occurrence, fire size and fire weather.Since our preliminary analyses showed that using only one set of fire parameters resulted in unrealistically low inter-annual variation in area burned, we employed six classes of fire parameters and thus had to simulate years individually.If LANDIS-II and DFSS included this modification, that would allow one TSLF map to develop over 1200 years of burning, and thus incorporate the influences of forest age on fire spread [52].Third, our TSLF maps did not assume censoring due to mortality, and thus unrealistically contained TSLF grids >801 years.We could have applied a censoring to TSLF points, e.g., [13] but we instead believe that fire history field researchers should examine tree-rings to determine if it is an open-grown post-fire stand, or a closed-grown successional stand [53].If the latter occurs then either soil charcoal should be used to get the TSLF, e.g., [54], or the TSLFs should be examined as a truncated distribution [7].
A fourth improvement would be to evaluate the influence of age-related changes in the hazard of burning on the required number of TSLF sample-points.Many but not all forest types exhibit an increase in the hazard of burning with increased age [55] which can cause increased positive spatial autocorrelation of TSLF [52].Since increasingly positive spatial autocorrelation requires an increased number of sample-points to obtain a spatially representative mean value of landscape conditions [56], forests that have age-related increase in the hazard of burning might require more TSLF sample-points.Indeed, in a preliminary evaluation of this issue using neutral landscape models, we found that increases in spatial clustering of forest ages, such as might occur due to age-related increases in hazard of burning, did result in an increase in the required number of TSLF sample-points (Supplementary Materials).
We suggest six different ways that future TSLF fire history studies could apply our results to estimate the FC in forests with a stand-replacing fire regime.First, and most importantly, researchers should always use a study area ≥3× the largest total annual area burned.If the study area is under fire suppression, then they should use the largest total annual area burned from pre-suppression TSLF records or from other fire history studies.Second, researchers could use the minimum number of TSLF sample-points that we suggest for the relevant study area scale.Third, researchers should increase the number of TSLF sample-points above our minimum to account for spatial and temporal heterogeneity, such as caused by age-or environment-related increases in the hazard of burning (Supplementary Materials).Fourth, researchers could divide our TSLF sample density requirement by the appropriate scaling factor related to largest total annual area burned.For example if a study area had the largest total annual area burned 1/10 that in our study area, then instead of 0.81 samples per 100 km 2 they would use 0.81 samples per 10 km 2 .This approach does, however, simplistically assume similar slopes of the total annual area burned distributions.Fifth, researchers could import their spatial environment into LANDIS-II and DFSS, and parameterize the simulations in a manner that would reflect that area's temporal variations in fire weather.They could then use our Python code to determine either the number of TSLF points required to reach an acceptable error, or determine the error their FC estimate will have given their constraints in study area and sample number.Since all study areas contain unique spatial patterns of environmental factors that influence the hazard of burning, and unique age-related changes in the hazard of burning, this approach will yield the most realistic indication of the required number of TSLF sample-points.Sixth, an ongoing TSLF fire history could begin with our minimum suggested sample density, and employ a Monte Carlo simulation while doing fieldwork to estimate the FC and the error in that estimate; if the error is too high then the researchers would gather additional TSLF samples until the error was acceptable.

Conclusions
The method of constructing and sampling TSLF maps that we developed allows guidelines to be tested in regard to how large of a study area and how many TSLF sample-points should be employed to estimate the FC using the TSLF method.This approach builds upon previous simulation studies which showed that increased sampling effort would reduce error in FC estimates.
Our results indicate that TSLF fire histories should be created for areas that are ≥3× the largest total annual area burned that occurs in the region within which the study area is located.If smaller study areas are employed, then areas with a short but unknown FR will likely obtain an overestimated FC, and areas with a long but unknown FR will likely obtain an underestimated FC.Since these biases will be the result of the study area being too small, they cannot be compensated for by employing more TSLF samples.
Our results further indicate that in forests whose fire regime is similar to our simulation, 187 TSLF sample-points is the minimum required to estimate the FC with an error of less than 10% in a study area 3× the largest total annual area burned.Larger areas would require more samples but at a lower spatial density.Since our suggestions are based on a spatially homogenous study area with a cyclic climate, study areas with progressively higher spatial and temporal heterogeneity due to factors such as age-related increases in the hazard of burning, should employ progressively more TSLF samples.To optimize the number of TSLF samples in a fire history study, we encourage researchers to develop site-specific applications of our simulation approach prior to conducting fieldwork.

Figure 2 .
Figure 2.An example of (a) the stacking process of annual fire-year maps for the six classes of total annual area burned to create (b) a simulated time since last fire (TSLF) map.In the simulated TSLF map, warm colors represent area burned in recent years, while cold colors show areas burned longer ago.

Figure 4 .
Figure 4.The rank-ordered distributions of the percent area of the landscape burned in the 44 years of historical fire records for Wood Buffalo National Park (WBNP) and the 1200 years of simulated data, and their exponential regression parameters.The y-axis is truncated at 0.1% as simulated data with fire-years lower than that all had zero area burned.

Figure 5 .
Figure 5.The cumulative TSLF distribution and regression model created for the 2100 × 2100 grid simulated landscape using all 1200 fire-year maps.The FC MA estimate from the slope is 61.7 years.

Figure 6 .
Figure 6.Three examples of how the FC MA estimate varied with the number of TSLF sample-points, using the TSLF map for the 2100 × 2100 grid simulated landscape that contained all 1200 fire-year maps.For example, the minimum required number of sample-points to obtain an absolute difference in the measured mean (ADMM) of 10% of the FC MA estimate is (a) 65, (b) 147 and (c) 571.Those values correspond to the minimum, median and maximum values in the distribution in Figure 7.

Figure 7 .
Figure 7.The frequency distribution for 1000 simulations of the minimum number of TSLF sample-points required to estimate FC MA with an ADMM of 10% when using the TSLF map for a 2100 × 2100 grid simulated landscape that included all 1200 fire-year maps.The minimum number of TSLF sample-points is 65, median is 147, the 95th percentile is 239, and maximum is 571.

Figure 8 .
Figure 8. Relations between simulated fire rotation (FR S ) for the whole 2100 × 2100 grid simulated landscape and FC MA estimates for each of the 201 TSLF maps.Of the nine scales at which these analyses were conducted, two are depicted here: 3× and 0.5×.Each point represents one of the 201 TSLF maps.The R 2 adj and root mean square error (RMSE) are, respectively, 0.97 and 7.7 for the 3× study area scale, and 0.03 and 54.4 for the 0.5× study area scale; these values provide four of the points used in Figure 9.

Figure 9 .
Figure 9.The R 2 adj and the RMSE between FR S for the whole 2100 × 2100 grid simulated landscape and the FC MA estimate for each of the 201 TSLF maps at each of the nine study area scales.

Figure 10 .
Figure 10.The R 2 adj and RMSE between FR S for the whole 2100 × 2100 grid simulated landscape and the standard deviation of TSLF ages for all 201 TSLF maps at each of the nine study area scales.

Figure 11 .
Figure 11.Relations between FR S for the whole 2100 × 2100 grid simulated landscape and FR S minus the FC MA estimate for each of the 201 TSLF maps at the relevant study area scale.Of the nine scales at which these examples were conducted, two are depicted here: 3X and 0.5X.

Figure 12 .
Figure 12.The minimum number of TSLF sample-points to estimate FC MA with an ADMM of 10%.Each sample-point pairs the FR S from the whole 2100 × 2100 grid simulated landscape with the 95th percentile estimate of the 1000 trials of how many TSLF points are required; this was done for all 201 TSLF maps.The minimum number of TSLF sample-points is expressed in two ways: the total number for that study area scale, and the density of sample-points (points/100 km 2 ).Of the nine scales at which these examples were conducted, two are depicted here: 3× and 0.5×.

Figure 13 .
Figure13.The minimum number of TSLF sample-points, and the density of TSLF sample-points per 100 km 2 , which were obtained for each of the nine study area scales.The number of TSLF points at each scale is the 95th percentile from the 201 TSLF maps, each of which in turn is represented by the 95th percentile of their 1000 trials of how many TSLF samples were required to meet an ADMM of 10%.

Table 1 .
Fire characteristics in six classes of total annual area burned (TAAB) for 44 years of historical records from Wood Buffalo National Park and 1200 years of simulated fires.

Table 2 .
Parameters of fire size distribution and annual ignition number for individual fires that occurred in each of the six classes of total annual burned (TAAB), including mean (µ) and standard deviation (σ) of individual fire size lognormal distribution, maximum individual fire size (α), and annual number of individual fires (η).Values are in hectares as that unit is employed in LANDIS-II.

Table 3 .
Linear regression results between the FR S for the whole 2100 × 2100 grid simulated landscape and FR S minus the FC MA for each of the 201 TSLF maps using the relevant study area scale.