Postfire recruitment failure in Scots pine forests of southern Siberia

Wildfire disturbances effect changes in vegetation communities that in turn influence climate. Such changes in boreal forest ecosystems can persist over decadal time scales or longer. In the ecotone between boreal forest and steppe in the region southeast of Lake Baikal in southern Siberia, shifts between the two vegetation types may be precipitated by variations in site specific conditions, as well as disturbance characteristics such as fire frequency and severity. Warmer, drier conditions in the region have been associated with a decrease in fire return intervals and greater burn severity that may, in turn, drive conversion of forests to steppe vegetation at a greater rate than has occurred prior to the onset of warming and drying. Stand-replacing fires in Pinus sylvestris stands in southern Siberia may lead to recruitment failure postfire, particularly on southwest to west-facing slopes, which are more often dominated by grasses. This study uses a combination of field data and remotely sensed indices of vege- tation and moisture to distinguish between recruitment pathways in southern Siberia, and to study the influence of factors related to soils, topography, fire severity and winter snow cover on these. We expected that recruitment success would be associated with lower burn severity (higher NBR), higher greenness (NDVI) and moisture (NDMI), and winter snow (NDSI) postfire. We also expected phenological characteristics to differ among recruitment paths. Prior to burning, our sites are broadly similar in terms of remotely sensed indices of moisture (NDMI), vegetation (NDVI), and winter fractional snow cover (NDSI), but recruitment failure sites are generally drier and less green postfire. Initial differences in greenness and moisture among sites characterized by abundant recruitment (AR), intermediate recruitment (IR) and recruitment failure (RF) become more pronounced over the initial decades postfire. The earliest separability of AR and RF sites using remotely sensed indices occurs in the winter months 3–4 years postfire, during which time NDSI is highest for AR sites and lowest for RF. Although seasonality was important with regard to distinguishing among AR, IR and RF index values, the timing of phenological events such as start and end of season did not differ significantly among the sites.


Introduction
Boreal forests store 40 Pg C (Thurner et al., 2014) and account for 20% of the global forest C sink (Pan et al., 2011). Eurasian and North American forests each account for about half of boreal carbon storage (Thurner et al., 2014) and the northern land sink (Gurney et al., 2002). Boreal forest carbon stocks are increasing due to its expansion into tundra and former agricultural areas  as well as increases in ecosystem productivity from longer growing seasons and warmer temperatures (White et al., 1999). Boreal carbon stocks are also vulnerable to losses from increased respiration (Quegan et al., 2011), droughts (Ma et al., 2012;Michaelian et al., 2011;Peng et al., 2011), and more frequent and severe disturbance events (Gustafson et al., 2010) including wildfires Settele et al., 2014).
Boreal fire/climate feedbacks and their magnitudes and time scales are still poorly understood. Disturbance cycles are associated with exacerbating and mitigating feedbacks to climate change (Johnstone et al., 2016;Rogers et al., 2013) (Settele et al., 2014). More frequent and more severe fires lead not only to substantial boreal forest carbon losses Miquelajauregui et al., 2018;Stephens et al., 2014) and possibly the reversal of a boreal carbon sink (Bradshaw and Warkentin, 2015), but also provoke shifts from forests to non-forested ecosystems (Scheffer et al., 2012), particularly in areas that have experienced climate change pressures for an extended period prior to disturbance (Seidl et al., 2017). Persistent forest losses may occur because increased fire frequency interrupts relay succession dynamics (Johnstone and Chapin, 2006a), or because conditions such as shading (Holmgren et al., 1997), soil moisture and temperature or organic layer depth are so modified that tree seedlings fail to become established and other plant functional types are more successful (Johnstone and Chapin, 2006b;Trugman et al., 2016).
Siberia has experienced the greatest temperature increase in the northern hemisphere (Groisman et al., 2013), and southern Eurasian boreal forests are contracting their range, even as they expand to higher elevations and further north (Koven, 2013), consistent with predicted shifts associated with climate change in the region (Tchebakova et al., 2009). In this region, subtle changes in temporal precipitation patterns together with longer fire seasons due to rising temperatures can trigger a tipping point in the ecological systems that leads to the demise of forests and the expansion of steppe vegetation (Tchebakova et al., 2011). Fire frequency and precipitation are the two primary drivers that maintain forest or grassland ecosystems according to observational data (Bowman et al., 2013;Hirota et al., 2011;Mayer and Khalyani, 2011), models of vegetation change (Brazhnik et al., 2017, Lehsten et al., 2015Tchebakova et al., 2011), and ice core reconstructions of vegetation and climate (Eichler et al., 2011). Southern Siberia is on a gradient between much cooler and wetter taiga and the drier and warmer Mongolian steppe. Fire return intervals (FRI) in the ecotone between taiga and steppe vegetation are on the order of 25 to 50 years (Chu et al., 2016;Furyaev, 1996;McRae et al., 2006;Swetnam, 1996;Wirth et al., 1999), considerably shorter than FRI further north (> 100 years [Wirth et al., 2002]) and longer than intervals to the south (7-16 years [Hessl et al., 2012]). The landscape in southern Siberia is generally a patchwork of forest and grasses, with forests found primarily on slopes, and the discontinuity of forest vegetation cover indicates a potential for ecological instability (Peterson, 2002).
The potential for shifts between forest and grassland-type ecosystems, or bi-stability between these two states have been observed globally (Hirota et al., 2011) in boreal regions of Eurasia (Scheffer et al., 2012;Tchebakova et al., 2009), tropical regions in South America (Betts et al., 2004;Malhi et al., 2008;Oyama and Nobre, 2003), Africa (Lehsten et al., 2015;Staver et al., 2017), and Australia (Bowman et al., 2001). Dominant plant functional type is largely governed by feedbacks K. Barrett, et al. Remote Sensing of Environment 237 (2020) 111539 by which tree growth is either promoted or hindered (Adamek et al., 2016), and fire is frequently the catalyst for such shifts (Bowman et al., 2013;Tarancon et al., 2014). In Russian boreal forests, as much as a third of burned area fails to regenerate postfire . Persistent forest losses postfire have been commonly observed in southern Siberia , where logging impacts fuel loads available for combustion and in turn fire severity (Kukavskaya et al., 2016). Tree mortality dynamics in the southern boreal forest are critical to understanding carbon cycling and climate feedbacks in the region (Bradshaw and Warkentin, 2015;Gustafson et al., 2010;Koven, 2013).
In this study, we seek to distinguish among areas in southern Siberia that recruit trees abundantly postfire, intermediate recruitment, and recruitment failure (AR, IR, and RF, respectively) using remotely sensed data and information from field observations. We use long time series of remotely sensed indices of greenness, moisture, snow cover, and fire severity from Landsat and phenological characteristics derived from MODIS indices to compare these recruitment types prior to burning and for up to 17 years postfire. We test hypotheses that differences in recruitment type are caused by higher burn severity and lower winter snow cover, and that variability in recruitment leads to differences in greenness and vegetation moisture as well as phenological characteristics of AR, IR, and RF sites.

Study area
Our study area is in the southern Eurasian boreal forest range, in the Zabaikalsky Krai and the Republic of Buryatia, immediately southeast of Lake Baikal (Fig. 1). The territory is characterized by highly-dissected and eroded mountainous terrain with mountain ridges elongated from south-west to north-east. Altitudes range from 300 m in the valleys of the Zabaikalsky Krai to almost 3500 m above sea level in the Eastern Sayan Mountains of the Republic of Buryatia. Mean monthly temperature varies from a low of −24 to −26°C in January to a high of +15 to +18°C in July. The average annual precipitation varies from 200 mm at southern forest, forest-steppe and steppe zones to 700-900 mm in the Stanovoy highlands and in the Eastern Sayan Mountains (Geniatulin, 2000;Gerasimov, 1965). Summer (mainly July and August) accounts for 60 to 70% of annual precipitation (Forest Plan of the Republic of Buryatia, 2008; Forest Plan of the Zabaikal Region, 2014). Forests in the region are mountain taiga, typically light conifer and deciduous stands, mainly consisting of larch (Larix gmelinii and L. sibirica), Scots pine (Pinus sylvestris), and to a lesser extent Siberian pine (Pinus sibirica), aspen (Populus spp.) and birch (Betula spp.). Pine stands are widespread on dry and mesic sandy and sandy-loam soils. Permafrost is primarily confined to the northern part of the region and to the upper part of north facing slopes, where pine is not found. Larch stands grow in the upper parts of north slopes and in the zone of continuous permafrost in the northern region. Mixed larch-pine stands with some component of deciduous species are typical on mesic and wet loamy soils (Buryak, 2015). In the foreststeppe, south-facing slopes are occupied by steppes, and north-facing slopes by Scots pine or birch (Betula pendula).
Fires have dominated Siberian boreal forests for the last 9000 years (Katamura et al., 2009), and, similar to the association with black spruce in Alaska, their prevalence coincided with the spread of boreal woodlands 11,000 to 10,000 years ago (Bezrukova et al., 2010). Siberian wildfires are typically lower severity surface fires (Rogers et al., 2015), and the proportion of stand-replacing fires in Siberia is lower in the south (Krylov et al., 2014). Climate and wildfire in the region are less directly linked than in the boreal forests of North America due to anthropogenic impacts including ignitions  and effects on available fuels through forest management (Kukavskaya et al., 2013) and through livestock grazing (Hessl et al., 2012). The disconnection between climate and fire is evidenced by low correlations between fire activity and weather/climate conditions (Eichler et al., 2011), plus limited synchrony in fire events (Hessl et al., 2012, Swetnam, 1996. Fires in the region are characterized by high density and high annual variability during a short fire season (Chuvieco et al., 2008), and compared with global fire regimes are rare, intense, and large (Archibald et al., 2013). Large-scale climate patterns such as the Arctic Oscillation can influence interannual fire variability (Balzter et al., 2007).
Scots pine spread through the region about 7000 years ago, corresponding to a shift towards a drier and colder climate similar to present day (Bezrukova et al., 2010). Scots pine forests in the region are characterized by the highest level of disturbance and fire hazard as they grow in the drier slopes with greater insolation and in the lower part of the slopes with a more developed road network (Buryak, 2015). Fire return intervals in Siberian Scots pine stands are on the order of 20-40 years (Ivanova et al., 2010;Wirth et al., 1999). Pine stands are managed for timber and firewood, and most burned forests that are accessible by road are normally also logged, and sometimes replanted postfire, as they are elsewhere in the Siberian boreal forest (Hu et al., 2015). Illegal logging in the region significantly increase the amount of fuels available to burn because of slash left on site, thus increasing potential fire severity and the probability of fire spreading to tree crowns (Kukavskaya et al., 2013).

Methods
Our sites mostly extend along the Trans-Siberian highway between Chita and Ulan Ude. We depended on expert knowledge of foresters and ecologists in the region to locate pine stands that burned after 2000, and prefire forest type was confirmed by identification of partially-burned stumps. Sites were generally restricted to south and sometimes west-facing slopes, where pine stands are typically found. Sites were characterized in situ according to an initial assessment of recruitment levels as AR (n = 14), IR (n = 37), and RF (n = 13). These characterizations were later confirmed as AR and RF sites were at the extremes of the variability in stocking density among the sites (Figs. 2 and 3). The stocking density for each site was converted to a "tall seedlings" (i.e., > 1.5 m) equivalent using Russian forestry standards (Russian Forest Restoration Rules, 2016), and compared with the required stocking density to support recruitment for a given site based on ecoregion, species and forest type and soil moisture. Sites that did not meet the stocking density of tall seedlings requirement were maintained as RF, sites that had more than the requirement but less than two times the requirement were maintained as IR, and sites with more than twice the required stocking density were maintained as AR. Field sites that did not correspond to the appropriate stocking density (n RF = 0, n IR = 18, n AR = 1) were omitted from the analysis.

Hypotheses
We tested the following hypotheses regarding the characteristics of AR, IR, and RF sites. First, that RF sites would exhibit the lowest postfire vegetation index (NDVI) values, followed by IR and AR sites due to differences in leaf area, green biomass, canopy structure and photosynthetic activity (Carlson and Ripley, 1997;Dong et al., 2003;Gamon et al., 1995). Second, that RF sites have experienced more severe fires (i.e., lower NBR) than IR and AR sites, which have led to RF through the combustion of seed stock and viable plant propagules (Flinn and Wein, 1977). Third, that RF sites are characterized by lower vegetation moisture levels (NDMI) than IR and AR sites due to differences in species and biomass levels. Fourth, that fractional winter snow cover (NDSI) is lower in RF sites than IR and AR sites, possibly leaving seedlings in RF sites vulnerable to colder temperatures and stronger winds in winter (Myers-Smith and Hik, 2013;Sturm et al., 2001). Finally, that the timing and magnitude of phenological events, which depend on plant functional type (Bonan et al., 2002) would vary among AR, IR, and RF sites due to differences in vegetation communities.

Field methods
Field sites were located in an area that appeared homogeneous in terms of stand characteristics such as stocking density and stand age for an area at least 250 × 250 m. Soil texture was characterized generally as light, medium, or sandy loam, or sand, and moisture as either "mesic" or "dry", where pine stands are most commonly found in the region (Wirth et al., 2002). A nested transect was used for comparison with both Landsat and MODIS imagery. A 200 m transect was sampled every 30 m for comparison with MODIS data, and additional subplots every 10 m for the first 30 m of the transect for comparison with Landsat data. Postfire recruitment was assessed in ten 1 m 2 quadrats in which stems of regenerating trees were counted. In these quadrats, total ground cover was estimated as a percentage, and fractional cover was estimated by species. Prefire vegetation was assessed for a larger area in three 10 m 2 circular plots, in which diameter at base height, as well as height and diameter at breast height (1.3 m) if more than a stump remained, was measured for all dead and live trees.

RS methods
Many of our sites had recently burned more than once, and determining time since fire was a challenge for sites that were known to have burned multiple times. We performed all of the analyses described using the first and then the last fire known to have occurred as the date of the burn, and we report the values derived using the most recent fire based on somewhat higher separability of AR and RF sites using these values (but see Appendix A for the results of the analysis using the first known fire to have occurred). The use of oldest fire year did not affect temporal patterns, direction or magnitude of differences among indices for AR, IR and RF sites. We assessed time series of remotely sensed indices of vegetation greenness, fire severity, moisture, and snow as a function of time since fire, as well as seasonal variability for the periods prior to burning ("prefire"), 1-2 years postfire ("initial assessment"), and > 12 years postfire ("extended assessment"). These periods are associated with early colonization by herbaceous species, followed by expansion of woody plants such as shrubs and tree species, and eventual dominance of trees if present (Chu et al., 2016;Frazier et al., 2015;Yang et al., 2017). Upon evaluating the yearly time series data of indices and phenological characteristics, we added a fourth, "intermediate assessment" period 3-4 years postfire, due to apparent nonlinearities in the recovery of indices. Observations were classified according to season: spring (doy 60-150) growing season (doy 150-240), autumn (doy 240-300), winter (doy > 300 or doy < 60).

Indices
The Landsat data used in the analysis were from the pre-Collection surface reflectance products combined from Landsat 5, 7 and 8 (sensors TM, ETM+, and OLI, respectively). We attempted the same analysis with Collection 1 data, but the number of observations was insufficient for use of the methods presented. We used indices of greenness (Normalized Difference Vegetation Index NDVI) (Tucker, 1979), moisture (Normalized Difference Moisture Index, NDMI) (Wilson and Sader, 2002), burn severity (Normalized Burn Index, NBR) (Key and Benson, 1999) and snow (Normalized Difference Snow Index, NDSI) (Riggs et al., 1994) to test the hypothesized relationships between recruitment failure and the post-fire environment (Section 3.1). Cloud, cloud shadow, snow and water pixels were masked from the data using the data flags associated with the surface reflectance product using the C Function of Mask (Foga et al., 2017), and only clear pixel values were maintained in the analysis.
While the MODIS data products are coarser in spatial resolution (500 m) than the Landsat images (30 m), MODIS indices are available daily, which is preferable for studying phenological characteristics when compared with the 16-day return time of Landsat. NDVI and NDMI were calculated using the version 6 MODIS Nadir Bidirectional Reflectance Distribution Function (BRDF) Adjusted Resolution (NBAR) bands (MCD43A4). The daily values are produced from the most representative observation in a 16-day moving window for both Aqua and Terra sensors, and assigned to the ninth day in the time window (Schaaf and Wang, 2015). These corrected MODIS bands were used to control for the strong BRDF effects at higher latitudes, that can impact even ratio-type indices such as those used in this analysis (Lucht and Lewis, 2000;Roman et al., 2011;Verbyla et al., 2008).

Phenology
Phenology characteristics for AR, IR, and RF sites may differ due to the distinct plant functional types that dominate the different trajectories and associated impacts on phenological patterns, a phenomenon observed in persistent postfire vegetation changes in tundra (Barrett et   K. Barrett, et al. Remote Sensing of Environment 237 (2020) 111539 length of season, peak value, position of peak, trough, position of trough, mean autumn, spring and growing season values) were derived from the BRDF-corrected indices for AR, IR, and RF sites as a function of year since fire. Annual metrics were used to control for the possibility that phenological characteristics change over time during recovery from wildfire disturbance. Phenology was modelled using the greenbrown package in R (Forkel and Wutzler, 2015) for each year in relation to the fire. The index values for the year were first smoothed and interpolated using a mean function and the temporal filling and gap smoothing spline method (Forkel and Wutzler, 2015). Annual phenological characteristics were calculated using a threshold of maximum and minimum index values for each site using the methods described in White et al. (1997) with the smoothed and interpolated NDVI values.

Separability
Separability of two distributions can be assessed using several different distance calculation methods including Bhattacharyya (Bhattacharyya, 1943;Kailath, 1967) and Jeffries-Matusita distance (Bruzzone et al., 1995;Swain and King, 1973), Divergence and Transformed Divergence indices (Singh, 1984;Swain and King, 1973). The separability metric used in this analysis is the M-statistic, originally developed for differentiating between reflectance values of forest and non-forest pixels in different AVHRR wavelengths (Kaufman and Remer, 1994). The M-statistic, calculated using the spatialEco package in R (Evans and Ram, 2017), was chosen for having a transparent critical threshold level above which separability is deemed acceptable. M separability between two distributions for groups A and B, is calculated as follows: M is > 1 when the mean of two normal distributions are separated by more than one standard deviation, or when the assignment of observations to a class (e.g., winter NDSI values of AR and RF sites 1 year postfire) are unambiguous for 84% of the observations (Fig. 4). Other measures of divergence were evaluated, and generally results were similar, but without a threshold for evaluating separability such as M > 1, it was difficult to compare outputs or to determine when separability became sufficient.
We focused upon the separability of Landsat-derived NDVI, NBR, NDMI, and NDSI for AR and RF sites (as the two endpoints of the recruitment level scale) on a seasonal basis as a function of time since fire.
The seasonal time step was a balance between having a fine enough temporal resolution to observe sub-annual changes, and a long enough time step to accumulate a sufficient observations for comparison. NDSI was evaluated as the annual averages of winter values. In addition to the separability using different indices, we assessed yearly separability using the phenological characteristics described above. The Shapiro-Wilk test was used to exclude dates for which index values were nonnormally distributed. Landsat data were available for a maximum of 16 years prefire, and a maximum of 17 years post fire, with an average of 10, 9, and 9 observations (min = 1, max = 26, 29, 28) per year for AR, IR, and RF sites, respectively. Daily MODIS data were available for a maximum of four years prior to the fire and up to 17 years postfire, with 365 direct or modelled observations per year.

Field data
Prefire stocking density was similar across AR, IR and RF sites (Fig. 3). RF sites had somewhat lower median stocking density prior to burning, possibly because stumps were consumed by the fire, but the differences among sites were not statistically significant. Soils in AR and IR sites were more often characterized as mesic than dry (10 out of 14, and 20 out of 37, respectively), whereas soils in RF sites were more likely to be dry (11 out of 13). Soil texture in AR sites was less likely to be characterized as sandy, and more often associated with loam, whereas RF sites were more likely to be sandy and IR sites had a mix of soil types, most often characterized as sandy loam (Fig. 5). RF sites had a thinner organic layer (mean = 0.89 cm, sd = 1.13) than AR (mean = 2.75 cm, sd = 2.18) and IR sites (mean = 1.71 cm, sd = 2.75 cm) (Fig. 5). RF sites were more likely to be found on westand southwest-facing slopes and in flat areas, whereas AR and IR sites were located across south-, southeast-, southwest-, and west-facing slopes (Fig. 5). Fractional ground cover in AR, IR and RF sites was very similar, around 0.45 for all. AR sites were more often dominated by grasses, while RF sites were more likely to be dominated by forbs (Fig. 6). Shrubs were only observed in intermediate recruitment sites. Fig. 7 shows the progression of Landsat-derived mean growing season NDVI and NDMI as a function of time since fire. Both indices Fig. 3. Stocking density of trees prefire (left), and postfire tree regeneration (right) for AR, IR, and RF sites. Prefire differences were not significant, while postfire differences were.

Remotely sensed data
K. Barrett, et al. Remote Sensing of Environment 237 (2020) 111539 show a substantial decrease in the year of the fire, and indices for AR sites exhibit a positive slope postfire. IR sites showed a shallower increase in postfire NDVI, while RF sites do not show a trend. NDMI for IR and RF sites exhibits a negative slope postfire. Table 1 shows the separability of seasonal index values for AR and RF sites for immediate, intermediate, and extended assessment periods. The earliest separability occurs in the winter intermediate assessment period, during which the separability NDMI, NBR, and NDSI values for AR and RF is > 1. All indices were sufficiently separable in the growing season in the extended assessment period. Fig. 8 shows the distributions of winter indices for AR and RF sites. In this period. NDVI values are uncharacteristically low for AR sites, and NDMI is higher than in RF sites. NDSI values in AR sites were substantially higher than those for RF sites in the intermediate assessment period. NDSI values were not significantly different during other assessment periods (Table 1, Fig. 9).
Most phenological parameters did not show a trend or a separability > 1 at any time. Fig. 10 shows only the separability of MODIS NDVI-derived phenological parameters that exhibit values > 1, as a function of time since fire. None of these parameters (mean spring value, mean autumn value, mean growing season value, and peak value) are related to the timing of phenological events, but rather reflect the importance of seasonal variability in NDVI. Separability or AR and RF sites reached a separability value of > 1 after 15 years postfire. Interestingly, separability of all phenological parameters displayed exhibits an initial increase the year of the fire, a subsequent decrease 3-4 years postfire and a sustained increase 5 years after the fire.

Hypotheses
Differences in index values became most apparent in the extended assessment period > 12 years post-fire. As expected, RF sites exhibit the lowest NDVI values, followed by IR and AR sites (Fig. 7). Low NDVI   K. Barrett, et al. Remote Sensing of Environment 237 (2020) 111539 values for RF sites are likely explained by the low stocking density and lack of canopy cover in these sites (Gamon et al., 1995), given understory cover is similar among AR, IR, and RF sites (Fig. 6). NDMI is lowest for RF sites, and highest for AR sites (Fig. 7), differences that first become distinguishable in winter during the intermediate assessment period (Table 1). We believe the low NDMI values for RF sites is associated with the lack of canopy vegetation and therefore vegetation moisture (Hardisky et al., 1983). NBR values are lower for RF sites in the growing season and autumn during the extended assessment period, and in the winter during the intermediate period (Table 1). This is consistent with our hypothesis that fire severity is higher in RF sites, but we would have expected to see such differences earlier in the recovery period, i.e., in the immediate assessment period 1-2 years post-fire. Winter NDSI is lowest for RF sites and highest for AR sites, but only during the intermediate assessment period (Fig. 9). This result is consistent with our hypothesis that winter snow cover may be an important consideration for post-fire recruitment success.
NDVI, and NDMI were generally more useful than NBR or phenology in discriminating between AR and RF sites. It is possible that NBR is not as strongly related to fire severity in the region particularly as the fires were all stand replacing, so considerations that affect surface reflectance characteristics such as fractional mortality are less useful severity metrics than greenness and moisture levels. Other studies have highlighted the shortcomings of NBR to measure ecologically significant variations in fire severity (Lentile et al., 2006), particularly in boreal forests (French et al., 2008). The differences in organic layer depth among AR, IR, and RF sites (Fig. 5) are statistically significant, however, and could reflect more severe or more frequent fires in RF sites. While seasonality was important with regard to the separability of AR and RF sites, the timing of phenological events, such as start of season and end of season, were not useful in separating AR and RF sites. The magnitude, but not the timing of NDVI phenology varies among AR, IR, and RF sites (Fig. 10). We hypothesized that differences in the timing of phenological events may result from differences in plant functional type, for example differences in grass and tree canopy green up (Archibald and Scholes, 2007). The lack of separability may be because RF sites have similar phenology to AR sites, or that the data used to derive the phenology were insufficient, or that the phenological model was not a good fit to the data. The Collection 6 MODIS product (MCD43A4) uses a 16-day compositing period, days that appear to be consecutive could be as much as 17 days apart, which can affect the accuracy of phenological models (Ahl et al., 2006;De Beurs and Henebry, 2010). Although the models showed good visual agreement with time series of index values, it is difficult to ascertain how well these performed overall as there is no metric of goodness of fit available. The results presented here are not sufficient to rule out the existence or importance of differences in phenology between AR and RF sites, but they highlight the usefulness of comparing seasonal variability in index levels.

Importance of winter intermediate assessment period
We prioritize earlier detection of recruitment trajectories because of associations with forest health and the potential for early intervention, because as much as two thirds of the forested area in the region is K. Barrett, et al. Remote Sensing of Environment 237 (2020) 111539 managed (Gauthier et al., 2015). Differences between recruitment levels are detected earliest when seasonal differences are considered, especially winter values during the intermediate assessment period. Values across all indices generally became more similar during the intermediate assessment period, with the exception of winter values, which became more disparate. During this period, AR sites exhibited lower NDVI and higher NDMI and NDSI than RF sites. NDVI and NDMI are positively correlated for all other seasons in all other periods of assessment, but negative during this time, possibly due to greater fractional snow cover (evidenced by substantially higher NDSI levels) in AR sites. Photosynthetic activity, though it may occur in winter even under snow (Starr and Oberbauer, 2003), is obscured from the sensor by snow, and the area appears wetter (i.e., exhibits higher NDMI), likely due to the snow moisture content. Precipitation in the region falls mainly as rain in the mid to late summer, so although snow is a limited resource in the region, it is important in protecting regenerating vegetation and soils from extreme cold in winter (Sturm et al., 2001), as well as winter desiccation (Kharuk et al., 2010a) wind shear (Sturm et al., 2001), and herbivory (Hamilton et al., 1980;Tape et al., 2010). Snowmelt also provides critically important moisture in the spring (Westerling et al., 2006;Buermann et al., 2018), and forest trees trap snow in winter months, leading to greater snow moisture (Onuchin et al., 2018). Dendroecological studies in Siberia have documented previous mortality events and subsequent recruitment failure in the 16th and 18th Century associated with lack of snow and summer water stress (Kharuk et al., 2010a). Lack of snow induces apical shoot mortality, and for branches above snow level, chlorosis in Siberian pine (Kharuk et al., 2010a). Tree mortality generally occurs during the snow-free period (Tremblay et al., 2007), and the positive effect of snow on the expansion of woody vegetation is well-documented in tundra, where shrubs trap snow that encourages their expansion into areas dominated through mechanisms of protection from physical elements and insulation from cold temperatures (Halliger et al., 2010;Myers-Smith et al., 2011, Tape et al., 2006. It is possible that a similar mechanism operates in terms of promoting abundant postfire recruitment in the mountain taiga of southern Siberia.
The increasing separability of AR and RF sites over time suggests that a feedback mechanism may reinforce differences between the sites over time, as opposed to differences in recruitment being caused entirely by differences in fire severity, for example. While the RF sites are slightly drier and less green prior to burning, differences in NDMI and NDVI only become significant postfire, and could be consequence of recruitment failure as opposed to the cause. The increase in snow cover in AR sites 3 to 4 years postfire occurs during a period in which these sites are less distinguishable from RF sites in all other seasons, possibly because differences in vegetation have not yet become established. NDSI values in AR sites during the intermediate assessment period are atypical of the postfire environment, and approach prefire NDSI values which were generally higher, possibly because the snow was trapped and shaded by trees and understory vegetation.
NDSI does not generally correlate with any volumetric measurement of snow (e.g., depth, snow water equivalent), but it can be used to infer the fraction of area covered with snow (Hall and Riggs, 2007;Nolin, 2010). We note that the NDSI values observed are lower than the threshold used to identify snow-covered pixels (0.1-0.5, Hall et al., 1995), which may be due to low winter precipitation in the region and subsequently low fractional snow cover. Thresholds for differentiating snow from non-snow NDSI values vary substantially among ecosystems and regions (Appel, 2018), and negative NDSI values have been observed for some values of fractional snow cover in Siberia and elsewhere Appel, 2004, 2006;Appel, 2018). We are therefore unable to derive information on the depth of snow for a given area, and NDSI values are here used only as a proxy for fractional snow cover. Furthermore, the radiometric resolution of Landsat 5 primarily and 7 sensors is not well-suited to snow discrimination, but with improvements to the Landsat 8 sensor, it is likely that the relationship between fire severity and snow cover will be the focus of future research (Wang et al., 2016). Landsat 8 data were not used for this analysis because of the lack of historical data record.
Prefire NDSI values are lower for AR sites, meaning it is also unlikely that these are located in topographic positions that are better suited to snow accumulation. RF sites, being located more often on southwest and west facing sites, may lose more snow due to greater insolation when the day is warmest. The differences in NDSI do not appear in any other period, highlighting the potential importance of the timing of the increased snow cover as a function of time since fire. It is important to note that the intermediate period postfire occurs in different years for different sites, and it is therefore not the result of two extreme snow years but rather a feature that occurs in AR sites at a critical period in the recovery from wildfire disturbance. This, and the fact that the separability between AR and RF sites increases over time, raise the possibility that AR and RF sites may have followed a similar recruitment trajectory for the initial postfire period and diverged in the intermediate period due to factors external to site level characteristics, such as weather conditions. Weather has been observed to be a driver of ecosystem shifts both in interactions with wildfire disturbance (Trugman et al., 2016) and in the absence of disturbance (Holmgren et al., 2006;Kharuk et al., 2013Kharuk et al., , 2017. If weather in the initial assessment period is a factor in the recruitment trajectory, detecting differences between AR and RF sites may only be possible after several years of regrowth.

Broader ecological significance
The results of this analysis suggest that climate feedbacks exist with boreal postfire recruitment trajectories, in addition to documented impacts of disturbance frequency and severity and ecosystem productivity and respiration. The interaction between precipitation and temperature during recovery from fire is also likely to be an influential factor because of the importance of precipitation falling as snow or as rain (Berghuijs et al., 2014;Stenseth et al., 2002;Trenberth, 2011), and evapotranspiration effects when temperatures are warmer (Seneviratne et al., 2010). The feedbacks between climate and postfire recruitment trajectories (Brazhnik et al., 2017) may contribute to those factors that influence climate drying through a rapid transformation of forests to open woodlands (Scheffer et al., 2012). The drier and browner conditions associated with recruitment failure sites correspond to model predictions regarding forest losses in the region (Tchebakova et al., 2009(Tchebakova et al., , 2011, which offset increased terrestrial carbon storage from boreal forest expansion upslope and northward into tundra. Persistent forest losses are more pronounced in the southern mountain taiga because of topographic conditions that result in distinct differences in climate and vegetation cover over short distances (Tchebakova et al., 2011).
Boreal postfire recruitment trajectories feed back to climate through modifications of ecosystem productivity and respiration (Reich et al., 2001), albedo (Randerson et al., 2006;Rogers et al., 2013), and evapotranspiration dynamics (Bond-Lamberty et al., 2009;Kang et al., 2006). These feedbacks are significant in the broader context of ecosystem resilience (Holling, 1973;Gunderson, 2000;Reyer et al., 2015) and the potential for re-organization to a new state as the result of disturbance (Johnstone et al., 2016). The objective of this study was to differentiate between categorical levels of stocking density of replacement vegetation that persist after wildfire disturbance, although we recognize these categories may contain multiple recruitment trajectories, particularly in the case of IR and AR sites. IR and AR sites will form the basis of elaborating which recruitment trajectories characterize the region.
Although forests frequently do not return to prefire characteristics (Reyer et al., 2015;Walker et al., 2017) and the potential for postfire shifts in plant functional type in boreal forests and other ecosystems is well-documented (Barrett et al., 2012;Scheffer et al., 2012;Tchebakova et al., 2009;Wolken et al., 2011;Zedler et al., 1983) incorporating such shifts into models of interactions between climate and vegetation dynamics has only begun to be explored (Koven, 2013, Quegan et al., 2011Zhu et al., 2015). It is important to incorporate such shifts as arising from meaningful ecosystem processes such as midperiod snow cover, as opposed to naive and often site-specific thresholding values (Koven, 2013;Zhu et al., 2015). We must also understand such processes in the context of legacy effects versus actual environmental changes (Bergeron et al., 2017, Houghton and Nassikas, 2018, Pongratz et al., 2014, the former of which were impossible to evaluate without a consistent and reliable burned area data product for the region. We restricted our analysis to pine stands as this was known to be the prefire vegetation type in most recruitment failure sites. Fire severity in pine stands may be higher than in other forest types because of their location on south-to west-facing slopes on well-drained sandy soils, and unlike larch stands, the fact that they burn throughout the growing season and during the late summer (Swetnam et al. 1996), when fires tend to be more severe (Turetsky et al., 2011). Pine is less constrained in terms of its dispersal dynamics than larch (Kharuk et al., 2010a), and is not complicated by masting events and the potential interaction of these with disturbances (Ascoli et al., 2015;Shi et al., 2000). In the Mongolian pine-steppe ecotone, postfire recruitment was lowest in pine and spruce and higher in larch and spruce stands (Otoda et al., 2013), possibly because larch (Larix sibirica) is more resilient (Gower and Richards, 1990) and tends to dominate in the forest-steppe environment (Dulamsuren et al., 2016), where pine is less well-suited. Summer water stress, in particular, is worse for pine than larch (Kharuk et al., 2010b).

Conclusion
Recruitment failure sites in southern Siberia are associated with lower greenness and vegetation moisture levels, higher fire severity, and lower snow index than intermediate or abundantly recruiting sites. These differences become more pronounced over time, likely as a reflection of feedback mechanisms that reinforce recruitment trajectories postfire. This research supports the broad literature on fire as a catalyst in broad ecosystem shifts, which may be moderated or exacerbated by disturbance characteristics and interactions, site-specific conditions, and climate/weather during the disturbance and recovery periods. Shifts between Scots pine and grass-dominated vegetation in the region may be responsive to winter snow cover in the intermediate assessment period 3-4 years postfire. The potential for fire to convert forest K. Barrett, et al. Remote Sensing of Environment 237 (2020) 111539 dominated stands to open grasslands will likely have a substantial impact on the net terrestrial carbon balance, albedo, and, subsequently, climate (Heimann and Reichstein, 2008;Tchebakova et al., 2011). It is important to re-create such processes explicitly in coupled climate and vegetation dynamics models to determine the presence and magnitude of feedback mechanisms. The feedbacks that promote or hinder tree recruitment may become established in the first five to ten years postfire (Johnstone et al., 2004), though our ability to detect such feedbacks using remotely sensed data may take longer, particularly if they are the result of lag effects or conditions that occur after the initial disturbance. Finally, incorporating information on postfire recruitment trajectories in the region with coupled vegetation and climate models will require an elaboration of all recruitment trajectories beyond abundant and poor recruitment, which will be the basis of our future work.

Author's responsibilities
Kirsten Barrett is the lead author and designed the analysis of field and remotely sensed data, Robert Baxter helped to design the field sampling method and contributed to the interpretation of remotely sensed data, Elena Kukavskaya oversaw field data collection and contributed to study design and data interpretation, Heiko Balzter contributed to the writing and editing of the manuscript, Evgeniy Shvetsov contributed to the description of remote sensing methods and data analysis, Ludmila Buryak was instrumental in identifying field sites and interpretation of the ecological significance of remotely sensed data.

Funding
This work was supported by the UK Natural Environment Research Council [grant number NE/N009495/1].

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.
Appendix A. This appendix contains the figures generated using remotely sensed data with the year of burn estimated as the earliest fire date on record Fig. 11. Time series of growing season NDVI (top) and NDMI (bottom) for AR, IR and RF sites, with trendlines superimposed for the postfire period. Values for AR sites increase postfire, whereas RF and IR sites either stabilize or decline post-fire. K. Barrett, et al. Remote Sensing of Environment 237 (2020) 111539 Table 2 Mean separability of AR and RF sites as a function of assessment period and season. The earliest separability of AR and RF sites occurs during the winter 3-4 years postfire.