Thirty years of forest productivity in a mountainous landscape: The Yin and Yang of topography

We measured light-related patterns of primary productivity within a topographically complex Oregon watershed over a 30-year period. Second-growth conifer densities were experimentally altered in 1981. Plots receiving at least 3434 MJ m − 2 over a 6-month growing season averaged 40% greater aboveground net primary productivity (ANPP) than those receiving less light ( p = 0.000). Unthinned stands potentially built enough LAI to compensate for low light, but risked mortality that exceeded resilience. The two light levels acted as basins of attraction for other physiological and ecological processes, including size – density relationships and limiting foliar nutrients. Initial (1981) LAI and the irradiation step (above or below 3434 MJ m − 2 ) explained 60% of variation in a 30-year ANPP. Irradiation within each light group did not affect ANPP. At high irradiation, foliar N/Ca and slope steepness (both negative) explained 58% of the variation in residuals from the initial models, while at low irradiation on north, east, and west aspects, 83% of residual variation was explained by foliar Mg (+), understory cover (+), and 30-year mortality ( − ). Light use efficiency (LUE) of fully stocked stands correlated with LAI and foliar N/K. Results suggest that understory influence on tree foliar N (+ or − ) enhances ANPP by regulating critical nutrient ratios. Mortality reduced or eliminated differences among thinning levels, which did


INTRODUCTION
The importance of forests in regulating global climate has been much discussed (e.g., Bastin et al., 2019;Domke et al., 2020).However, significant gaps exist in our understanding of the factors influencing terrestrial carbon sinks (Baldocchi & Penuelas, 2019).Remote imagery provides indispensable large-scale measurements of the carbon cycle and at least some factors that influence carbon dynamics (Schimel et al., 2019), but does not consistently correspond with data from field studies, particularly in complex terrain or disturbed forests (Berner et al., 2017;Girardin et al., 2016;Majasalmi et al., 2017).
The global network of eddy covariance sites provides crucial tools for documenting carbon fluxes (Baldocchi, 2020) but is error prone in complex terrain (Jocher et al., 2017;Madani et al., 2017).Dynamic global vegetation models have large uncertainties that require data to resolve, particularly local heterogeneity (Bonan & Doney, 2018;Hararuk et al., 2019;Rollinson et al., 2017).Reichstein et al. (2014) argued that a deeper understanding of the factors influencing terrestrial carbon fluxes will require localized on-the-ground monitoring and research that can be incorporated into models, a point echoed by others (e.g., Fisher et al., 2014;Girardin et al., 2016).Local effects-for example, topography, soils, and disturbance-are more important determinants of plant functional traits than macroecology (Bruelheide et al., 2018), especially in the complex topography of mountainous landscapes (Nicklen et al., 2019).Ma et al. (2017) found that including topographic variables improved estimates of forest biomass, and interaction between local topography and regional climate affects the growth and occurrence of individual species (Whitbeck et al., 2016).
Moreover, and particularly in 21st-century climates, understanding dynamics of disturbed and recovering landscapes is critical for the accuracy of predictive models.Natural and managed disturbances superimpose variations in stand stocking levels on topoedaphic variables.Wildfires, insect and pathogen infestations, and wind and ice storms take a toll (e.g., Coomes et al., 2012), and that is almost certainly going to increase.In mountains, wildfires are often of mixed severity, and vary according to topography (Perry et al., 2011).Given the significant increase in wildfires (Artés et al., 2019) and the increasing prevalence of density management to reduce fire hazard (Stephens et al., 2020), long-term studies documenting the productivity and resilience of partially stocked forests becomes an important component of modeling future carbon sinks (Pretzsch et al., 2019).
An additional layer is superimposed on topoedaphics and stand structure: plasticity.Flexibility in plant physiology and allocation is an important driver of carbon flux variation (Zhou et al., 2017), which raises questions about how plants adapt to environmental diversity in complex landscapes.The newly emerging field of functional biogeography extends the question of local effects to examine functional plasticity within a given species.While some aspects of intraspecies functional variability have been known for decades (e.g., aboveground carbohydrate allocation, specific leaf area), the set of rich possibilities has been little studied (Westerband et al., 2021).How species respond to changing climates will almost certainly depend on their functional diversity, which may include not only allocation priorities but also physiological adaptations and nutrient selectivity to support adaptive physiology.Knowledge of functional variability across complex landscapes gives a more nuanced picture of the factors that will shape future carbon sinks.
Herein we report 30 years of aboveground net primary productivity (ANPP) in conifer plantations within a topographically complex watershed of the western Cascade Mountains, Oregon, USA, exploring how interactions among irradiation, topography, nutrients, and experimentally manipulated stand density influence ANPP.We had four primary objectives: 1.To test how well initial LAI and variables related to topography (irradiation, aspect) predict 30-year aboveground productivity.2. To test how well variables measured later in stand life (e.g., mortality, foliar nutrients) improved models of 30-year productivity.We hypothesized that foliar N and mortality would be significant factors.3. To test the relation between light use efficiency (LUE) and foliar N. Following evidence in the literature, we hypothesized the two would correlate positively.4. To quantify how mortality affected long-term ANPP differences among stands with different initial stocking levels.We hypothesized that heavier mortality in more densely stocked stands would narrow differences.

Study area
The study was conducted in the HJ Andrews Experimental Forest, centered at 44.21 N, 122.41 W on the west slopes of the Cascades Mountains in Oregon, USA.The climate of the Experimental Forest is wet and fairly mild in winter and warm and dry in summer (Bierlamaier & McKee, 1989).From 1979 to 2015, daily air temperature averaged 9.0 C, and yearly precipitation averaged 2194 mm (HJA Primet), 80% of which fell during winter.Annual precipitation may exceed 3500 mm on higher ridges (Dyrness et al., 1974).The climate moisture index (CMI) for the general area ranges from 100 to 200 cm year −1 (Berner et al., 2017).Soils are formed from andesitic rocks in the Order Inceptisol, Suborder Udept, Soil Series Kenney gravely loam (UC Davis SoilWeb https://casoilresource.lawr.ucdavis.edu/soilweb/).They are highly porous and store 30-40 cm of water (Bierlamaier & McKee, 1989).Study plots are in the Western Hemlock Zone (Franklin & Dyrness, 1973); however, the dominant forest type prior to logging was old-growth Douglas-fir, with a subcanopy dominated by western hemlock (Tsuga heterophylla), western red cedar (Thuja plicata), and several hardwood species.Vine maple (Acer circinatum) is common on many plots, as is a shrub/fern flora similar to, and almost certainly sprouting legacies of, the preceding old-growth stands.Currently, second-growth Douglas-fir dominates all study sites, with minor admixtures of western hemlock, western red cedar, and various hardwoods.
The study was installed in 1981 in four plantations, two on each side of Lookout Creek, the main drainage of the Experimental Forest.Lookout Creek in the area of the study drains east-west, bounded on the north by SE facing Blue River Ridge and on the south by north facing Lookout Ridge.Topography on the spur ridges that hold the study sites is complex.On the north side of Lookout Creek, 19 of 23 plots were facing in a southeast to southwest quadrant (150 clockwise to 210 ), while on the south side 21 of 24 plots were facing in a east-north-west quadrant (140 counterclockwise to 260 ).Slopes ranged from 12% to 75% on the south side of the watershed, and 12% to 70% on the north side.Table 1 shows topoedaphic variables.
Plantations were established after clearcutting and broadcast burning.One was logged in 1958 and planted in 1963, one logged in 1966 and planted in 1967, one logged in 1959 and planted in 1960, and one logged in 1959 and planted in 1963.Presumably, the broadcast burning following logging destroyed all advanced regeneration.The amount and timing of volunteer tree establishment after burning is unknown, so the oldest trees may have dated from within a year of logging.Therefore, at the time the experiment was initiated in 1981, the maximum age of stands could have ranged from 15 to 23 years old, but later planting and continuing natural fill-in would have resulted in varying amounts of younger trees that grew into the overstory.Except for one small area of grand fir (Abies grandis), which apparently did not survive, all were planted with 2-year-old Douglas-fir (Pseudotsuga menziesii) seedlings at initial densities ranging from 1313 to 1413 seedlings ha −1 .By 1981, after mortality and fill-in by volunteers, the pre-thinning stocking density averaged 1279 and 1300 trees ha −1 in the two plantations on the north side of Lookout Creek, and 3412 and 2412 trees ha −1 in the two plantations on the south side.(Table 2 gives stocking levels expressed as relative densities.)Stocking of unthinned plots on the south side (mostly northerly facing) was highly skewed, mostly due to a few locales on one site (107) with exceptionally high stocking due to volunteers.Overstory composition in 1981 averaged 90% and 99% Douglas-fir on the south and north side of Lookout Creek respectively, with the balance split fairly evenly between western red cedar and western hemlock.Hardwoods (mostly maples) were abundant in the mid-canopy but did not reach the overstory.

Experimental design and measurements
The design was split-plot, with thinning levels as main plots and subplots consisting of four treatments: pruning, fertilization with slow release tabs, fertilization plus pruning, and control (no treatment).Velasquez-Martinez et al. (1992) describe subtreatments as well as analytical procedures for soil and foliar nutrients.Main plots were located randomly within each plantation, and subplot T A B L E 1 Topoedaphic variables for the study sites.
One plot was installed in the approximate center of each subtreatment to measure variables described below.These "measurement plots" were buffered from adjacent sub-and main plots by at least 10 m.One subtreatment plot within an unthinned main plot was subsequently dropped because of insufficient buffering, leaving 47 plots in total.Measurement plot size was based on number of trees rather than area.Fifty trees per plot was the target, but that was not always possible while maintaining adequate buffers, and the actual number ranged from 35 to 50.Measurement plot sizes averaged 0.02 ha (SD = 0.01) for unthinned plots, 0.07 ha (SD = 0.02) for moderately thinned plots, and 0.16 ha (SD = 0.06) for heavily thinned plots.
Aboveground biomass was calculated from periodic diameter measurements of each living tree applied to the allometric model of Jenkins (2004).Following standard practice, diameters were measured at 1.37 m.We estimated leaf area by applying the foliar biomass equations of Gholz et al. (1979) to the Douglas-fir SLA value of Berner and Law (2016), and GPP using the equation of Waring et al. (2016): where Pb is the belowground proportion of GPP.After Waring et al. (2016) and Berner et al. (2017), we set c = 0.5 and Pb = 0.46.We calculated GPP carbon as 0.5 × dry mass.We calculated absorbed photosynthetic irradiation (APAR) by applying Beers Law to LAI, with the k value for Douglas-fir = 0.48 (Ungs, 1981), restricting that to stands with average LAI between 2001 and 2010 above 4.To calculate LUE for the period 2001-2010, we used mean APAR for 2001-2010.
Soils for macronutrient, pH, rockiness, and anaerobic mineralizable N (AN) analyses were collected in 1987 to a mineral soil depth of 20 cm from three random points per plot.(Rockiness was distributed unevenly among thinning levels so we do not include it in productivity models.)Analytic procedures are described in Velazquez-Martinez and Perry (1997).Foliage for nutrient analyses was collected in late summers of 1987 and 1996 from five random trees per plot.Branches were collected from midcrown at two opposite points per tree.Analytic procedures were as described for 1987 collections by Velasquez-Martinez et al. (1992), except in 1996, we sampled only current year needles.We use 1996 values in productivity models because they were the midpoint of the study period.
We estimated microbial biomass C using our AN measurements and a conversion equation derived from data in Myrold et al. (1989), who sampled a range of forested sites in western Oregon.In their data with two outliers removed, a linear regression between AN and microbial biomass C (determined by chloroform fumigation) had r 2 = 0.96 (p = 0.004, F = 249), with a slope of 6 mg kg −1 microbial C per mg kg −1 of AN.AN is a well-known proxy for microbial biomass (Binkley & Hart, 1989).
Understory vegetation cover was estimated visually in 2003 at four random points per subplot, and recorded in three height classes: "low" = up to 0.5 m; "medium" = 0.5-1.5 m, and "tall" = 1.5 m to at least 10% of lower overstory tree crowns.In 2010, cover of major understory by species was visually estimated within 1 m radius subplots in one-half of the plots.
Growing season solar irradiation (in megajoules per square meter) was calculated on each subplot with a hemispherical viewshed algorithm (Fu & Rich, 2002) based on a 1-m resolution LiDAR-derived digital elevation model of the ground surface.The ESRI ArcMap Solar We assumed photosynthetically active radiation (PAR) was 50% of the total (e.g., Waring et al., 2016), acknowledging that the proportion is likely higher when, or where, indirect irradiation comprises a relatively large proportion of PAR (Mercado et al., 2009).Daily temperature records from an HJ Andrews weather station at a similar elevation as our study sites showed that, between 1993 and 2016, freeze-free days ranged from 133 to 209 days year −1 and averaged 170 days year −1 .At the same weather station, the number of days that high temperatures exceeded freezing ranged from 346 to 365 days year −1 , and averaged 358 days year −1 (http:// andlter.forestry.oregonstate.edu/data/abstract.aspx?dbcode=MS001), which implies that temperatures are suitable and water available for photosynthesis during a substantial portion of the year (Waring & Franklin, 1979), and that our assumed 6-month "growing period" is conservative.

Statistical analyses
We analyzed main and subplot effects on ANPP with a general linear model, not including covariates because those of interest (topography, irradiation) had a highly nonlinear relationship with ANPP (see Results).For ANPP and other variables of interest that showed non-linearities in scatterplot, we used the R packages randomForest (vesion 4.6-14), and randomForestExplainer (https://github.com/ModelOriented/randomForestExplainer) to identify important explanatory variables, and Ipred and rpart (bootstrapped regression trees [RTs]) to quantify values of important explanatory variables at tree nodes.
We tested residuals for autocorrelation with the Durbin-Watson statistic and normality of residuals with the Shapiro-Wilk test.Autocorrelation within individual variables was tested with Moran's I.All reported regressions were free of autocorrelation and had normally distributed residuals.
Thirty-year ANPP increased up to an initial (1981 post-thinning) leaf area of ~4 m 2 m −2 , and then plateaued (Figure 2).Significantly lower ANPP than expected given their 1981 post-thinning LAI which occurred on two plots with mortality greater than 80%, both unthinned stands on 70% or steeper slopes whose mortality appeared to occur in a pulsed, catastrophic event (an ice storm) rather than through self-thinning.These plots may or may not be included in subsequent analyses (we will indicate), although we emphasize that they represent significant data points regarding risk associated with high stocking densities in vulnerable locations.Biomass increased approximately linearly over the 30-year period, while LAI increased up to 2001 then plateaued (Figure 3).The level at which LAI plateaued differed significantly among thinning levels, indicating that light available to individuals determined when wood became the primary sink for ANPP.
A split-plot analysis of variance (not including the two plots with catastrophic mortality) showed that ANPP of thinning levels varied significantly, while subtreatments did not (Table 3).(Including catastrophic mortality did not alter the basic results.)Differences among thinning levels depended on irradiation and are discussed in a following section.

Interactions between irradiation and topography
The strongly nonlinear scatterplot of ANPP versus irradiation (Figure 4a), was shown by bootstrapped RT to split into above or below a growing season irradiation of 3434 MJ m −2 , corresponding mostly but not solely to the two sides of the watershed.On the north side (broadly southerly slopes) local topographic shading resulted in four plots with irradiation levels low enough that RT lumped them with plots on the shaded side of the watershed.All plots on the south side of the watershed were below 3434 MJ m −2 .Based on those results, we divided the data into two irradiation groups, high and low (or sunlit and shaded) with 3434 MJ m −2 as the boundary.
Over the 30-year period, plots on low irradiation sites averaged 214 Mg ha −1 , while those receiving high irradiation averaged 40% higher, 300 Mg ha −1 (p = 0.000, t test) (Figure 4b).When adjusted for irradiation-related differences in 1981 post-thinning LAI, the mean difference was reduced but still highly significant (217 and 295 Mg ha −1 , p = 0.000).Despite a large range within the low irradiation group (2270-3400 MJ m −2 over 6 months), the irradiation we measured played no role in ANPP of plots within that group ( p = 0.59), nor did it in the high irradiation group which, despite occupying approximately the same area, had a much smaller range (3434-3564 MJ m −2 ).
Two of the five foliar macronutrients differed between irradiation levels.Ca was higher in low than in high irradiation (0.56% and 0.45%, p = 0.003), while K had the opposite trend (0.84 in low, 0.92 in high, p = 0.009).Exchangeable Ca did not differ between the irradiation levels, while soil K was slightly higher on high irradiation sites.

Predicting long-term productivity from initial conditions
We produced four general linear models of 30-year ANPP versus a set of initial conditions (1981 post-thin LAI, and either irradiation or aspect).Two models differed in how the irradiation term was included, either as a linear or a step function (below or above 3434 MJ m −2 over a 6-month growing season).A third substituted aspect for irradiation, and the fourth substituted facing (north or south side of the watershed) for aspect.We used log (1981 LAI) in all.Figure 5 shows how the three variables relate.Notably, the relation between aspect and irradiation differed significantly between the two sides of the watershed, south aspects having the highest radiation on the south-facing side of the watershed, and lowest irradiation on the north-facing side.All models produced statistically significant results, with small differences in their predictive ability (Table 4).R 2 values ranged from 0.55 to 0.61, with only minor differences in the Akaike information criterion.Expectedly, 1981 LAI (log transformed) was highly significant in all.The strongest models included either the irradiation step function (above or below 3434 MJ m −2 ) or aspect.

Beyond initial: Modeling residuals
To investigate how additional variables (including those measured later in stand growth) influenced 30-year ANPP, we modeled residuals from the initial LAI + IrradiationStep model shown above (to avoid confusion we term these "residsI").We focused on irradiation because it is the most generalizable, though as we will show below, aspect could not be ignored.Inputs to general linear models included foliar nutrients measured in 1996, 30-year mortality, cover of medium and low height understory (cover of subcanopy trees was confounded with thinning treatment), soil pH, and anaerobic mineralizable N (shown in the literature to correlate with conifer growth; see Discussion).
Models that included 45 plots (all but the two with catastrophic mortality) performed poorly, either having non-normal residuals, very low explaining power, or both.Accordingly, we grouped the watershed by irradiance, using the split at 3434 MJ m −2 over a 6-month growing season.On sunlit, plots (≥3434 MJ m −2 ) residsI correlated negatively with both foliar N/Ca and slope steepness (Table 5).The negative correlation with N/Ca was particularly strong.N/Ca, in turn, correlated negatively with the cover of shrubs and ferns, increasing cover associated with N/Ca more favorable to ANPP (Figure 6).
On low irradiation plots (<3434 MJ m −2 ), we obtained a significant model only when restricted to those with broadly north, east, and west aspects (adjusted aspect <140), which excluded six plots on southerly aspects but with low irradiation due to local topographic shading.In order to get a model with normally distributed residuals it was also necessary to exclude the two plots with catastrophic mortality.The resultant model, which included 20 plots, explained 83% of the variation in residsI with three significant variables: foliar Mg (p = 0.000) and medium height understory contributing positively (p = 0.000), and the interaction between thinning level and mortality negatively (p = 0.000).

LUE (ANPP/APAR)
We calculated the commonly used LUE (ANPP/APAR) for plots with mean LAI above 4 m 2 m −2 for the period 2001-2010, which comprised stands absorbing at least 85% of PAR.Twenty-one plots fit that criterion, 15 on the south-facing side of the watershed (all unthinned or moderately thinned plots) and 6 on the north-facing side (5 unthinned and 1 moderately thinned).LUE followed a Gamma distribution with a mean (aboveground dry matter) of 1.0 g MJ −1 (range = 0.50-1.76).Expressed as GPP aboveground carbon, the mean was 1.85 g C MJ −1 .Contrary to our hypothesis, LUE correlated poorly with foliar N, which explained only 4% of the variation in LUE (p = 0.37).However, LUE correlated significantly and positively with foliar N/K (r 2 = 0.34, F = 9.2, p = 0.007, log-log) (Figure 7).The general linear model log(LUE) = c + LAI + logN/K explained 75% of the variation in LUE (F = 25.9, p = 0.000).In the regression and the GLM, we excluded one outlying plot (>2.5 SD from the mean) because a model including it did not have normally distributed residuals.

LUE on an area basis (ANPP/PAR)
We asked how efficiently stands utilized the irradiation falling on a given area (LUEarea), which allowed us to explore LUE across a range of stocking levels and avoid assumptions about foliar absorption.Within a light group, thinning had surprisingly little effect on the area-based efficiency of light gathering over the 30-year period.Thinning made no difference in the LUEarea of low-light plots, while at high light, only heavily thinned and unthinned plots differed (Table 6).Between light groups, only moderately thinned stands differed, those at low light having significantly lower LUEarea and higher mortality than those at high light (Figure 8).Moderately thinned stands at low light also had a lower initial (1981) relative density (RD) than at high light (Table 2), which over time could have accounted for the lower LUEarea.To test that possibility, we used a general linear model to separate the effects of initial RD from those of mortality.Fifty-six percent of variation in LUEarea of moderately thinned stands was explained by the interaction between irradiation and mortality ( p = 0.001).Initial RD did not enter the model.

Mortality and ANPP
Mortality's effect on ANPP varied with both irradiation and thinning level (Figure 9).On unthinned, sunny plots greater than 40% mortality sharply reduced ANPP, while on unthinned shaded plots ANPP increased with mortality up to 75%, and then decreased sharply on the three plots with mortality over 80%.On the former, survival correlated strongly with understory shrub/fern cover, mortality dropping from an average 60% for understory cover <10% to an average 27% for cover >30% (r 2 = 0.77, p = 0.01).Shrub/fern cover, in turn, correlated negatively with foliar N/Ca, higher cover being associated with N/Ca more favorable to ANPP (r 2 = 0.57, p = 0.05).

Size-density lines
Across all thinning levels, the size-density line differed significantly between the light groups (Figure 10).At any given stocking density, mean tree biomass (aboveground) was lower in low irradiation than high ( p = 0.000), while slopes did not differ.

DISCUSSION
Objective 1: How well can 30-year production be predicted solely from initial leaf area and variables related to topography (irradiation, local aspect, large-scale aspect)?
Models incorporating only initial leaf area and various topography/irradiation variables were highly significant, explaining from 55% to 61% of 30-year productivity.There was a large overlap between plots in the topography and irradiation groupings (only six of the 45 plots changed membership depending on partitioning scheme), and models were relatively insensitive to whether irradiation or topography was the independent variable.
Topography's effect on plant community types and diversity is well known (Moeslund et al., 2013), but, except for water balance, its effect on ecological processes where ABS is the absolute value.Plots on the south side of the drainage occur on a main ridge that faces north, while plots on the north side of the creek occur on a main ridge that faces south.Of the two S plots that don't fit the general S pattern, one is in a narrow creek bottom shaded on both the east and west sides, the other is on a steep west-facing spur ridge shaded on the east side.
T A B L E 4 General linear models of 30-year aboveground net primary productivity as a function of initial (1981) LAI and either aspect or irradiation excluded.has received less attention.On the HJ Andrews, old-growth forests attain their highest biomass on southerly aspects (Seidl et al., 2012), as we found for secondary forests, while in the drier environment of Colorado's highest forest biomass is on northerly aspects (Swetnam et al., 2017).Topography explained soil moisture in the Alps better than indices related to vegetation type (Raduła et al., 2018), and terrain explained a significant proportion of climate-related greening on the Tibetan Plateau (An et al., 2018;Teng et al., 2021).Although remote sensing of irradiation and ANPP in mountains has a long history (e.g., Peterson & Waring, 1994), we are not aware of previous ground research that directly relates ecological processes to irradiation patterns in mountains.

Model
RandomForest-ranked aspect as having greater mean relative importance than irradiation for ANPP.Aspect almost certainly has significant legacies related to irradiation that are not captured by the slice-of-time measure we made; for example, fire history and erosional patterns that influence soil processes.Aspect also integrates direct irradiation effects over decades, capturing yearly variability in cloud cover and growing season length that was unaccounted for in our one-year irradiation calculations.While we focused on light because of its greater generality, and for the most part, aspect and irradiation gave the F I G U R E 6 The relation between foliar N/Ca and residuals from the aboveground net primary productivity versus 1981 LAI + step irradiation model.Points are grouped by the percent cover of understory shrubs and ferns ("Medveg23").The 23% cover cutoff was determined by a regression tree model of the relation between residuals and understory cover.Foliar nutrients were measured in 1996.
same information, a few plots were hybrids, having both low irradiation and southerly slopes.For these, neither irradiation alone nor aspect alone accounted for patterns, adding yet another layer of landscape complexity.
Objective 2: How well do biologic and ecologic variables developing later in stand life explain residuals from the productivity-initial conditions equation?
We hypothesized that residuals would correlate significantly with mortality, and, drawing on results from an earlier study of these sites, foliar N, K, and Mg.These hypotheses were partially true, but overly simple because they ignored the effects of topographic variation within the watershed.Not surprisingly, variables (beyond initial) that influenced productivity depended on irradiation levels.Foliar Mg correlated positively with ANPP, but only on shaded sites.Given Mg's central role in light capture (review by de Bang et al., 2021), its importance on topographically shaded sites is not surprising.On shaded plots, gathering enough light to stay above the light compensation point is a top priority.Increasing the number of chloroplasts is one pathway for accomplishing that (Lichtenthaler & Babani, 2004), which would increase demand for both N and Mg.However, foliar N did not correlate with productivity on low light plots, nor did N/Mg, suggesting that foliage had ample N to meet the needs of chloroplasts and was limited by Mg.An imbalance between foliar N and Mg may result from relative availability.While N is supplied primarily by cycling, the majority of Mg and other nutrients come from rock weathering (Bormann & Likens, 2012), and rock weathering, in turn, ties strongly to biological activity (Bormann et al., 1998).On low irradiation sites, lack of energy may limit rock weathering, hence the supply of Mg relative to N. Excess N could then be stored or converted to other uses (Ripullone et al., 2003).On sunlight sites, the importance of foliar N was contingent on foliar Ca, ANPP attaining its highest values when N/Ca was between 2 and 2.3, and lowest when N/Ca was above 2.3.The low values of N/Ca were, in turn, associated with low understory cover, which may be due to the understory regulating foliar N, or signaling some unknown factor that influences N/Ca.We did not find any significant correlations between understory cover and site factors that we measured.We did find a strong unimodal relation between shrub/fern cover and foliar N (D. A. Perry and D. R. Oetter, unpublished manuscript).If shrubs/ferns affect N supply to trees, which is reasonable, understory competition for N would not necessarily decrease, and may increase overstory growth by regulating nutrient ratios.
Our results suggesting a role for understory should be extrapolated carefully, and likely depend on the particular mix of understory species.For example, one of the heavily thinned plots in our study developed a high cover of salal (Gaultheria shallon), and trees had high mortality (D. A. Perry and D. R. Oetter, personal observations).
However, the fact remains that understory cover and species composition can significantly influence overstory health and growth, and not always as expected.
The significance of foliar N/Ca on sunlit sites is consistent with Mainwaring et al. (2014), who found that Douglas-fir in the western Cascades did not respond to N fertilization if soil Ca to N ratios were too low (all ours were well above their threshold).Ca may be involved in growth in multiple ways (reviewed by de Bang et al., 2021;Wang et al., 2019).It is an essential component of photosystem II, and 20 photosynthesis-related proteins are associated with Ca.Perhaps most relevant to our high irradiation plots, Ca plays an essential role in dissipating excess energy, regulating electrical signaling that optimizes the balance between PSII quantum efficiency and excess energy dissipation through non-photochemical quenching (Białasek et al., 2017).
Our study stands in sharp contrast to others that have shown a strong link between ANPP and foliar N alone (e.g., Reich, 2012;Smith et al., 2002).However, although foliar N frequently correlates closely with primary productivity, there is no physiological reason that should always be so.Warren and Adams (2001) found that Pinus pinaster consistently over-invested in RUBISCO, leading to poor correlation between productivity and N.There are at least three possible reasons N alone did not explain ANPP in our study.First, our method of measuring foliar N differs from others'.We sampled from midcrown, which may not reflect changes in wholecrown N, especially if there are changes in leaf biomass not reflected in allometric equations.That seems likely.In our earlier study of these sites (Velasquez-Martinez et al., 1992), foliar N increased on fertilized plots only if they were also pruned, strong evidence that trees prioritize whole tree foliar retention over foliar N concentration.However, it is not clear why that would affect N but not the cations.Secondly, the western Cascades differ environmentally, ecologically, and geologically from studies showing strong N-only links to ANPP.As we discussed earlier, ample moisture and long-growing season produce relatively fast decomposition and N-cycling, and high cover of N-fixing shrubs following fire provide a history of large N inputs to soil (Youngberg & Wollum, 1976).These factors suggest that N supply may be sufficient to make it co-limiting, with physiologically essential cations input by rock weathering.Thirdly, we are not aware of other N-productivity studies that have documented the range of irradiation, which, as we found, would affect the demand for nutrients in addition to N. While all of these factors could be in play, our findings concerning the importance of cations in ANPP have clear physiological bases and, in the case of Ca, are supported by fertilization experiments in the western Cascades.It is not surprising that growth responds to more than a single element.The importance of nutrient ratios for tree growth is well known (Ingestad, 1979;Knecht & Göransson, 2004), and the DRISS system of nutrient evaluation is based on ratios (Bangroo et al., 2010).
The importance of Ca and concerns about its loss have been discussed in areas with base-poor soils (e.g., McLaughlin & Wimmer, 1999;Ouimet & Duchesne, 2005).The Andesites of our study are rich in cations relative to felsic bedrocks, however, abundance in rocks does not necessarily imply supply can keep up with that of N, especially in systems with rapid N cycling.The physiological demand for nutrient balances requires models that account for these supply-demand factors.

Light thresholds
Although the relation between irradiation and ANPP could be modeled linearly, two lines of evidence argue The relation between 30-year mortality and 30-year aboveground net primary productivity (ANPP), by irradiation and thinning level (HT, heavy thinning; MT, moderate thinning; UT, unthinned).
that a step model is the more accurate description.First, irradiation and ANPP were unrelated within each irradiation group.Second, the size density lines were different at a high level of probability, illustrating that light-related tree size differences crossed all diameter classes and stocking levels.
In our study, the light step reflected, in part at least, the risk of building sufficient leaf area to compensate for low light.That is clear in the scattergram of ANPP versus irradiation.Thinned stands on shaded plots did not build enough leaf area to compensate for low light.Unthinned, shaded stands potentially built enough leaf area to compensate for low light, but at the risk of mortality that exceeded resilience, resulting in wide variance and lower mean productivity than unthinned stands at high light.
Risk may not be the only factor contributing to a threshold light response.At the level of leaf physiology, the FvCB model (Farquhar et al., 1980) predicts an irradiation-related threshold transition between limiting physiological processes.Models that include costs associated with increasing leaf-level light compensation points, what Givnish (1988) terms the ecological compensation point, commonly show nonlinear responses to irradiation (Falster et al., 2018;Givnish, 1988).The "Kok Effect," inhibition of respiration by light (Kok, 1949;Lee et al., 2020), is possible on sunny plots.As Tcherkez et al. (2017) point out, "even minimal changes in leaf respiration may have a significant impact on plant C budget."Seasonal patterns in temperate conifer forests show that inhibition of respiration rises sharply with increasing radiation, then plateaus before eventually dropping off (Keenan et al., 2019), which is quite similar to the ANPP pattern we saw in space.However, if the Kok effect was a factor on sunny plots that should be reflected in higher LUE compared to shaded plots.We did not see that, but our calculation of LUE had too few shaded plots to make a definitive statement.
Trees in low light may have priorities other than growth.Huang et al. (2021) found that spruce trees under photosynthetic deprivation allocated C preferentially to storage at the expense of growth.Energy could also be devoted to relatively high leaf construction costs in the shade.The model by Kikuzawa (1991) and leaf economics analysis by Shipley et al. (2006) are relevant here.Kikuzawa's model predicts that leaf longevity correlates negatively with absorbed photosynthetic radiation (APAR), which, if true in our stands, would imply needles on shaded trees live longer.The analysis by Shipley et al. shows, in turn, that leaf longevity correlates F I G U R E 1 0 Log-log plot of mean tree biomass (aboveground) versus stocking density in 2010, by irradiation level.All thinning levels are included, excluding two plots with catastrophic mortality.Gray bands are 95% confidence intervals.High irradiation: adjusted r 2 = 0.89, SE = 0.13, F = 193.0;low irradiation: r 2 = 0.80, SE = 0.17, F = 99.8.From the Chow test, the null hypothesis that the two curves describe the same population is rejected (F = 13.9, p = 0.002).
positively with the proportion of leaf interior occupied by cell walls, which, because of Ca's importance in cell walls and membranes (Marschner, 1995), could explain the large Ca component in shaded plots of our site.As all our samples were current needles, leaf construction to extend life would have to begin in the first year, which seems reasonable.The hypothesis that the difference in foliar Ca between light levels stemmed from Ca supply (exchangeable) was rejected (t test p = 0.58; D. A. Perry and D. R. Oetter, unpublished data).

ANPP variation within irradiation groups not related to measured irradiation
Variation within the two ANPP groups on either side of 3434 MJ m −2 was relatively large and did not correlate with the irradiation we measured.That is not surprising for the high light group, which, despite encompassing approximately the same area as the low light group, had a much smaller range of irradiation.It may seem surprising that the low light group showed no response to its wide range of light, however, as discussed in preceding sections, models predict and field research confirms that under low light substantial fixed carbon is allocated to processes other than growth, hence little or no correlation with primary productivity is plausible.The primary driver of photosynthesis in topographically shaded locales may well be diffuse irradiation that we did not measure (e.g., beneath the canopy), which has been shown to contribute significantly to primary productivity (Bi et al., 2022).Xie and Li (2020) found that diffuse radiation remained constant across all aspects in a mountainous terrain, which would explain why stand-level diffuse irradiation would not correlate with our measured irradiation.
How does our primary productivity compare?
Accounting for difference in stand age, our mean yearly primary productivity was comparable to other studies in the western Cascades (Berner et al., 2017;Turner et al., 2007;Van Tuyl et al., 2005).Our GPP carbon was very close to the maximum reported from FLUXNET tower data for evergreen needle leaf forests (ENF) (Madani et al., 2017).
Objective 3: Does LUE correlate positively with foliar N?
LUE has been correlated with foliar N in several studies (Kergoat et al., 2008;Madani et al., 2014;Reich, 2012), but not all.Schwalm et al. (2006) found no connection between N and LUE at 11 FLUXNET sites in Canada.In our watershed, foliar N alone explained only 4% of the variation in LUE, while foliar N/K accounted for 34% and foliar N/K plus LAI explained 75%.As with ANPP, a nutrient ratio was more predictive of LUE than a single element (Ingestad, 1979;Knecht & Göransson, 2004).Sardans and Peñuelas (2021) reviewed the importance of K in plants.Like Ca, it has many functions, including supporting LUE and formation of Rubisco.Unlike the divalent cations, K is cycled much like N (Tripler et al., 2006), which affects its availability compared with Ca and Mg, and may explain why foliage contains excess N compared to the divalent cations, but not relative to monovalent K.In lab-grown seedlings, the optimal N/K ratio for coniferous plants is ~2, however N/K of most conifer forests ranges from ~1 to 1.2 (Knecht & Göransson, 2004).Ours ranged from 1 to 1.6, with the highest values on shaded plots.It was the high values on shaded plots that drove the positive N/K-LUE relationship, which confounds a strict causal relation for N/K.The physiological requirements of LUE are highly likely to differ between sunny and shaded sites, perhaps resulting in different N/K requirements.If true, our results would not reflect an N/K limitation so much as the relation between leaf physiology and light gradients.
LUE also correlated positively with LAI, likely reflecting a failure of Beer's Law to describe the light-gathering ability of densely stocked stands.Binkley et al. (2013) raised that concern, pointing out that Beer's law could not account for changes in crown structure that facilitated light capture.
The one outlying plot excluded from the LUE model is an interesting case study.That plot had the highest soil N in the study, among the highest mineralizable Ns, and the highest ANPP through 1996. However, between 2001and 2010, the period for which we calculated LUE, its productivity dropped dramatically and was the lowest in the study.That triggered us to look for a general relationship between fast early growth and lower later growth.For unthinned plots (but not thinned), growth from 2001 to 2010 correlated negatively with growth from 1981 to1996 (p = 0.08) (D. A. Perry and D. R. Oetter, unpublished manuscript).That decline did not correlate with mortality and appears to be a tortoise and hare phenomenon in which early fast growth cannot be maintained.Perhaps intense tree-tree competition within densely stocked stands is a lose-lose proposition that results in early resource limitations (more on this below).
How does our LUE compare?
Our mean LUE (ANPP/APAR) of stands with LAI above 4 m 2 m −2 , 1.0 g MJ −1 , was on the low end of the majority of studied tree plantations, which cluster between 1 and 1.5 g MJ −1 (Waring et al., 1998).Expressed as GPP carbon, 1.85 g MJ −1 , our LUE equaled FLUXNET's maximum for ENF.Albaugh et al. (2016) found that the aboveground LUE of nine plantations (distributed among tropical, temperate, and boreal sites) was a relatively uniform 1.51 g MJ −1 .Our lower value reflects at least two things.First, the Albaugh et al. data was for experimental plantations fertilized and watered to achieve optimum growth, while ours reflected real-world limitations.Second, whereas Albaugh et al. (2016) used productivity and irradiation in a single year, we used 10-year mean productivity and irradiation measured in 1 year, assuming irradiation in that year was representative of what was received over the 10-year period, thereby introducing an unknown error.Moreover, as Albaugh et al. (2016) point out, LUE based on a single year's growth, and LUE based on many year's accumulation are not the same.A major reason for that is mortality.A single year may or may not have mortality.Over the 10 years encompassed by our LUE, mortality ranged up to nearly 30%.Given these factors, the uniformity of the Albaugh et al finding (1.51 g MJ −1 ) implies that the combined effects of mortality and nutrient-water limitations reduced ANPP on our fully stocked stands by ~0.5 g MJ −1 year −1 .
Finally, some studies showing N to be a significant factor in LUE and, by extension, productivity, had no data from our region.Both Kergoat et al. (2008) and Madani et al. (2014) used eddy-covariance installations, which included no sites from the western Cascades (Madani et al.).Kergot et al. included one site in a comparable environment (western Canada), but it had no N data.
Objective 4: Does mortality reduce long-term ANPP differences among stands with different initial stocking levels?We hypothesized that, because of heavier mortality in more densely stocked stands, 30-year ANPP would not differ among thinning levels Our hypothesis was partially true but too simple.Mortality was certainly one factor in reducing or eliminating long-term productivity differences among initial stocking levels.At low irradiation, ANPP did not differ statistically among thinning levels.At high irradiation, unthinned stands maintained higher productivity than heavily thinned, but did not differ from moderately thinned.However, while much mortality related to stand density, not all did, and the complex effect of mortality on long-term productivity varied with topography and irradiation as well as stand density.Among our unthinned plots, mortality exceeding resilience was determined by landscape position (steep slopes on north aspects).Stocking density and consequent small tree size probably played a role, but we cannot make a definitive statement about that.Among moderately thinned plots, mortality was higher and ANPP lower on shaded than sunlit plots, even though shaded plots had lower initial stocking.
Does light level influence tree-tree interactions?
It is not clear why moderately thinned stands should have higher mortality and lower resilience at low than at high light.Initial stocking levels do not explain it; as pointed out above, moderately thinned stands at low light had lower initial stocking than at high light, which, if anything, should have lowered mortality.Heavy competition for light may have cost the winners as well as the losers.If diffuse light was an important component of irradiation on shaded plots, tree-tree competition may be quite different than on sunlit sites, more three dimensional and symmetric than top down and asymmetric.

Tortoises and hares: What you see may not be what you get
Site occupation is a commonly assumed primary factor influencing NPP (Yuan et al., 2018).However, less than full site occupancy is becoming more common in global forests due to natural and human disturbance.Natural wildfires often burn with mixed severity that leaves a mosaic of forest cover commonly related to topography (Perry et al., 2011).Full stocking today may leave forests vulnerable to extensive dieback in future climates (Jump et al., 2017), and management often involves thinning to increase resilience to fire, drought, and insects (e.g., Hessburg et al., 2021).
What does a view that incorporates long-term forest health imply for carbon-sink potential?Our study shows that, depending on topography, high stand density does not necessarily increase productivity over the long term, and may decrease it if mortality exceeds system resilience (cf.Anderegg et al., 2020).Because of these factors, and especially in reduced light, relatively low stocking density may well be the tortoise to high density's hare.
Less than full site occupation by trees provides balance with other values, such as biodiversity and stand health (Franklin et al., 2018).In the long view, with planning that considers the landscapes of risk and resource supply, and a management focus on maintaining resilience (Hessburg et al., 2021;Prichard et al., 2021), these values may be attained with minimal long-term sacrifice of carbon sinks.The environment-specific benefits, costs, and risks of forest density across complex landscapes are critical to predicting the future and warrant further study (Franklin et al., 2018), including better understanding of the potential for a web of connections via shared mycorrhizal fungi and how that affects stand-level behavior (e.g., Birch et al., 2021;Zhu et al., 2024).

Allometry
Since we used the same allometric equation for aboveground biomass, all we can say with certainty about aboveground growth is what we measured-bole growth, though the correspondence between our data and other studies' data gives confidence in the allometric equations and parameters we used.Nonetheless, there is little doubt that allometric equations are overly coarse.In one of the few studies to document local variation in allometry of mature conifers, Espinosa Bancalri and Perry (1987) found that three adjacent Douglas-fir stands of similar age and minor differences in topography allocated biomass differently.
Trees may allocate less to stems and more to foliage at low light levels, but allocation changes are small compared to those in leaf morphology (Poorter et al., 2012).In 1987, we calculated leaf area on a subset of trees using the pipe model for Douglas-fir (Waring et al., 1982) and found no difference among sites or thinning levels (Waring et al., unpublished data).Even if there were allocation differences we did not detect because wood constitutes the majority of aboveground biomass (>90% in the allometric equations we used), any differences in allocation would be unlikely to balance out the ANPP differences between irradiation levels, though they could well influence processes.

Size-density curves
We found that the intercept, but not the slope, differed between irradiation levels, which is consistent with Lonsdale and Watkinson (1982), who found that deep shade lowered the size-density intercept of the perennial grass, Lolium perenne.The intercept is likely determined by the whole tree light compensation point (WTLCP), which, to remain positive in reduced light requires reducing respiration, hence less bole wood relative to gathering surface.Perry's (1984) model of physiological/allometric factors in the size-density line predicted that the intercept depends on ratio of crown to bole wood, essentially a proxy for the respiration-photosynthesis balance.

Implications
A drop from the ocean: who would have thought this infinitely little too much Robinson Jeffers.(from his poem, Science) Topography has the potential to limit stimulating effects of global warming on carbon sinks, an effect accentuated by low sun angle.Moreover, potential effects emerge from complex ecosystem dynamics and are not linear: there are points on the landscape where the surrounding topography, underlying edaphics, and vegetation structure interact to produce nonlinear dynamics in which small changes in driving factors can have large consequences.If our models smooth over thresholds we leave ourselves open to surprise.
The specific patterns we found are almost certainly unique to the topography and latitude of the study sites.However, the fact such thresholds exist on our sites suggests they likely do elsewhere, although the details will vary with topography, latitude, soils, resource limitations, and species-specific functional plasticity.However, all mountains create patterns of light and shadow with sharp boundaries.All have complex edaphics in one form or another.All mountain forests adapt structurally and biochemically to the mosaic of environments they encounter, adaptations that, among other things, may require shifts in leaf economics and supporting nutrients.

More than mountains
Our findings have implications for high latitude forests whether considered mountainous or not.Several studies have found that neither higher atmospheric CO 2 nor faster N cycling due to warmer soils increase the growth of boreal tree species (Girardin et al., 2016).Various factors might account for that, including limiting photosynthetic energy and nutrient constraints on adaptability.Along with temperature and its effect on growing season, sunlight limits productivity at high latitudes (Running et al., 2004), andZhang et al. (2020) showed that light limits response to longer Autumn growing seasons above 40 N. Total irradiation drops sharply at latitudes north of ~35 (FLUXNET2015 cited in Baldocchi & Penuelas, 2019), and north of 60 virtually all irradiation levels recorded by FLUXNET2015 were below the 3434 MJ m −2 threshold we found.The details of thresholds may well change with species and adaptive pressures, but their existence and the adaptive responses they elicit are highly likely.Mountainous terrain accentuates light limitations, but at high latitudes even relatively low topographic relief creates aspect-related ecological signatures (Whitbeck et al., 2016).

Models
Our findings question the suitability of parametric models commonly used to describe nature's patterns.Others have raised that concern.Breiman (2001) provides a statistician's view: "Usually, simple parametric models imposed on data generated by complex systems … result in a loss of accuracy and information."In our data, linear regression gave statistically significant fits between irradiation and ANPP, but was blind to a stair step that was obvious in scatterplot and captured by RT.As growth emerges from interactions among a large set of processes that respond in complex ways to light and other driving variables, and ecosystems are far from thermodynamic equilibrium, it is not surprising that nonlinearities would emerge (Perry, 1995).Loss of ecologically relevant information is often the price of broad generalization and may or may not have implications for future responses to climate warming.When lost information is about non-linear threshold changes, it may well matter (e.g., the following section).More than prediction is at stake.Management that recognizes and adapts to complex system behavior requires data and analyses that describe complexity (Puettmann et al., 2012).

Geoengineering
We echo others who have raised concerns about large-scale geoengineering projects (Board Ocean Studies, 2015;McCormack et al., 2016;Russell et al., 2012;Zarnetske et al., 2021).Light-related thresholds increase the probability of unexpected responses to reduced sunlight, which has been suggested as a geoengineering option to combat global warming.Ecological changes out of proportion to lowered irradiation are a distinct possibility, including sharp reductions in terrestrial carbon sinks and alterations in feedbacks from land to climate.
First two principal axes showing the relation between aboveground net primary productivity and site variables.The first two axes explained 50.5% of variation in listed site variables.Adjaspect, 180 − ABS (180 − aspect); MinN, anaerobic mineralizable N (correlates highly with microbial biomass).
The relation between 1981 post-thinning LAI and 30-year aboveground net primary productivity (ANPP), by thinning level (HT, heavy thinning; MT, moderate thinning; UT, unthinned).The two unthinned plots with 1981 ANPP below 100 Mg ha −1 suffered greater than 70% mortality (most of which occurred in a single ice storm) and were on slopes >60%.

F
I G U R E 3 (a) Change over time in aboveground biomass and (b) change over time in LAI, by thinning level (HT, heavy thinning; MT, moderate thinning; UT, no thinning).Shaded areas represent 95% CIs.Year 0 = 1981 post thinning.See text for 1981 post-thinning densities.T A B L E 3 General linear model of 30-year aboveground net primary productivity.

F
I G U R E 4 (a) The relation between growing season irradiation and 30-year aboveground net primary productivity (ANPP), by thinning level (HT, heavy thinning; MT, moderate thinning; UT, unthinned).(b) Boxplot comparing 30-year ANPP between irradiation groups.The height of the box shows the interquartile range (IQR), which is the middle 50% of the data.The midline is the median, where the position of the median line shows skewness; if in the middle of the box there is no skew, while above (below) the middle denotes a left (right) skew.The whiskers show the range of plots 1.5 above or below the IQR.Points outside the whiskers are outliers.Comparing the upper and lower graphs, note that outliers are unthinned plots.The irradiation groups differ at p = 0.000 (t test).See text for 1981 post-thinning densities.

F
I G U R E 8 (a) The relation between LUEarea (mean yearly aboveground net primary productivity [ANPP]/PAR) and irradiation, by thinning level (HT, heavy thinning; MT, moderate thinning; UT, unthinned).Mean ANPP calculated for the entire 30-year period.(b) The relation between thinning level and 30-year mortality, by irradiation level.Two plots with catastrophic mortality were not included.See Figure 4 caption for explanation of the boxplot features.
Downloaded from https://esajournals.onlinelibrary.wiley.com/doi/10.1002/ecs2.4865,Wiley Online Library on [16/07/2024].See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions)on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License Factors explaining residuals from the aboveground net primary productivity = logLAI81 + Irradstep models of Table4.Irradiation cutoff: 3434 MJ m −2 per 6-month growing season.Plots with <3434 MJ m −2 were restricted to those with N, NE or NW aspects (adjusted aspect <140).Both models are free of autocorrelation (Durbin-Watson test) and have normally distributed residuals (Shapiro-Wilk test). Note: