Variability in Mixed Conifer Spatial Structure Changes Understory Light Environments

In fire-adapted conifer forests of the Western U.S., changing land use has led to increased forest densities and fuel conditions partly responsible for increasing the extent of high-severity wildfires in the region. In response, land managers often use mechanical thinning treatments to reduce fuels and increase overstory structural complexity, which can help improve stand resilience and restore complex spatial patterns that once characterized these stands. The outcomes of these treatments can vary greatly, resulting in a large gradient in aggregation of residual overstory trees. However, there is limited information on how a range of spatial outcomes from restoration treatments can influence structural complexity and tree regeneration dynamics in mixed conifer stands. In this study, we model understory light levels across a range of forest density in a stem-mapped dry mixed conifer forest and apply this model to simulated stem maps that are similar in residual basal area yet vary in degree of spatial complexity. We found that light availability was best modeled by residual stand density index and that consideration of forest structure at multiple spatial scales is important for predicting light availability. Second, we found that restoration treatments differing in spatial pattern may differ markedly in their achievement of objectives such as density reduction, maintenance of horizontal and tree size complexity, and creation of microsite conditions favorable to shade-intolerant species, with several notable tradeoffs. These conditions in turn have cascading effects on regeneration dynamics, treatment longevity, fire behavior, and resilience to disturbances. In our study, treatments with high aggregation of residual trees best balanced multiple objectives typically used in ponderosa pine and dry mixed conifer forests. Simulation studies that consider a wide range of possible spatial patterns can complement field studies and provide predictions of the impacts of mechanical treatments on a large range of potential ecological effects.


Introduction
In many fire-adapted conifer forests of the Western U.S., changes in land use activities and disturbance regimes such as logging, grazing, and fire suppression have led to increased forest densities and altered overstory composition [1]. In ponderosa pine-dominated (Pinus ponderosa Livermore, California, USA; ±10 mm horizontal). From the reference locations, a Pentax PCS-515 laser total station (TI Asahi Co., Saitama, Japan; ±3 mm horizontal) was geolocated and referenced to true north and used to establish a square 20 × 20 m set of grid cells across the site. Grid cells were then used for locating all overstory trees and sampling locations. The perimeter of each cell created by the grid was framed with measuring tapes and used to record the location, species, diameter at breast height (dbh) for over 15,000 live trees exceeding 1 cm dbh. No tree was ever estimated more than 10 m from the framing tapes, resulting in an estimated RMSE of ±15 cm of the true horizontal location. Because treatment objectives in dry conifer forests of the region often include explicit goals for increasing spa

Modeling Light Availability
To estimate understory light availability across a density gradient, we selected 56 plots from a previously established set of approximately 400 plots from a related study (Tinkham, unpublished data). We selected plots within the dry mixed conifer portion of the stand (Figure 2) and stratified a gradient of basal area (ranging from 3 to 41 m 2 ha −1 as measured from a 15 m fixed radius). In July 2017, we captured canopy photographs in the center of each plot using a 24 megapixel mirrorless DSLR camera with 180 • hemispherical lens mounted on a self-leveling gyroscope and tripod (Compact OMount, Regent Instruments, Inc,. Quebec, Canada). We took photographs at dawn or dusk to produce a clear contrast between canopy and sky, with the camera mounted at 1.5 m in height and oriented toward true north. We classified each image as canopy and sky pixels using threshold grey values in Gap Light Analyzer [45]. We transferred classified photos to Winscanopy software (Regent Instruments, Quebec, Canada) and modeled understory solar radiation as photosynthetically active mean growing season photon flux density (PPFD; µmol m −2 s −1 ). To parameterize Winscanopy, we modeled understory radiation during the growing season (May-September) at an average elevation of 2810 m, and we assumed a solar constant of 1353 µmol m −2 s −1 [46], an atmospheric transmissivity of 0.622 (mean value for Colorado Springs, CO between May and September) [47], and a standard overcast sky diffuse light model [48]. To facilitate comparison with other studies, we converted all PPFD measurements to proportion of full sun (FS) by dividing radiation values by the maximum radiation modeled for the site of 577.4 µmol m −2 s −1 . Example simulation output demonstrating overstory tree spatial patterns for all treatments described in Table 1. Yellow shaded portions represent south-facing aspects dominated by dry-mixed conifer cover included in the analysis. Single trees (10%); 2-4 trees (30%); 5-9 tree (35%); 10-15 trees (15%); 16+ trees (10%) Figure 2. Example simulation output demonstrating overstory tree spatial patterns for all treatments described in Table 1. Yellow shaded portions represent south-facing aspects dominated by dry-mixed conifer cover included in the analysis. Understory light availability depends on a number of forest structural attributes such as canopy cover, tree density, tree height and vertical stratification, canopy architecture, and tree spatial pattern [31,32,49]. To model how forest structure impacts understory light availability, we used a model selection approach to describe how understory light availability, as modeled through canopy photography, varied as a function of forest structural attributes at a range of spatial scales. We considered five basic overstory structural metrics including basal area (BA; m 2 ha −1 ), tree density (all stems > 12.7 cm, D 13 ; stems ha −1 ), large tree density (stems > 25.4 cm, D 25 ; stems ha −1 ), stand density index (SDI; 25 cm stem equivalents ha −1 ) using the Shaw [50] method, and neighborhood index (NI). Neighborhood index and similar metrics have been used as proxies for tree competition (e.g., [51]), and the metric is based on tree size (BA), but trees are weighted as a function of distance from the observation point (plot center), implicitly including a spatial component in this metric. We calculated the neighborhood index for each plot center using the following equation: where NI r is the neighborhood index of a plot with radius r, BA i is the basal area of the ith tree within radius r, d i is the distance (in m) of the ith tree from the plot center for all n trees within radius r, and A r is the area (in ha) of the plot. For all plots, we calculated each of the five overstory metrics at multiple scales, ranging from 3 to 30 m in radius (in 2 m increments). Next, we modeled light availability as a function of overstory structural metrics at a range of scales through a model selection procedure to determine what metric and scales best describe understory light availability. Structural metrics compared included BA, D 13 , D 25 , SDI, NI. We modeled light availability using each structural metric measured at a range of scales, and included both single-and dual-scale models. For single-scale linear models, light availability took the following form: where m s is each structural metric (BA, D 13 , D 25 , SDI, NI) at scale s. Because light availability is typically non-linearly related to tree density [52,53] and because our calculation of light availability is a proportion bounded from 0 to 1, we modeled light availability using a beta regression with a natural log link using the betareg package in R with parameter optimization via the BFGS gradient method [54]. Spatial patterns such as edge effects may influence light availability, so we also generated dual-scale models that considered each structural metric as measured at two scales. Dual-scale linear models took the following form: where m s is each structural metric as measured at scale s, m s is each structural metric at scale s for all scales where s < s . Thus, these models included each forest metric at a small scale, s, and a larger scale s . After combining five structural metrics (BA, D 13 , D 25 , SDI, NI), two model forms (single-and dual-scale), and scales between 2 to 30 m (every 2 m), we generated 75 single-scale models (5 structural metrics × 15 scales), and 525 dual-scale models (5 metrics × 105 combinations of s and s ). To select the combination of scales and metrics that best predicted understory light availability, we calculated and compared Akaike information criterion (AIC) for all 600 models [55].

Simulating Restoration Treatments
To assess how simulated treatments varying in spatial complexity changed the distribution of light availability, we generated stem maps simulating a range of treatment outcomes and applied the best model of light availability obtained in the above analysis. To simulate treatments varying in spatial pattern, we used the untreated stem map to simulate treatments that were similar in overall residual basal area, yet differing in tree spatial patterns including thinning from below, thinning overstory trees randomly, and restoration treatments varying in tree group size distributions ( Figure 2, Table 1). Previous studies identify groups of trees as those with the potential for trees to develop interlocking crowns and continuous canopy, which in several studies of ponderosa pine is often represented as those with less than approximately 6 m spacing (e.g., [56,57]). In this study, we also used an intertree distance of 6 m based on the approximate 2.6-3.4 m crown radius range predicted from allometric equations linking tree diameter to crown radius [57] when trees reach the larger range of sizes noted at the site of 40-60 cm, and in the upper range of sizes noted in historical reconstructions of the region [2]. In addition, we use the term "group size" to refer to the number of individual trees in a group, rather than the areal coverage of tree groups which is more commonly referred to as "tree patch size" (e.g., [20]).
For each mechanical thinning simulation, we used a target basal area of 11.5 m 2 ha −1 (Table 1). These values were chosen based on guidelines recommended for the region [19] and a silvicultural prescription developed for the site. We conducted 10 simulations for each stochastic thinning treatment. Example simulations for each treatment are shown in Figure 2. To simulate random thinning, we thinned all small trees (< 12.7 cm, < 5 in.) to a density of 40 ha −1 and then randomly selected larger trees for removal to reach the basal area target. To simulate thinning from below, we thinned starting with the smallest trees and proceeded to larger trees until the basal area target was reached.
To simulate more spatially complex restoration treatments, we first thinned small trees randomly to 40 ha −1 , but overstory trees were thinned with varying degrees of spatial aggregation to a range of target tree group sizes (Table 1). We conducted 10 restoration simulations each with low, moderate, and high target levels of tree aggregation (Table 1), following Tinkham et al. [29]. We developed an algorithm in R that thins trees to the desired basal area, while reaching tree group size targets. The algorithm selects a tree in the stem map to initiate each group. The tree is randomly selected from a weighted probability of dbh in cm raised to the 1.6 power to simulate anchoring trees off of large existing trees (sensu [58,59]). Trees with boles within 6 m of the initial tree bole are added to the group, and proximal trees are added iteratively until a target group size is reached. The intertree distance used here represents a 3 m crown radius delineated to represent the potential for interlocking and adjacent crowns as trees mature.
When a tree group is formed, stems within 6 m of any group members are removed, and another group is initiated until the target stand basal area and target group size distributions are reached. See https://github.com/jbcannon/ICO_thinning_simulation (accessed on 12 November 2019) for example program and Figure S1 for example animation.

Assessment of Structural Complexity and Light Availability
To assess how simulated restoration treatments achieved desired conditions of decreasing density while maintaining horizontal and tree size complexity, we calculated basal area, stem density, quadratic mean diameter, and tree group size distributions. We used Shannon's equitability (E H ) to assess changes in horizontal complexity. This metric is often used to assess evenness in species assemblages but can be applied to other aspects of forest structure to assess complexity [60]. We determined tree density in 20 × 20 m grid cells and binned density into five classes (0-100, 100-200, 200-300, 300-400, and > 400 trees ha −1 ). We then calculated E H using proportions of the stand in each density class using the following equation: where p i is the proportion of each stand in density class i, resulting in an estimate of horizontal complexity ranging from 0 to 1, with 1 representing equal representation among the five density classes [61,62]. To assess changes in tree size complexity, we calculated a stand density index ratio using the Reinke [63] and Shaw [50] methods, a proxy for tree size complexity [64,65] using the following equation: where, SC is size complexity, D i is the diameter of the ith tree, TPH i is the number of trees represented by the ith tree, Dq is the stand-level quadratic mean diameter, and TPH is the stand-level stem density [65]. Note that we subtract from 1 the equation presented in Shaw [65] so that size complexity approaches 0 for even-aged stands, and increases for uneven-aged stands as size variability increases. We also present diameter distributions of treatment simulations to assess changes in tree size complexity.
To infer how simulated restoration treatments alter light availability, variability, and high light conditions, we applied the best predictive model of light availability from above to each stem map to generate raster data predicting light availability. A small edge buffer (6 m) was applied to the stem map to eliminate edge effects. Summary statistics were calculated to compare treatment effects on mean light availability (mean FS), light availability coefficient of variation (FS CV), and the 75th percentile light availability.

Modeling Overstory Structural Effects on Light Availability
All of the top 10 best performing models included either stand density index or basal area as the metric capable of explaining the most variation in light availability and had pseudo−r 2 values ranging from 0.4217 to 0.3512 ( Table 2). The top three models (with ∆AIC < 2) consistently included a smaller scale of 10 m, and a larger scale between 22 and 24 m ( Table 2). Although the models differ in the structural metric used (SDI versus BA), these metrics are closely correlated (r > 0.988) at all spatial scales. The best model of light availability included stand density index measured at 10 and 22 m (Table 3), which we used to predict light availability in subsequent simulation analysis. This model suggested that in dry mixed conifer stands, predicted light availability varies from 20% to 100% full sun, and that light availability increases rapidly when stand density index as measured at 10 m radius is reduced to < 300 to 500 ha −1 (Figure 3).

Effects of Treatment Simulations on Forest Structure and Complexity
Treatment simulations resulted in stands with a constant basal area yet varying in other aspects of stand spatial structure such as density, tree group size distribution, horizontal complexity, and tree size complexity (Table 4). Simulations reduced basal area from an initial 20.0 m 2 ha −1 to 11.5-11.6 m 2 ha −1 in agreement with simulation constraints (Table 4). However, overall stem density varied among simulations. Simulated treatments reduced average tree density between 57% and 83%, with thinning from below resulting in the largest reduction (501 to 86 ha −1 ), random thinning resulting in the smallest reduction (501 to 214 ha −1 ), and aggregated treatments with moderate reductions (from 501 to approximately 195 ha −1 , Table 4). We found that changes in tree group size distributions followed simulation parameters. Thinning from below the concentrated basal area in smaller tree group size classes, whereas random thinning led to greater evenness among tree group size classes (Table 4). Consistent with simulation constraints, we found that among the aggregated treatments, the low aggregation treatment exhibited the highest proportion of trees in smaller groups (1-4 trees per group), while the highly aggregated dry mixed conifer stands exhibited the greatest proportion of trees in larger groups (≥ 5 trees per group; Table 4). Simulations varied little in tree species composition and generally reduced species abundances proportionally to their initial abundance ( Figure 4A). However, one exception was that thinning from below essentially eliminated species primarily represented by smaller diameter individuals such as quaking aspen ( Figure 4A, Figure 5C).

Effects of Treatment Simulations on Forest Structure and Complexity
Treatment simulations resulted in stands with a constant basal area yet varying in other aspects of stand spatial structure such as density, tree group size distribution, horizontal complexity, and tree size complexity (Table 4). Simulations reduced basal area from an initial 20.0 m 2 ha −1 to 11.5-11.6 m 2 ha −1 in agreement with simulation constraints (Table 4). However, overall stem density varied among simulations. Simulated treatments reduced average tree density between 57% and 83%, with thinning from below resulting in the largest reduction (501 to 86 ha −1 ), random thinning resulting in the smallest reduction (501 to 214 ha −1 ), and aggregated treatments with moderate reductions (from 501 to approximately 195 ha −1 , Table 4). We found that changes in tree group size distributions followed simulation parameters. Thinning from below the concentrated basal area in smaller tree group size classes, whereas random thinning led to greater evenness among tree group size classes (Table 4). Consistent with simulation constraints, we found that among the aggregated treatments, the low aggregation treatment exhibited the highest proportion of trees in smaller groups (1-4 trees per group), while the highly aggregated dry mixed conifer stands exhibited the greatest proportion of trees in larger groups (≥ 5 trees per group; Table 4). Simulations varied little in tree species composition and generally reduced species abundances proportionally to their initial abundance ( Figure 4A). However, one exception was that thinning from below essentially eliminated species primarily represented by smaller diameter individuals such as quaking aspen ( Figure 4A, Figure 5C).   Table 4. Results of treatment simulations on metrics of stand-scale forest structure. Metrics include basal area (BA), stem density, quadratic mean diameter (D q ), and tree group size distribution (single, 2-4 trees, 5-9 trees, 10-15 trees, 16+ trees). Horizontal complexity (E H , Equation 4, [60]) and tree size complexity as measured by SDI ratio (Equation 5, [62]). For both E H and size complexity, values of 1 indicate greater horizontal or size complexity, respectively, while values closer to 0 indicate less complexity. Values shown are mean (sd) of ten stochastic simulations. Untreated stands had only one replicate and thus no variability metrics shown. Note some proportions do not sum to unity due to rounding.

Metric
Units Untreated Random Thin Below

Effects of Treatment Simulations on Light Availability and Variability
Increases in mean light availability were relatively similar across treatments, but treatments exhibited differences in within-stand variability in light levels (Table 5, Figure 6A). Mean light availability increased from 0.592 FS in the untreated stand to between 0.676 FS in the thin from below treatment to 0.665 in the low-aggregation treatment ( Table 5). Variability of light levels paralleled results for within-stand horizontal complexity, where untreated stands had light levels primarily concentrated in the 0.40 to 0.60 FS range, whereas thinning from below concentrated a large area in the 0.70 to 1.0 FS range ( Figure 6A). Among the aggregated treatments, the highly aggregated stand Treatments varied in their ability to maintain or enhance aspects of horizontal and tree size complexity with thinning from below resulting in the largest reductions in complexity, while random and highly aggregated treatments best maintained both aspects of complexity (Table 4). Both untreated stands and thinning from below had the lowest horizontal complexity, with both stands having density equitability < 0.70 (Table 4) and > 60% of the total area with stem densities < 100 ha −1 or > 400 ha −1 , respectively ( Figure 4B). Random thinning had the highest horizontal complexity (0.93), and aggregated treatments had moderate horizontal complexity (from 0.78 to 0.91). Among the aggregated treatments, the highly aggregated treatment preserved the most area in the highest and lowest density classes (< 100 and > 400 ha −1 ), whereas the low-and moderately aggregated treatments resulted in a greater abundance of moderately dense classes (200-400 ha −1 , Figure 4).
As expected, diameter distributions differed dramatically among untreated and treated stands with large reductions in density of stems < 30 cm in all simulated treatments ( Figure 5) with the result that all treatments reduced tree size complexity (Table 4). However, the thin from below treatment especially reduced size complexity from an initial value of 0.11 in untreated stands to 0.01, whereas all other treatments only slightly reduced size complexity to between 0.07 and 0.08 (Table 4). Random thinning treatments retained more mid-sized stems (10-30 cm dbh) relative to all other treatments ( Figure 5B). Thinning from below eliminated most stems < 30 cm ( Figure 5C). Aggregated treatments were relatively similar in tree size complexity (ranging from 0.7 to 0.8). However, the high aggregation led to greater retention of mid-sized stems relative to the low-and moderately aggregated treatments and was thus most similar to the random thinning treatment ( Figure 5D-F).

Effects of Treatment Simulations on Light Availability and Variability
Increases in mean light availability were relatively similar across treatments, but treatments exhibited differences in within-stand variability in light levels (Table 5, Figure 6A). Mean light availability increased from 0.592 FS in the untreated stand to between 0.676 FS in the thin from below treatment to 0.665 in the low-aggregation treatment ( Table 5). Variability of light levels paralleled results for within-stand horizontal complexity, where untreated stands had light levels primarily concentrated in the 0.40 to 0.60 FS range, whereas thinning from below concentrated a large area in the 0.70 to 1.0 FS range ( Figure 6A). Among the aggregated treatments, the highly aggregated stand had the greatest variability in light availability (Table 5). Highly aggregated stands had more total area concentrated in the lowest (0.50) and highest (>0.80) light ranges compared to the other aggregated treatments. Conversely, the low-and moderate-aggregation treatments had a greater proportion of area concentrated in moderate light levels (0.50-0.80, Figure 6A). Changes in high-light availability increased in all treatments, with the 75th percentile FS increasing from 0.689 in untreated stands to between 0.745 to 0.778 (Table 5). Among the treated stands, 75 th percentile FS was highest in the high-thin from below (0.778) and high-aggregation treatment (0.773), lowest in the low-aggregation treatment (0.745), and moderate in the random thinning and moderate-aggregation treatment (0.755) ( Table 5).

Discussion
This study highlights how residual overstory spatial patterns have cascading impacts in mixed conifer forests such as altering the light environment and components of overstory structure such density, composition, and tree size complexity. The use of multi-scale models in predicting understory light transmittance from different forest structural metrics reveals an impactful spatial dependence. The 10 best performing models from our study highlight the importance of considering forest structure at two scales in heterogeneous forests. The larger scales (22-24 m) correspond approximately to the dominant overstory height in these stands, indicating that in spite of the high level of heterogeneity of forest structure in dry mixed conifer forests, light availability is in part driven by interactions between canopy height and latitude as has been found in temperate forests [49]. This spatial scale provides a logical approximation of the interaction between vegetation height and solar angle for representing the attenuation of light, which is approximated by basal area or stand density. Some studies have shown that a single radius (approximating vegetation height) can be used to model light availability when stands are homogeneous [66]. Conversely, in our study, a smaller scale (10 m radius) was also important to represent the highly variable fine-scale local forest density found within our study site. The natural clustering pattern found in most ponderosa pine sites [13,21] indicates that light availability in these forests is determined by both site-level density of trees, as

Discussion
This study highlights how residual overstory spatial patterns have cascading impacts in mixed conifer forests such as altering the light environment and components of overstory structure such density, composition, and tree size complexity. The use of multi-scale models in predicting understory light transmittance from different forest structural metrics reveals an impactful spatial dependence. The 10 best performing models from our study highlight the importance of considering forest structure at two scales in heterogeneous forests. The larger scales (22-24 m) correspond approximately to the dominant overstory height in these stands, indicating that in spite of the high level of heterogeneity of forest structure in dry mixed conifer forests, light availability is in part driven by interactions between canopy height and latitude as has been found in temperate forests [49]. This spatial scale provides a logical approximation of the interaction between vegetation height and solar angle for representing the attenuation of light, which is approximated by basal area or stand density. Some studies have shown that a single radius (approximating vegetation height) can be used to model light availability when stands are homogeneous [66]. Conversely, in our study, a smaller scale (10 m radius) was also important to represent the highly variable fine-scale local forest density found within our study site. The natural clustering pattern found in most ponderosa pine sites [13,21] indicates that light availability in these forests is determined by both site-level density of trees, as well as fine scale differences in tree locations. Future modeling efforts could investigate the use of both distance-weighted metrics like the neighborhood index along with directional effects, as a tree along the solar incidence path should have a significantly stronger effect on light attenuation.
Our study demonstrated that while similar levels of density reduction can result in similar mean light availability, the range and variability of this resource differs due to the differences in overstory spatial pattern ( Figure 6B). Increasing the degree of spatial aggregation within a treatment tended to increase variability in understory light availability. This is supported by both Martens et al. [67] and Battaglia et al. [31] who found that light variability increases when residual canopy structure is aggregated for both piñon pine (Pinus edulis Engelm.) and longleaf pine (Pinus palustris Mill.) forests, respectively. Our low-aggregation simulations extend this finding and indicate that when mechanical treatments decrease structural complexity, light variability can also be decreased relative to untreated stands. Such results highlight that a large range of ecological outcomes may occur across gradients of spatial pattern.
We found a non-linear regression best explained the relationship between basal area and light availability in dry mixed conifer stands (Figures 3 and 6B), which has important implications for achieving ecological and management objectives related to structural variability and treatment longevity. A study by Vyse et al. [38] demonstrated greater survival of ponderosa pine seedlings over more shade-tolerant species when light availability exceeded approximately 65% full sun (crossover point irradiance (sensu [36])), and found that this crossover point occurred at 5.0 m 2 ha −1 of basal area compared to an inflection point corresponding to between 10 and 20 m 2 ha −1 in our study. Optimal light environments for ponderosa pine regeneration following traditional silvicultural regeneration treatments have been reported as 9-27.5 m diameter group-selection openings [68] or shelterwood treatments with 6-14 m 2 ha −1 residual basal area (or not less than 4 m 2 ha −1 [69]) in comparison with more continuous cover from alternative treatments [70]. Moreover, aggregated overstory retention is generally recommended to improve seedling growth of ponderosa pine [70], where openings provide adequate resource availability for ponderosa which can otherwise be significantly negatively affected by competition from saplings and larger canopy trees [33]. Yet in some cases, fuel hazard reduction treatments fall short of achieving the requisite high levels of overstory aggregation [20,28,30], and instead, result in more uniform overstory retention which may favor shade-tolerant species such as Douglas-fir [38]. Uniform overstory retention may thereby reduce the longevity of treatments with an objective to reduce ladder fuels [41]. Although exact light intensities and gap structures necessary to favor shade-intolerant species will differ by region, forest type, and species under consideration, the general finding that large basal area reductions are necessary to favor ponderosa pine indicate that careful consideration of the placement, size, and frequency of overstory gaps (i.e., fine-scale spatial pattern) is critical to increasing the proportion of stands that favor ponderosa pine over more shade-tolerant species.
Increasing heterogeneity in overstory and light conditions gives rise to cascading effects on many abiotic and biotic parameters creating variable microsite conditions that may impact plant diversity, ponderosa pine regeneration, and stand resilience. In our study, the high aggregation treatments generated the largest area of high light conditions, with the highest 75th percentile light availability among all treatments. The relatively high-light environments within low-density stands or openings have been linked to high understory productivity, as well as high herbaceous diversity in stands dominated by ponderosa pine [18,71]. Furthermore, high microclimatic variability characterizes many mixed conifer forests [72] and has been implicated in fostering high understory herbaceous diversity [73]. Ponderosa pine regeneration will be driven, in part, by microenvironments with biophysical conditions that support favorable moisture availability over both space and time [74][75][76][77][78][79]. Light environments that have been shown to favor ponderosa density and survival include open areas of few to no trees, and areas with moderate to low canopy cover, especially from intermediate canopy trees (e.g., 6-30 cm dbh) [33], where temperature extremes are alleviated, moisture is retained, and vegetation competition is moderated relative to resource requirements of ponderosa. In our study, the highly aggregated treatment simulation led to high overstory structural complexity and variable light environment that may better support ponderosa pine regeneration. Previous studies report that light environments defined by canopy cover and structure may drive fine-scale biotic conditions relevant to regeneration such as litterfall distribution and depths, vegetation type and structure, rainfall interception and throughfall, temperature and wind patterns, seedfall patterns, seed-and tree-predatory wildlife habitat, and moisture and nutrient uptake [39]. Variability in light conditions is, therefore, indicative of variable resource quantity and distribution [32,80] for tree regeneration. Variation in light and corresponding resource environments can increase stand forest resilience by providing a diverse range of microenvironments (spatially and temporally) to support complex regeneration dynamics [81]. In addition, variable microsite habitats may support or limit different life-stages of trees, and lead to facilitative interactions for young trees and competitive interactions for older individuals [75,82,83].
Beyond ponderosa pine, many pine-dominated systems, especially those shaped by fire, are characterized by high levels of heterogeneity that can impact resource availability and stand regeneration processes. Old-growth longleaf pine systems exhibit patchy regeneration driven by intra-specific competition, fire, and hurricane mortality [84]. Gap size and position drive light availability and regeneration success in longleaf stands [85,86], and aggregated harvest patterns have been shown to increase resource availability at both local and stand scales [35]. Although our study found relatively similar mean-light levels at the stand scale across treatments examined, we found large variability in light availability following aggregated retention (e.g., more area with high light, and more area with low light). Similarly, in mature red pine (Pinus resinosa Aiton) stands, both light availability and variability increased with aggregated harvests [87]. Further, within-gap variation in light availability in moisture led to differences in photosynthetic rates among three northern Pinus species, indicating the potential for gap portioning in these systems [88]. In moisture limited Scots pine (Pinus sylvestris L.) systems, successful regeneration is often associated with small gaps that provide moderate shade and improve summer moisture conditions [89], highlighting that fine-scale microclimate variability can impact stand-dynamics.
Direct field measurements of changing light conditions are useful for understanding changes in light availability and variability under relatively simple mechanical prescriptions such as individual or group selection studies [31]. However, given the range of target conditions and outcomes of restoration in ponderosa pine and mixed conifer forests, simulation studies such as the one presented here provide a framework for understanding the impacts of a large range of restoration outcomes on an important abiotic gradient. Nevertheless, there are several limitations to use of treatment simulations. To focus analysis and interpretation on spatial patterns of overstory retention, our simulation model was based strictly on retention of specific overstory patterns and size structure and did not explicitly consider overstory composition, although results may differ if stands contain a large component of spatially aggregated species such as quaking aspen. Furthermore, our simple spatial simulation model was based on specific targets for retaining a distribution of tree group sizes. Although this model successfully simulated fine-scale spatial structure (Table 4), it did not include larger-scale considerations such as creation of large openings or retention of groups of trees greater than 20 individuals. Future work can improve upon existing efforts by realistically simulating typical tree marking practices such as anchoring from large trees, identifying and expanding existing gaps, and explicit incorporation of larger openings, or by comparing changes in light variability following actual treatments where spatial patterns were specified in the prescription. Fine-scale simulations can be integrated with coarser-scale simulation models under development [90] to address nuances of treatment implementation and outcomes at multiple scales. Furthermore, our study was opportunistically conducted in a stand where extensive stem mapping had taken place, and used an empirical model of light availability based on stand structure, which may vary in other ecosystems and locations that vary greatly in latitude, stand composition, and tree height [49].

Conclusions
In ponderosa pine and dry mixed conifer forests, management goals for restoration typically emphasize structural changes to increase resistance and resilience to disturbances. Some of these objectives include (1) reducing overstory density, (2) increasing horizontal complexity, (3) maintaining tree size complexity, and (4) reducing the potential for regeneration by shade-tolerant species, such as Douglas-fir, that can act as ladder fuels and reduce treatment longevity [91]. Our results suggest that even with similar reductions in basal area, treatment prescriptions that differ in spatial pattern may differ markedly in their congruence with each of these goals, exhibiting important trade-offs among objectives related to tree density, structural complexity, and light environment. Thinning from below resulted in the largest reduction in overstory density, yet led to the largest reductions in structural complexity due to removal of all small stems. Conversely, random thinning maintained much of the existing horizontal and tree size complexity, yet resulted in modest density reductions and provided less with high light availability. Because the aggregation simulation was designed to anchor the formation of tree groups around large trees, the low-aggregation treatment was similar to an even-spaced thinning treatment and was dominated by single trees and small groups of trees, and was thus relatively uniform, resulting in relatively low horizontal complexity, light variability, and the least area with high-light availability. Conversely, the high-aggregation treatment after anchoring on larger trees retained many of the smaller surrounding trees creating spatial complexity while also maintaining tree size complexity. Thus, the high aggregation treatment balanced several objectives resulting in moderate density, relatively high horizontal and size complexity, and among the highest light variability and area with high-light availability. Treatments intentionally designed to create a wide distribution of tree group sizes may best meet restoration objectives and improve resilience in ponderosa pine forests.