Climate has a larger effect than stand basal area on wood density in Pinus ponderosa var. scopulorum in the southwestern USA

Abstract

Key message

Stand basal area of ponderosa pine ( Pinus ponderosa var. scopulorum Engelm.) in the US Southwest has little effect on the density of the wood produced, but climatic fluctuations have a strong effect. Wood density increases during drought, particularly if the drought occurs in late winter/early spring. Future droughts, as are predicted to increase in the US Southwest, may lead to production of smaller radial increments of higher density wood in ponderosa pine.

Context

Forest restoration treatments in the US Southwest are generating large quantities of small-diameter logs. Due to negative perceptions about ponderosa pine wood quality, this material is often seen as a “waste disposal” problem rather than a high-value resource.

Aims

Our objective was to understand more about variation in southwestern US ponderosa pine wood density, an important indicator of wood quality. Specifically, we investigated the effect of stand basal area on wood density, and the effect of annual and quarterly climatic variation on wood density.

Methods

We collected samples from 54 trees grown at six different basal area levels from a replicated stand density experiment. Pith-to-bark strips were used in an X-ray densitometer to obtain annual density and growth measurements from 1919 to the present.

Results

Stand density had a strong effect on growth rate, but little effect on wood density. However, climatic variation did influence wood density, which increased in drought years before quickly returning to pre-drought levels.

Conclusion

Stand basal area is not a good indicator of wood density for foresters planning to utilize material from timber harvests in the southwestern USA. Future droughts, as are predicted to increase in the region, will likely reduce wood volume production but may increase wood density in ponderosa pine.

1 Introduction

Overstocked forests in the southwestern USA have prompted land managers to respond with landscape-scale mechanical thinning treatments designed to reduce the threat from catastrophic wildfires and improve forest health (Covington et al. 1997; Waltz et al. 2014; Kalies and Yocom Kent 2016). These treatments produce large volumes of mainly small-diameter (< 40.6 cm) ponderosa pine (Pinus ponderosa var. scopulorum Engelm.) stems and woody by-products (Hampton et al. 2008; Lowell et al. 2008). The markets available, including pallet stock, molding, post and pole, and clean or dirty chips (Lucas and Kim 2016), rarely cover the cost of timber harvest and transportation (Hjerpe and Kim 2008; Lucas et al. 2017). Additionally, the highly variable nature of the harvested stems (particularly in regard to stem form and diameter), and the often low density of ponderosa pine wood may limit its potential for use in higher-value wood products. Thus, a triple threat is formed: land managers struggle to pay for restoration treatments without state or federal subsidies, wood “waste” accumulates in the forests in the form of burn piles, and meanwhile, it is necessary to import wood products (State Import Data 2018).

Increased utilization of the forest resource can begin to solve these problems, but a major barrier is the lack of knowledge of the fundamental properties of ponderosa pine grown in the southwestern USA. Density is an important wood property that is often used as an indicator of wood mechanical properties (Kretschmann 2010), to quantify potential yields of biomass products (Taeroe et al. 2015), and to estimate carbon storage (Flores and Coomes 2011). Wood density is highly variable at multiple scales; species-level oven-dry density values in US conifers (green volume basis) range from 290 to 680 kg m⁻³ with ponderosa pine averaging 380 kg m⁻³ (Miles and Smith 2009). Within a tree, density values vary both within rings and from pith-to-bark due to the contrast between corewood (juvenile wood) and outerwood (mature wood; Burdon et al. 2004). For example, average density of an annual ring in Scots pine (oven-dry weight and green volume) ranged from 274 to 697 kg m⁻³ (Auty et al. 2014), a range greater than the species-level means of all US softwoods. A greater understanding of this variation in southwestern ponderosa pine is critical for selecting appropriate end-uses for the material, thus increasing potential revenue generated from restoration treatments and possibly making them more economically viable.

The goal of this study was to understand sources of wood density variation in Northern Arizona’s ponderosa pine resource. Outside of the changes from earlywood to latewood, the next largest source of variation is from pith-to-bark within an individual tree, due to changing hydraulic and mechanical needs as the tree ages, described as the “typical radial pattern” (Lachenbruch et al. 2011). For example, young trees tend to produce earlywood with narrow tracheids that are more resistant to embolism than large tracheids (Lachenbruch et al. 2011). Because of the tradeoff between water transport and wood density (Sperry et al. 2006), this causes a tree’s earlywood (EW) near the pith to be of higher density than earlywood near the bark. Whether or not mean ring density declines in the same fashion depends on the ratio of earlywood to latewood (LW); an increase in LW proportion can counter these EW density trends and cause a pith-to-bark increase in average ring density. Ring width generally decreases from pith to bark due to the geometric constraint of adding consecutive layers of wood to an enlarging stem and thus ring area (i.e., basal area increment) is a more reliable representation of tree growth than ring width (Gartner et al. 2002).

Stand density and competition for resources may influence wood density, but this effect depends on whether the species has an abrupt or gradual transition from EW to LW. In Sitka spruce (Picea sitchensis [Bong.] Carr.) (Gardiner et al. 2011) and Norway spruce (Picea abies L.) (Dutilleul et al. 1998), both gradual-transition species, wood density was negatively correlated with stand growth rate in spacing trials. This suggests that in these species, high-density stands may produce trees with superior structural wood quality. In abrupt-transition species (“hard pines,” such as ponderosa pine), results are less conclusive. Some studies have shown a negative correlation between growth rate and wood density in loblolly pine (Pinus taeda L.) (Jordan et al. 2008) and radiata pine (Pinus radiata D. Don) (Nicholls and Wright 1976; Bannister and Vine 1981). Meanwhile, no such correlation has been found in loblolly pine (Megraw 1985) and ponderosa pine (Myers 1960; Voorhies 1969).

The southwestern USA is subject to periodic droughts that have the effect of reducing crown growth (Adams and Kolb 2005) which can increase latewood proportion and therefore ring density (Larson 1969). Additionally, turgor pressure is the mechanism that drives cell expansion (Hsiao 1973; Rathgeber et al. 2016; Rodriguez-Zaccaro and Groover 2019). Therefore, in addition to indirectly affecting ring density through suppression of crown growth, drought could directly increase ring density by limiting the expansion of earlywood cells (Hsiao 1973; Rodriguez-Zaccaro and Groover 2019). Ponderosa pine growth is most responsive to the previous winter’s precipitation (Adams and Kolb 2005; Kerhoulas et al. 2013), suggesting that drought could cause a narrower earlywood band and increase the overall ring density. Many dendroclimatic studies have been conducted in the US Southwest, but they have typically focused on the effects of climate on annual ring width, while few have examined climatic effects on wood density components or ring density profiles. One exception from northern Mexico is Pompa-García and Venegas-González (2016), who found a strong positive correlation between winter precipitation and maximum latewood density in Cooper pine (Pinus cooperi Blanco).

Many factors can influence the wood density of ponderosa pine, including low-frequency variation associated with the typical radial pattern and stand density, and high-frequency variation superimposed by management actions and yearly climatic fluctuations. Here, we focus on understanding the importance of stand density, management history, and historical climate variability as influences on the radial profile of wood density in southwestern ponderosa pine. Specifically, we address the following questions:

• Does stand density affect ring density, latewood proportion, or other wood density components?

• Do climatic variables affect intra-ring wood density components? If so, are some seasons more influential than others?

2 Materials and methods

2.1 Study site

The study was located at Taylor Woods (35° 16′ 11″ N, 11° 44′ 30″ W), a replicated stand density “levels-of-growing-stock” experiment just outside of Flagstaff, AZ, USA. Taylor Woods is a naturally regenerated even-aged stand of ponderosa pine predominantly originating in 1919 (Schubert 1971). The site consists of a ponderosa pine overstory with scattered patches of New Mexico locust (Robinia neomexicana A. Gray) and an understory of Arizona fescue (Festuca arizonica Vasey). Slopes are less than 4% (Bailey 2008), and elevation averages 2266 m. Soils are productive for the region: relatively deep, well-drained Typic argiboroll over fractured bedrock (Meurisse 1971). The site index (base age 100) is 22.3 m (Bailey 2008). The area experiences a bimodal pattern of precipitation, with peaks in the winter months (November–March) and during the summer monsoon (July–August). Mean annual precipitation from 1919 to 2017 was 566.8 mm (Fig. 1) and mean annual temperature for the same time period was 6.3 °C (National Centers for Environmental Information, National Oceanic and Atmospheric Administration).

Fig. 1

Total yearly precipitation (trendline shown in blue) from 1918 to 2017 at Flagstaff Pulliam Airport, with inset showing 30-year normal monthly precipitation (NOAA). Figure produced in R with some limited post-processing in Inkscape

The experiment was established in 1962, when trees were approximately 43 years old. At this time, treatment units were thinned from approximately 47.9 m² ha⁻¹ to a basal area determined by their respective growing stock level (GSL). The GSL is the basal area that the residual stand will have when the mean tree diameter is 25.4 cm (Myers 1967). Thus, once stands reach a mean diameter of 25.4 cm, GSL is synonymous with basal area. Three treatment units were established for each GSL of 6.9, 13.8, 18.4, 23.0, 27.5, and 34.4 m² ha⁻¹ and have been subsequently thinned approximately every 10 years to maintain the target stand density (Fig. 2, Table 1). The most recent thinning occurred in 2017 and provided the material for this study. All thinnings were preceded by stand inventory/marking and favored the retention of dominant/codominant trees and the removal of trees with mistletoe and/or porcupine damage, poor form, excessive limbiness, and poor vigor (Ronco et al. 1985). Taylor Woods has an extensive history of research, and more information about stand history, stocking levels, and other related GSL studies can be found in Myers (1967), Schubert (1971), Ronco et al. (1985), Bailey (2008), and Uzoh and Oliver (2008).

Fig. 2

Basal area growth at Taylor Woods since 1962. Corresponding horizontal dashed lines show the desired basal area for the GSL level, while the solid lines show the observed values at 5-year intervals. Thinnings occurred in 1962, 1972, 1982, 1992, 2003, and 2017. Figure produced in R with some limited post-processing in Inkscape

Table 1 Summary statistics for GSL treatment and study trees. The treatment mean column shows GSL-wide average DBH, and the sample tree mean DBH column shows the average DBH of the nine trees sampled from each treatment. For treatment mean and sample tree mean DBH, standard deviation is given in parentheses

2.2 Plot/tree measurements and sample collection

Before the recent thinning in 2017, three 0.04-ha fixed-radius plots were randomly installed in each of the 18 treatment units at Taylor Woods. Diameter at breast height (DBH; 1.37 m above the ground level) was measured for all trees in these subplots to determine stand basal area. An inventory list for all trees scheduled for removal at Taylor Woods site was then used to randomly select three trees for destructive sampling from each of the 18 treatment units, for a total of 54 trees with 9 trees from each GSL treatment. Before felling, tree total height and the base of the live crown were measured. After felling, 2.54-cm-thick cross-sectional disks were collected every 2.4 m from ground level to a height of 7.32 m, with an additional sample taken at breast height. In total, 267 disks from 54 trees were collected.

2.3 X-ray strip processing and testing

To produce X-ray densitometry samples from the disks, a pith-to-bark strip was cut from the north side of each sample. The north side was chosen to avoid compression wood developing as a result of the region’s prevailing winds arising from the west/southwest. Strips were repeatedly soaked in acetone until the solution ran clear (typically two soaks). This was done to remove extractives, which can influence wood density but do not contribute to wood mechanical properties (Panshin and de Zeeuw 1980; Eberhardt and Samuelson 2015). Additionally, air-dry density of each sample was measured to aid in calibrating the densitometer. The radial strips were then cut to approximately 5 mm in the tangential direction, mounted on hardwood strips, and the assembly was cut to approximately 2.3 mm in the longitudinal direction. The samples were conditioned at 20 °C and 29% relative humidity to bring them to a testing moisture content of around 6%. The samples were tested on a Quintek QTRS-01X Tree Ring Scanner (Quintek Measurement Systems, Knoxville, TN) with a step size of 25 μm and the X-ray beam passing through the sample on the transverse face (Jacquin et al. 2017).

After obtaining the raw data from the densitometer, the next step was to delineate annual ring boundaries, cross-date samples, and determine the EW to LW transition. Rings were initially delineated at the latewood-earlywood boundary using the threshold-based ring boundary assignment given by the QMS Tree Ring Scanner software, with the threshold defined as the average density of each sample. The initial ring delineation required correction due to missing rings, false rings, and other unusual density patterns. Thus, each ring boundary was corrected by visual cross-dating, and ring boundaries were statistically cross-validated using the R package dplR (Bunn 2010). The cross-validating helped find errors in visual cross-dating and gave a quantitative method to screen out questionable samples. A total of 22 scans had a mean inter-series correlation below 0.35 (Adams and Kolb 2005) and were removed from the study, leaving us with 245 samples from 53 trees. The first ten rings from each sample were excluded because of high ring curvature near the pith and higher incidence of compression wood. Finally, the latewood boundary was assigned to the point in the ring where density reached 80% of the difference between the minimum and maximum values (Lundqvist et al. 2018). Data were summarized to produce response variables of basal area increment (cm²), ring density (kg m⁻³), LW proportion, average EW density (kg m⁻³), average LW density (kg m⁻³), and maximum LW density (kg m⁻³).

2.4 Climate data and chronology development

Total precipitation (PRCP, mm) and monthly average temperatures (TAVG, °C) were obtained from the Global Summaries of the Month, produced by the National Centers for Environmental Information (NCEI, National Oceanic and Atmospheric Administration [NOAA]). Data from the Flagstaff Pulliam Airport weather station were used except in years 1941 through 1947, which were missing. For these, data from the Fort Valley Experimental Forest weather station (approximately 11 km away) were substituted. Palmer drought severity index (PDSI) was downloaded from the National Climatic Data Center (NCDC, NOAA), using Arizona Division 2. PDSI is an index of drought severity calculated from precipitation, temperature, and potential evapotranspiration (Palmer 1965).

For each of the six response variables, mean chronologies of the de-trended series were produced. Cubic splines (using a 50% frequency response cutoff with 0.667 series length) were used to de-trend each series, removing low-frequency variation such as the typical radial profile while preserving high-frequency variation due to climatic fluctuation (Cook 1981, 1985). A dimensionless index was calculated (with mean of approximately one) using the ratio between the observed value and the detrended series. Autoregressive models were then fit to the indices to remove autocorrelation, a step known as “prewhitening” (Cook 1985; Bunn 2008). Next, chronologies were developed by averaging the indices using their biweight robust mean, which is a mean value produced by assigning higher weights to observations closer to the arithmetic mean, in order to reduce the effects of outliers and enhance the common signal (Cook and Holmes 1984; Cook 1985). The R package dplR (Bunn 2008) was used to accomplish all the steps in chronology development.

2.5 Data analyses

Two approaches were used to test for an effect of GSL: whole-sample averages and models that included covariates representing annual variability. The whole-sample averages included years 1963 to 2016 (all years after the initial thinning). Mixed-effects models were fit with a random effect for tree and fixed-effects of sample height, GSL, the sample height × GSL interaction, and the 5-year pre-1962 average. If type III ANOVA tables showed a significant effect of GSL (alpha = 0.05), Tukey-adjusted pairwise comparisons were investigated. The annual variability models had a similar structure but included terms for year of ring formation and PDSI. Autocorrelation was modeled with a first-order continuous autoregressive term. Again, type III ANOVA tables were used to test for significance of GSL and followed up with pairwise comparisons.

The R package treeclim (Zang and Biondi 2015) was used to assess the influence of climate on the chronologies. Response functions were calculated for the correlation between density components and quarterly climatic variables. These response functions are designed to deal with multicollinearity in the predictors and are calculated by regressing the response against principal components of the climate data to produce a “response coefficient” (Zang and Biondi 2013). Response coefficients have a possible range of values between − 1 and 1 and indicate strength of correlations in a similar manner to Pearson’s R values. To improve estimates of confidence intervals, stationary bootstrapping with 1000 resamples was used. The influence of climatic variables was investigated by quarter for the water year and at an annual scale; the first quarter of the water year was defined as October–December of the previous year, the second quarter as January–March of the current year, the third quarter as April–June, and the fourth quarter as July–September. If the confidence interval produced by the bootstrapped response function did not overlap zero, the effect of the climate variable was considered significant. To allow for differing responses among the GSLs, GSLs were grouped into low (6.9 and 13.8 m² ha⁻¹), mid (18.4 and 23 m² ha⁻¹), and high (27.5 and 34.4 m² ha⁻¹) levels.

3 Results

Higher GSLs were associated with an increase in both stand basal area and trees per hectare, and a decrease in mean tree diameter (Table 1). Of the 16,844 rings analyzed, the mean values for LW proportion, ring density, EW density, LW density, and maximum LW density were 12%, 446 kg m⁻³, 414 kg m⁻³, 715 kg m⁻³, and 763 kg m⁻³, respectively. The density components over time for the breast height samples, averaged over the six GSLs, are shown in Fig. 3. Variation in ring density was mostly explained by EW density variation (P < 0.0001, r = 0.964). LW density and maximum LW density also correlated very strongly (P < 0.0001, r = 0.970); thus, maximum LW density was used in the climate analysis for consistency with many dendroclimatology studies (Davi et al. 2003; Büntgen et al. 2010; Pompa-García and Venegas-González 2016). Drought years caused spikes in all the density components except LW density. The 1962 thinning to establish the experimental plots caused a sustained increase in LW density lasting for the duration of the study and reduced the influence of drought on LW proportion.

Fig. 3

PDSI and four measures of wood density at breast height, averaged over the six GSLs, vs. year (1940–2016). The wood density measures are ring density (RD; kg m⁻³), EW density (EWD; kg m⁻³), LW proportion (LWP), and LW density (LWD, kg m⁻³). Red solid lines indicate years of thinning (1962, 1972, 1982, 1992, 2003), and gray dashed lines indicate drought years commonly used as marker years in the Southwest (1951, 1956, 1963, 1977, 1989, 1996, 2002). Figure produced in R with some limited post-processing in Inkscape

3.1 Effects of growing stock level

After averaging values of our variables of interest over 1963 to 2016, GSL strongly influenced basal area increment (P < 0.0001) but had no significant effect on LW proportion (P = 0.07) or ring density (P = 0.12). Several pairwise comparisons between GSL levels on basal area increment were significant, and the effect of GSL was reduced at greater sample heights (Fig. 4). Regarding annual variability, year of ring formation was significant in most models and PDSI was significant in all models. Similar to the whole-sample analysis, GSL had no effect on wood density components, but did influence basal area increment (Table 2). The one exception was that GSL had a significant effect on LW proportion in the annual variability models (Table 2), but none of the pairwise comparisons were significant. Results are summarized in Fig. 5, which shows that GSL strongly influenced the long-term trend in basal area increment but had no effect on the long-term trend for ring density. A small effect on LW proportion is also evident; the lowest three GSLs appear to separate from the higher three GSLs, a trend that became increasingly apparent in recent decades.

Fig. 4

Sample-level mean values from 1963 to 2016 of basal area increment (BAI, top), ring density (RD, middle), and LW proportion (LWP, bottom) for the six GSLs at three sample heights. Letters indicate significant differences in BAI groupings at the indicated sample height (no differences in RD or LWP were significant). Figure produced in R with some limited post-processing in Inkscape

Table 2 Significance of fixed effects terms and model fit statistics. The sample height × GSL interaction, PDSI, and pre-1962 values were significant at < 0.0001 for all models, so for clarity, they are not included in the table

Fig. 5

Long-term trends in basal area increment (BAI, top), ring density (RD, middle), and LW proportion (LWP, bottom) for breast height samples. Lines were produced by fitting cubic smoothing splines to the data with degrees of freedom restricted to 5. Figure produced in R with some limited post-processing in Inkscape

3.2 Effects of climate

The analysis in treeclim revealed that precipitation had stronger effects on the response variables than temperature (Fig. 6). Total precipitation had a significant effect on ring density and EW density in the second (current year January–March) and third (current year April–June) quarters of the water year, as well as at the annual level. These effects were similar between the three GSL groupings. Second quarter precipitation correlated negatively with LW proportion (low coef = − 0.18, mid coef = − 0.17, and high coef = − 0.15), leading to a significant annual effect at the highest two GSL classes (mid coef = − 0.11 and high coef = − 0.14). First quarter precipitation had a positive correlation with maximum LW density in all three GSL groupings (low coef = 0.26, mid coef = 0.27, and high coef = 0.19), while fourth quarter (current year July–August) precipitation in the lowest GSL correlated negatively with maximum LW density (coef = − 0.21). There were fewer significant correlations with temperature, but the most notable was between third quarter average temperature and both ring density (mean group coef = 0.21) and EW density (mean group coef = 0.19).

Fig. 6

Response coefficients representing correlations between quarterly temperature (top row) and precipitation totals (bottom row) with breast-height wood density components at low (6.9 and 13.8 m² ha⁻¹), mid (18.4 and 23 m² ha⁻¹), and high (27.5 and 34.4 m² ha⁻¹) GSLs. Wood density components are ring density (RD), LW proportion (LWP), EW density (EWD), and maximum LW density (MXD). Significant correlations (confidence interval not overlapping zero) are indicated with asterisks. Figure produced in R with some limited post-processing in Inkscape

4 Discussion

4.1 Suitability of ponderosa pine for wood products

The density values for southwestern US ponderosa pine in this study are comparable to many other commercially important conifer species in the western USA. The overall mean density at 6% moisture content was 446 kg m⁻³, which is equivalent to 464 kg m⁻³ at 12% moisture content. This value is comparable to density of other interior west species (reported at 12% moisture content), such as lodgepole pine (Pinus contorta Douglas ex Loud.) (465 kg m⁻³), Engelmann spruce (Picea engelmannii Parry ex Engelm.) (368 kg m⁻³), and white fir (Abies concolor [Gord. & Glend.] Lindl. ex Hildebr.) (417 kg m⁻³; Alden 1997), but well below interior-west Douglas-fir (Pseudotsuga menziesii [Mirb.] Franco) (500 kg m⁻³). There may be other challenges with processing southwestern US ponderosa pine related to tree form and knot size, but wood density should not be a limiting factor to utilization. Although our results are obtained from a single site, they are similar to (but slightly above) previous published average values of wood density (12% moisture content) for ponderosa pine (449 kg m⁻³; Alden 1997).

4.2 Effects of stand basal area

In contrast to studies of wood density in gradual-transition species such as in the spruce genus (Dutilleul et al. 1998; Gardiner et al. 2011), increased stand density did not significantly affect any wood density components of ponderosa pine trees in this study, though it was associated with a significant but small increase in latewood proportion. These findings are similar to those of previous studies of abrupt-transition species such as Scots pine (Pinus sylvestris L.) (Peltola et al. 2007) and loblolly pine (Megraw 1985), which showed no effect of growth rate on wood density. From a management perspective, this indicates that dense slow-growing stands of ponderosa pine in the US Southwest do not necessarily produce wood with high density. Similarly, faster growing stands resulting from thinning operations will not necessarily have low wood density. Because trees in this study were 43 years old when initially thinned, our conclusions about the effects of stand density extend only to outerwood (mature wood). Thinning stands during the corewood (juvenile) phase may cause an expansion of the corewood and negatively impact wood quality (Larson et al. 2001).

4.3 Effects of climate

Climate, and particularly precipitation, had a strong effect on wood density and growth in this study. We found that drought years were associated with a sharp increase in ring density and earlywood density, but had a smaller effect on latewood proportion and maximum latewood density (Fig. 3). Maximum latewood density is often used in temperature reconstructions (Davi et al. 2003), but our findings suggest that earlywood density may be more useful in reconstructing past droughts in the Southwest, as measured by PDSI. Drought acts to increase earlywood density either through reduced cell expansion due to low turgor pressure, or through increased support of the hydraulic conduits to better withstand strong tension caused by water stress (Rodriguez-Zaccaro and Groover 2019). These two hypotheses are not mutually exclusive; the former suggests a mechanism underlying the response, while the latter suggests an adaptive response of trees to drought stress. It is worth noting the recent study by Candel-Pérez et al. (2018), in which they found that drought decreased ring density of Scots pine in northern Spain. This decrease was driven mostly by a decrease in latewood density, which we also observed during some drought years in the present study (Fig. 3). The increased earlywood density that we observed was absent in the Candel-Pérez study, which helps to explain the different response in overall ring density.

When analyzed by water year quarter, the most important season of precipitation was the second quarter (current year January–March). Other studies in southwestern US ponderosa pine have shown this to be one of the most significant time periods influencing stem radial growth (Adams and Kolb 2005; Kerhoulas et al. 2013). Spring precipitation (defined as March–May in this case) and vapor pressure deficit (VPD) are also influential at a stand level; they predict deviations from the historical range of variability in stand density of southwestern US ponderosa pine forests (Rodman et al. 2017). In the present study, increased precipitation during the second quarter was associated with decreases in both latewood proportion and earlywood density (Fig. 6), which had the combined effect of reducing the overall ring density. The decreased latewood proportion likely resulted from increased earlywood width due to more soil water availability earlier in the growing season. The decreased earlywood density was likely caused by sufficient turgor pressure to facilitate cell expansion. The other strong seasonal influence of precipitation was the correlation between first-quarter precipitation (October–December of previous year) and maximum latewood density. These results are similar to the those of Pompa-García and Venegas-González (2016), who found a positive association between maximum latewood density and December–February precipitation in Cooper pine.

We found that warmer pre-monsoon periods (April–June) were associated with increased ring density, while warmer monsoon periods (July–August) were associated with decreased ring density (although the latter was not statistically significant; Fig. 6). Warmer temperatures increase VPD and amplify water stress (Eamus et al. 2013). This is particularly important for southwestern US ponderosa pine in the spring and early summer, the most intense seasonal period of water stress in most years due to low spring precipitation and low humidity (Kolb et al. 1998; Gaylord et al. 2007). This could increase earlywood density by low turgor pressure and cell expansion. Effects of high temperature on VPD during the monsoon season, however, may be mitigated by the increased air humidity and greater soil water availability for evapotranspiration. Warmer temperatures likely have smaller effects on tree turgor, cell expansion, and earlywood density later in the summer than in the pre-monsoon season. However, warmer temperatures may increase respiratory demand (Ryan et al. 1995), reducing available carbon for cell wall thickening.

The southwestern USA is projected to become increasingly arid throughout the twenty-first century (Seager et al. 2007). Our results suggest that future droughts will decrease ponderosa pine wood volume production, but the wood produced during drought will be denser. From a wood products standpoint, individual trees will produce a lower volume of wood, but mechanical properties associated with density will likely improve. Regional wood supply, however, will likely be reduced by drought-driven disturbances, such as wildfire and bark beetle outbreaks that can rapidly kill many trees (Allen et al. 2010; McDowell et al. 2015; Kolb et al. 2016). Regarding estimation of forest carbon storage, our results underscore the need to account for changes in wood density due to drought. Moreover, results are valid only for the range of climate conditions explored in this study. Extreme and prolonged droughts in the future may not have the same effects on wood density as the relatively short-term and episodic droughts we investigated in this study between years 1919 and 2016.

4.4 Study limitations

One caution in extending these individual-tree results to the GSL level is that the DBH of sampled study trees was 3 to 6 cm smaller than the average in all GSLs except the 6.9 m² ha⁻¹ level. This is likely because the trees were marked for thinning before the sample trees were selected. Dominant trees are favored for retention in GSL studies, so the trees available for study were smaller in most of the treatments. However, the trees we selected survived five previous thinnings, so we believe that they are representative of the treatment levels. Furthermore, study trees are more representative than residual trees of the material that would in practice be removed in restoration treatments.

Although we had an adequate amount of data for the study, a larger sample size would have been necessary to answer questions regarding the interaction between GSL and climate, and about the short-term effects of stand thinning on the responses. Figure 6 shows some evidence of a GSL × climate interaction, but incorporating this interaction into models would have used up more degrees of freedom than we could accommodate with the data. Also, testing for an effect of thinning on individual samples would have required control trees against which to make comparisons. Because no trees in the control plot at Taylor Woods were cut in the 2017 thinning, such comparisons were not possible. This provides opportunities for future research, which would not necessarily require a replicated density study, but simply paired thinned and unthinned sites. Finally, we acknowledge that soil water storage, and its potential to buffer against drought, could have introduced some unexplained error to the climatic analysis. In the absence of any soil moisture holdover between years, the confidence intervals presented in Fig. 6 would have likely been narrower. To fully address the magnitude of these effects, we would need soil water content measurements for the study years, which were unavailable.

5 Conclusions

In this study, we investigated whether long-term management of stand basal area and short-term climate fluctuations affected growth and wood density components in a replicated ponderosa pine stand density experiment in Northern Arizona. We found that increased stand basal area did not affect any of the density responses but did have a strong negative effect on tree growth and a small but significant positive effect on latewood proportion of annual rings. Climate strongly affected the responses, mainly by increasing earlywood density in drought years. The average density of ponderosa pine wood at our southwestern US study site was similar to other commercial species in the western USA. Future droughts, at least within the ranges of duration and severity explored in this study, may lead to decreased tree volume growth but higher wood density in ponderosa pine.