Climate has a considerable influence on tree growth. Forest managers benefit from the empirical study of the historic relationship between climatic variables and tree growth to support forest management frameworks that are to be applied under scenarios of climate change. Through this research, we have utilized long-term permanent sample plot records, historic climate data sets, and linear mixed modelling techniques to evaluate the historic influence of climatic variables on the growth rates of major boreal tree species in Newfoundland and Labrador, Canada. For the commercially significant spruce and fir forests of the province, we found growing degree-days (GDD) to negatively correlate with tree productivity in warmer regions, such as much of Newfoundland (±1350 GDD), but positively correlate with growth in cooler regions, such as those in Labrador (±750 GDD). With respect to precipitation, environmental moisture was not on average a limiting factor to species productivity in the province. These dynamics have implications for the productivity of the spruce–fir forests of the study area when considered alongside contemporary climate projections for the region, which generally entail both a warmer and wetter growing environment.
Le climat exerce une influence considérable sur la croissance des arbres. Les gestionnaires forestiers bénéficient de l’étude empirique de la relation historique entre les variables climatiques et la croissance des arbres pour supporter le choix des pratiques d’aménagement forestier à appliquer dans le contexte du changement climatique. Dans le cadre de cette recherche, nous avons utilisé des données à long terme provenant de places-échantillons permanentes, des données sur le climat passé et des techniques de modélisation linéaire mixte pour évaluer l’influence des variables climatiques sur le taux de croissance des principales essences boréales arborescentes à Terre-Neuve-et-Labrador, au Canada. Dans le cas des forêts d’épicéa et de sapin de la province ayant une importance économique, nous avons trouvé que les degrés-jours de croissance (DJC) sont corrélés négativement (±1350 DJC) avec la productivité des arbres dans les régions plus chaudes, telles que dans la plupart des endroits à Terre-Neuve, mais sont corrélés positivement (±750 DJC) avec la croissance dans les régions plus froides telles qu’au Labrador. En ce qui a trait aux précipitations, l’humidité du milieu n’est pas en moyenne un facteur limitant pour la productivité des différentes essences dans la province. Ces dynamiques ont des implications pour la productivité des forêts d’épicéa et de sapin dans la zone à l’étude lorsqu’on tient compte des prédictions contemporaines du climat qui annoncent des conditions de croissance à la fois plus chaudes et plus humides pour la région. [Traduit par la Rédaction]
Accurate predictions of forest growth are a fundamental component of any forest management framework (Peng 2000). Empirical models of forest growth and yield (G&Y) have been used as the basis of forest management frameworks for over two centuries (Paulsen 1795; Subedi and Sharma 2011; Sullivan and Clutter 1972). Such models are based on historic records of growth and outputs are often given as timber yield, or basal area, as a function of diameter or volume (Porté and Bartelink 2002). Simple parameterizations and high predictive accuracy have led to the gradual establishment of the empirical approach to G&Y as the operational standard in modern forestry management applications (Vanclay 1994). Despite their current prominence, the use of many empirical models rests on an underlying assumption of a constant growing environment (Vanclay 1994; Yaussy 2000; Kimmins 2004). As a result, empirical models often have limited predictive capability when long-term changes occur in growing conditions.
Climate has a direct impact on tree growth and overall forest productivity (Morison and Morecroft 2006; Taylor et al. 2017), as well as indirect impacts on forest functions and processes, including disturbance (Boulanger et al. 2013, 2014) and competitive interactions (Oboite and Comeau 2019; Adler et al. 2012). As such, there is growing scientific consensus that the long-term changes associated with certain emission scenarios pose significant risk for abrupt and irreversible impacts on the composition, structure, and function of forest ecosystems around the world (IPCC 2014). Indeed, North American forests are already showing signs of complex structural change under the influence of climate change (Hogg et al. 2017; Tremblay and Boudreau 2011; Michaelian et al. 2011). The impact of projected climatic change on forest growth and composition is expected to vary regionally, contingent on local climate and existing forest conditions. Quantifying the climate impacts on forest growth and composition has therefore become a shared challenge for ecologists and environmentalists internationally in recent time (McMahon et al. 2011).
Forest managers often maintain databases containing extensive records of measurements taken at periodic intervals on static sites across the management area. These records, commonly referred to as permanent sample plot (PSP) inventories, generally include measurements of tree diameter at breast height (DBH), tree height, live crown, and (or) tree age. PSP inventories can provide valuable insight into tree growth over time and across diverse geographies and are used as key inputs in G&Y models (Sullivan and Clutter 1972; Subedi and Sharma 2011).
In Newfoundland and Labrador (NL), Canada, the province’s PSP data set contains records of tree growth across a diversity of conditions, ranging from highlands to coastal and subtundra forest ecosystems. The objective of this study was to evaluate the potential impacts of climate variability on the growth response of major boreal species in the province based on the historic empirical relationships between species growth and climatic terms. NL’s long-term PSP inventory can be used in combination with spatially explicit records of historic climate (McKenney et al. 2011, 2013) to empirically investigate the historic response of species to temperature and day and precipitation regimes. A linear mixed effect (LME) modelling approach was used to evaluate species-specific growth response to numerous climatic variables. While a number of fixed-effect environmental parameters were included within the models in an effort to minimize random error, this study is centred around those growing environment variables that are expected to shift under scenarios of climate change, namely annual cumulative growing degree-days (GDDs) and cumulative annual precipitation (Finnis and Daraio 2018).
Materials and methods
Our study area was the province of Newfoundland and Labrador (NL) (Fig. 1). Newfoundland (NF) has a land area of 11.1 million ha, of which nearly half is forest area. Labrador (LD) is home to 18 million ha forest, which makes up 60% of its overall land area (NLDNR 2014). There are 21 tree species native to the province. Balsam fir (Abies balsamea (L.) Mill.) is the most abundant species in NF and can be found frequently in pure stands in the western areas of the island. Repeated fires and highly acidic, rocky, moist soils have also established black spruce (Picea mariana (Mill.) B.S.P.) as a dominant species in central NF. On nutrient-rich, moist upland soils and where prevailing climate is moderate, white spruce (Picea glauca (Moench) Voss) can be prevalent (NLDNR 2014).
In LD, the combination of a cooler, dryer climate, and the fine-textured soils of lowland areas allows black spruce to prosper in association with balsam fir, white birch (Betula papyrifera Marsh), balsam poplar (Populus balsamifera L.), trembling aspen (Populus tremuloides Michx.), and white spruce. Elsewhere in LD, growing conditions are limited by poor drainage, harsh climate, and thin, nutrient-poor soils. These landscapes are characterized by subarctic forests of pure black spruce and mixtures of black spruce and balsam fir (NLDNR 2014). While deciduous species have not formed major pure stands in the province, they are often significant components of mixed-wood stands. The full list of species within the scope of this study are detailed in Table 1.
Regional climate is summarized in Table 2. Precipitation rates vary across NF (750 to 2000 mm·year−1; Table 2). Annual precipitation tends to be less in plots located in LD compared with those in NF (725 to 1250 mm·year−1; Table 2). In both regions, coastal settings benefit from the moderating effects of the ocean on ambient temperatures, but generally experience greater precipitation volumes per annum than inland environments do, particularly through winter months. January mean daily temperatures range from −4 to −10 °C across NF and from −16 to −25 °C across LD. July mean daily temperature is ±16 °C across NF and ranges from 6 to 13 °C in LD. Average daily extreme temperature in both seasons vary about 5 °C from the daily mean, with inland regions often experiencing greater variability. Summary climate information was obtained from Environment and Climate Change Canada (Government of Canada 2021).
Permanent sample plot data
The NL Department of Fisheries and Land Resources provided their PSP data, which has coverage of much of the province (Fig. 1). The database spans a period from 1985 to 2017 and includes a network of 1130 plots and 640 000 tree measurements. As with many PSP inventory programs, plot re-measurement does not occur on an annual interval in NL, but rather at periods typically ranging from 3 to 7 years. Because of the sporadic interval of plot re-measurement and because growth can only be measured from the second plot measurement onward, 360 000 growth measurements were derived from the PSP inventory.
The database structure separates mature and immature plots. While tree age is widely recognized to be an important influence to tree growth (Johnson and Abrams 2009; Sillett et al. 2010; Stephenson et al. 2014), in this research, mature and immature plot data sets were combined. Tree observations associated with dead trees, or where no DBH measurement was recorded (±150 000 observations combined), were also removed from the pool of viable tree measurements.
This study assumed net tree growth to be continuous throughout tree life (Stephenson et al. 2014). Tree observations found to have negative or no average annual growth between re-measurements were assumed to be a consequence either of observational or recording error and were removed from the pool of viable growth measurements accordingly. Tree observations immediately preceding, and subsequent to, such observations were also removed. The rationale for the latter filter being that one error in a sequence inhibits a viable growth value being determined between preceding and subsequent observations where growth is calculated as the difference between observations. Approximately 62 500 tree observations were removed through this filter.
Plots where disturbance regimes (i.e., moose browsing, fire, pests, and windthrow) had documented influence on resident tree growth were also removed from the pool of viable plots in this analysis. Through this filter, 259 plots and ±600 plot establishment and re-measurement observations were removed from the study. The pool of viable tree observations was ultimately reduced to 130 000 growth observations following all data set manipulations. The distribution of the NL PSP network is portrayed in Fig. 1 and summarized in Table 1. A box–whisker plot detailing median tree growth observation by species and region, as well as interquartile range and minimum and maximum growth values, is provided as Fig. 2.
Plots range in size from 20 m2 to 1000 m2 in NL’s inventory and include an array of intermediate sizes. Further, plot size is not fixed in NL’s inventory and can increase as a plot matures and stem density naturally declines. Stand-level predictors in this study were therefore calculated as per-hectare averages to accommodate variable sizes across plot measurements. Spatial variables (i.e., slope and potential solar radiation) were calculated as the mean value within a 30 m radius from the plot centre.
Our LME models included three fixed-effect climatic predictors, two fixed-effect tree-level predictors, up to four fixed-effect stand-level predictors, and random effects for plot with trees nested within plot. Interaction terms were included in two of the models to identify species-specific and regional variations in growth response to the precipitation, GDD, and (or) slope predictors. A complete list of the variables, including information on variable type and class, is provided through Supplementary Table S11.
Historic monthly cumulative precipitation, monthly mean daily minimum and mean daily maximum temperatures were estimated for every PSP from 1985 to 2010 (McKenney et al. 2011, 2013). The data sets provided were generated using thin-plate smoothing splines, as implemented by ANUSPLIN (Hutchinson 1995). ANUSPLIN is a widely used “nonparametric” surface fitting method that can generate estimates of certain climatic variables at any location where the independent variables are known, in this case longitude, latitude, and elevation. The historical monthly models have 95% confidence limits of approximately ±1.1 °C for maximum temperature, ±1.3 °C for minimum temperature, and 10%–40% for monthly precipitation (McKenney et al. 2011). Larger margins of error are associated with minimum temperature estimates versus maximum, typically because station data coverage does not accurately adjust for effects such as cold air drainage. Limited station coverage also affects the ability of the models to capture spatial and temporal variation in precipitation estimates.
Cumulative annual precipitation, cumulative annual GDDs, and cumulative annual potential incoming solar radiation were included as climatic variables in the LME models. Various methods for the estimation of annual cumulative GDD from monthly mean minimum and mean maximum temperatures were evaluated (Erbs et al. 1983; Schoenau and Kehrig 1990; Thom 1954). This study used the GDD estimation approach employed in the JABOWA (Botkin 1993) family of forest gap models. The JABOWA estimation method was selected because JABOWA estimates cumulative annual GDD from monthly averages of mean daily mean minimum and mean maximum temperatures only, through an assumed sinusoidal annual temperature profile. This study had only daily mean minimum and mean maximum monthly temperatures on which to base its estimation of GDD on:
where DEGD is the estimated cumulative annual GDD, Tjuly is the mean monthly temperature in July, and Tjan is the mean monthly temperature in January, and 4.4(°C) is the base temperature below which growth is assumed to be negligible. In our study, monthly mean daily temperature was estimated as the mean of the monthly average daily minimum and maximum temperatures (see Botkin 1993).
We used a geographic information system to obtain elevation data for each plot in the PSP inventory (ESRI 2018). Elevation data was collected in GeoTiff format from the Shuttle Radar Topography Mission (SRTM) one arc-second digital elevation model (DEM) (Farr et al. 2007). The SRTM global DEM has a ±30 m resolution with near global coverage. Potential solar radiation was calculated for each PSP based on methods from the hemispherical viewshed algorithm developed by Rich et al. (1994) and further improved by Fu and Rich (2000). Through their approach, the total amount of radiation calculated for a particular location or area is given as global radiation. The calculations of direct, diffuse, and global insolation are repeated for each feature location or every location on the topographic surface, thereby producing insolation maps for the entire area.
Tree-level fixed-effect variables
Data for tree-level predictors are field measurements recorded by the NL government as part of their forest inventory program. DBH was the only continuous variable at the tree level, and tree species was the only nominal variable. The tree species term (SPC) was coded into nine levels, and balsam fir was established as the reference level (Table 1). Certain species were not explicitly documented in NL’s PSP database and were instead coded in the inventory under one of three general tree classes. These classes were removed from this analysis, as it had to be assumed these groups contained trees of more than one species, and therefore, growth response may differ to climatic terms within the class; ±800 tree observations were removed by this means.
Stand-level fixed-effect variables
Stand-level variables were used to assess the influence of stand-level competitive dynamics within models as well as the significance of an individual trees relative competitive position within a stand. The tree stand index (TSI) is a percentage ratio: the ordered position of a tree within all trees in a plot by DBH, to the total number of trees in the plot. Index values range from zero to 100; the largest diameter tree(s) within a plot scoring an index of 100 and the smallest diameter tree(s) an index of zero. While the TSI term captures only the relative DBH of a tree as compared with other trees in a plot, tree DBH has been assumed in the northeastern boreal study area to generally correspond to tree height. This is a particularly reasonable generalization because species diversity is low within NL. Trees with greater DBH and, therefore, presumed greater individual height would have greater access to sunlight. Species-specific variations in typical height–diameter relationships were inherently not accounted for through this approach. However, where tree height was only recorded on a sample tree basis through NL’s PSP program, DBH measurements offer the most feasible means of accommodating a trees competitive position within models.
Stand-level variables based on basal area calculations were included in models to reflect the influences of stand crowding on tree growth trends (BAS). As plot size varied within NL, basal area calculations were adjusted to cm2·ha−1 measurements. We calculated percent deciduous basal area (CHW) as an additional metric to accommodate the influence of stand composition on tree growth. CHW being the proportion of deciduous individuals within the plot. As with TSI, the CHW variable scaled from 0 (no deciduous competition in plot) to 100 (no conifer competition in plot). The CHW variable, by design, did not attempt to accommodate important and complex species-specific interactions; LME methods would not generally be the most suited approach for sophisticated analysis of interspecies competitive interactions.
A slope geoprocessing operation was run against the DEM to derive a continuous percentage slope surface with full coverage of the study area. From this layer the median percentage slope value was extracted for each plot in the inventory within a diameter corresponding to the width of the plot. As with the potential solar radiation predictor, the percentage slope predictor was based on the SRTM DEM. The study area includes both plots within the island of NF and inland LD. We included the plot region (REG; i.e., NF or LD) as a predictor to reflect differences in long-term climate normals within the larger study area of NL.
Growth increment (ΔBA) in this study refers to the mean annual change in basal area between two measurements of a tree taken at different points in time. Where plot re-measurement (and by extension tree re-measurement) occurred at inconsistent intervals through the PSP program, growth was averaged per annum across the number of years lapsed since the previous measurement. ΔBA observations greater than 25.0 cm2·year−1 were considered to be observation and (or) recording error and were removed from the analysis accordingly (±50 observations).
A natural log transformation was applied to the dependent growth variable to meet the assumption of homogeneity in residual variance.
where ΔBA is the average annual basal area growth increment, with a mean of μ and a variance of σ2, and where a ln(ΔBA) dependent variable approximately follows a normal residual distribution.
In many cases, the coefficient of variation would be quite small. Where the coefficient of variation is small:
With a one-unit increase in a predictor term, and provided all other predictor terms remain unchanged, the response in the log scale increases, on average, by
where βi would approximate zero. Given this, variable estimates for a function fitted to ln(ΔBA) can be interpreted as a percentage change in the ΔBA dependent per unit change in a predictor. Log-normal annual ΔBA was initially evaluated as a function of stand level, tree level, and climatic variables (i.e., the “baseline” model):
where ΔBA is the average annual basal area ΔBA, β0 is the intercept, β1 through β8 are the coefficients for the corresponding predictor terms, SPC modifies the intercept by an estimate associated with tree species, and j indexes a plot-level variable. The random effects of plot (PLT, n = 1130) and trees (TRE, n = 117 395) nested within plot were used to account for spatial and temporal dependencies. LME models were used in determining coefficients for all initial models (Bates et al. 2015). To test for regional differences in climate response regional interaction (REG; Supplementary Table S11) was then introduced into eq. 6, yielding eq. 7:
where subscript k indexes a variable by region (i.e., NL or LD, the “regional” model). Regional interaction was only applied to climatic terms, as the scope of this study was centred around growth response to PCP and GDD variables. A regional interaction term was applied to the SLP variable, as it was anticipated that the term may adjust for environmental moisture conditions rather than for relief exclusively. Other predictors (i.e., DBH, BAS, TSI, CHW, and INS) may too have differing growth response between regions, but these dynamics were not considered within the study scope. To further evaluate differences in growth response between species and regional settings, species interaction, s, was introduced to the GDD and PCP predictors, rather than as a nominal variable in its own right (i.e., the “species-regional” model), yielding eq. 8:
Spearman rank correlation (Supplementary Fig. S11) and variance inflation factors (Supplementary Table S21) were used to test the independence of model covariates prior to model analysis. Generalized variance inflation factors for each independent were <3 (Supplementary Table S21); values not generally indicative of multicollinearity between covariates (Menard 1995). Spearman rank correlation values between covariates used in this analysis would not generally be considered indicative of significant correlation (Supplementary Fig. S11). LME models were prepared using the lme4 and lmerTest packages (Bates et al. 2015; Kuznetsova et al. 2017) in R (R Core Team 2019).
Marginal and conditional coefficient of determination were used to evaluate the proportion of variance explained by the models, as well as the proportion of variance explained through the plot and tree random effect terms (Moore et al. 2013). Restriction maximum-likelihood (REML) criterion at convergence, residual mean-squared error (RMSE), Bayesian information criterion (BIC), and Akaike information criterion (AIC) were used to evaluate and compare model fit. Coefficient estimates from LME models were used to evaluate the nature and magnitude of species-specific relationships to climatic variables, as well as the influence of the regional term. Student’s t test was used to evaluate the strength of the relationship identified.
Marginal R2 exceeded 0.45 in each of the models prepared, whereas conditional R2 exceeded 0.85 (Supplementary Table S31). The explained variance for the plot and nest tree random effect terms in all models was relatively large as compared with residual variance (Table 3, Supplementary Tables S4–S51). Considering (i) the separation between the conditional and marginal R2, and (ii) the variance explained by the random effect terms as compared with the residual variance, we can conclude that the random effect terms explain a considerable proportion of the variance in our models and were justified in their inclusion (Supplementary Table S31).
RMSE was 3.1 cm2 in all models, which is relatively low given the range of ΔBA values in the PSP data set (Fig. 2). Use of a log-normal growth dependent was found to improve the distribution of residuals considerably as compared with models where no transformation was applied to the dependent (Supplementary Fig. S21). Of the models prepared, the species-regional model was found to have optimal BIC and AIC values and would be considered the best-fitted model accordingly; although the difference in fit between all models as determined through the aforementioned indicators was quite small.
Stand-level variables were found to be highly significant to species-specific growth response in all regressions. The BAS predictor was found to have a negative relationship to individual tree growth; −1.9% per cm2·ha−1 and a t value of <−50.0 in all models. TSI was found to positively influence growth at ±1.5% per index rank and was significant (t value > 90.0) across all models. CHW was also found to be significant (t value < −10.0) across all models and negatively influenced tree growth at a rate of approximately −1.3% per 1% deciduous composition in a stand (Supplementary Table S41).
In terms of base productivity, black spruce was the only species found to have an average rate of growth faster than the balsam fir reference level in NL (+8.9% ΔBA·year−1). Most deciduous species, however, including white birch (−22.0% ΔBA·year−1), red maple (−28.7% absolute ΔBA·year−1), and trembling aspen (−11.9% ΔBA·year−1), had an average growth rate less than the balsam fir reference level (Supplementary Table S41).
Through the regional model, GDD was found to have a positive and significant mean relationship to growth in LD, at a rate of 0.02% annual increase in productivity per additional GDD within a year, but negative and significant relationship in NF, at a rate of −0.02% decrease in annual productivity per additional GDD (Supplementary Table S51). Conversely, PCP was found to have a significant and negative relationship to productivity in LD (−0.06% annual decrease in productivity·mm−1·year−1) and a significant and positive relationship to productivity in NF (0.004% annual decrease in productivity·mm−1·year−1; Supplementary Table S51). Where no regional interaction was included in the baseline model (Supplementary Table S41), the contrasting nature of regional growth response to climatic terms either cancelled out entirely rendering the variable insignificant (i.e., PCP) or otherwise favoured the growth response in NF (i.e., GDD). NF is home to a much greater proportion of plots constituting the provincial PSP inventory, and thereby has greater influence on overall model fit (Table 1).
The regional GDD response of black spruce, balsam fir, and white spruce was found to be positive in LD (0.01%, 0.02%, and 0.06% annual increase in productivity·GDD−1·year−1, respectively) and negative in NF (−0.03%, −0.005%, and −0.05% annual decrease in productivity·GDD−1·year−1, respectively) through the species-regional model (Table 3; Fig. 3). White birch was the only species found to have a positive growth response in both regions: 0.04% and 0.01% annual increase in productivity·GDD−1·year−1 in LD and NF, respectively (Table 3; Fig. 3). Species with less growth observations, including red maple, trembling aspen, white pine, and yellow birch (Table 1), yielded estimates for growth response in both NL and LD that were less significant or otherwise altogether insignificant as compared with species with greater representation (Table 3; Fig. 3). Relatedly, the relationship between white pine and trembling aspen to the PCP variable was insignificant, both in NF and LD.
Of those species included in this research scope (Table 1), three species were found to have a positive growth relationship to the PCP variable in NF: black spruce (0.03% annual increase in productivity·mm−1·year−1) and white spruce (0.03% annual increase in productivity·mm−1·year−1; Table 3; Fig. 3). No species was found to have a positive relationship to PCP in LD. Of the boreal species included in this study, white spruce was found to have a negative growth relationship of the largest magnitude to PCP in LD (−0.14% annual decrease in productivity·mm−1·year−1). Black spruce was also found to have a negative growth relationship to PCP in LD (−0.07% annual decrease in productivity·mm−1·year−1), but a positive growth relationship to PCP in NF. Of those species included in this study, yellow birch had a negative growth relationship to GDD of the greatest magnitude in NF (−0.1% decrease in productivity·GDD−1). Balsam fir, white birch, and red maple were found to have a negative growing response to PCP in both the NF and LD regions (Table 3; Fig. 3).
INS was found to be insignificant across all models. The SLP variable was observed to relate positively with growth in NF (Supplementary Table S51), but had no significant relationship in LD. With no regional interaction in the baseline model (Supplementary Table S41), the SLP variable was altogether insignificant.
The negative relationship between tree productivity and precipitation for the majority of the boreal species included in this study suggests that either environmental moisture is not generally a limiting factor to tree growth in the province or that the precipitation covariate is capturing the effects of highly collinear variables, such as cloud cover. In either case, excess precipitation or variables highly collinear to precipitation, such as cloud cover, seem to have a negative influence on tree productivity for most species in the study area. This characterization was particularly true in LD, where growth response to precipitation was found to be negative for all species.
LD’s generally cooler growing environment may have historically acted as a constraining factor to evapotranspiration potential as compared with similar latitudes in NF, where ambient temperatures are relatively warmer. Reduced evapotranspiration potential in LD as compared with NF could lend to a surplus of moisture in the growing environment in LD, but a deficit in NF. Further, the GDD covariate used in this study would be collinear and capture the influence of other important temperature-related variables including seasonal high temperatures (Ettl and Peterson 1995; Lloyd and Fastie 2002). Such dynamics help to rationalize the positive response to GDD in NF, but negative response in LD. There is some expectation that potential evapotranspiration regimes will be subject to change under the influence of 21st century climate projections (Chattopadhyay and Hulme 1997; Kingston et al. 2009). This would be particularly true in NL, where regional projections forecast both a warmer and wetter growing environment (Finnis and Daraio 2018), and where the GDD growth response of conifer species in LD was generally identified to be positive through this work. Such conditions could very well lend to increased evapotranspiration rates in LD and would likely diverge from the historical growth relationships to precipitation identified through this work. Furthermore, many of NL’s native species including black spruce, balsam fir, tamarack, and trembling aspen are tolerant of a diversity of soil moisture conditions and would generally be considered to be well adapted to those present in LD (Burns and Honkala 1990a).
In NF, commercially significant softwoods of the province (i.e., black spruce, balsam fir, and white spruce) were observed to have a negative growth relationship to GDD, and in the case of balsam fir, a negative growth relationship to precipitation as well. Where regional climate projections generally entail warmer and wetter conditions across the province through to 2100 (Finnis and Daraio 2018), our findings may support spruce, which were observed to have a positive growing relationship to precipitation, as having a slight ecophysiologic edge over fir under these projections. This would be particularly true for black spruce, which was found to have the fastest base rate of growth in the NF region (15.9% faster than that of balsam fir; Supplementary Table S51).
Still, the interrelated nature of climatic variables to tree growth through processes such as evapotranspiration (Kişi 2006) means changes in climatic variables are rarely as simple as a cause and effect outcome (Morison and Morecroft 2006; Taylor et al. 2017). Long-term changes to numerous climatic variables, including the cumulative annual precipitation, ambient temperature, permafrost, and the frequency and intensity of storm events (IPCC 2014), will undoubtedly influence the climate-tree growth relationships identified through this research. So, while increased GDD and precipitation may favourably position spruce under current climatic conditions and through current forest dynamics, other simultaneous changes in growing variables anticipated under scenarios of climate change make inference into net change in forest growth outcomes less certain.
The observed regional difference in GDD response remains significant, however, both in that it confirms GDD as a constraining factor to species-specific growth trends within the province and that these climatic terms might have markedly different influence over relatively small spatial scales within a wider global context (<1000 km in this study). While some studies have found the relationship between growth and GDD to be quadratically related (Lapointe-Garant et al. 2010), growth response varied in this study from one species to the next, and the specific nonlinear relationship might be best determined at the species level. The GDD growing conditions captured through NL’s PSP database is not reflective of the complete ecologic niche of any of the boreal species studied and therefore could not inform a nonlinear fit that would be valuable beyond applications within the study area and indeed beyond those conditions that currently exist within the study area. This holds true for those contrasting precipitation–growth responses identified for black and white spruce which were moderated by regional setting.
With respect to species’ growing niche, some studies have found that an expanded fundamental niche of temperate tree species resultant of contemporary climate forecasts may lend to the retreat of dominant conifer forests in forest transition zones in lower-latitude climates, such as the Atlantic Maritime provinces (Boulanger et al. 2017; Steenberg et al. 2011; Taylor et al. 2017). Concerns around boreal retreat stem from dramatic changes in regional competition dynamics between warm-adapted species like red maple (Abrams 1998) and cold-adapted species like balsam fir and black spruce (Ashraf et al. 2015). Generally, warm-adapted species such as red maple have greater productivity in warmer climates; NL being toward the northern extent of their North American realized ecologic niche (Burns and Honkala 1990b). Given the negative growth response to the precipitation variable observed through this study, and given the insignificant growth response to the GDD variable, we found no evidence of climate–growth relationships currently that would support red maple or other fast-growing temperate species having a ecophysiological advantage under scenarios of climate change versus the established balsam fir and black spruce forests of the province.
Warm-adapted species like red maple currently constitute a marginal portion of NL’s northeastern boreal forest (Table 1). Given they represent a small component of NL’s current forest composition, these species would have limited capacity to rapidly outpace and replace overwhelmingly conifer dominant forests under the typical horizon of most forest management planning frameworks, even under conditions which may better align with a more temperate species niche. This is not to suggest present day competitive stand dynamics will remain status quo within NL; shifts in forest competitive dynamics resulting from changes in climatic variables change may, and likely will, have implications on the growth relationships identified herein, as discussed.
Indeed, competitive interactions within boreal forest regimes are a significant factor in net forest productivity (Vose et al. 2012; Adler et al. 2012; Boulanger et al. 2017; Taylor et al. 2017). We found the metrics of competition significantly influenced growth rate. With reference to the stand density predictor, greater densities would be indicative of more stand crowding and greater competition for environmental resources. Accordingly, we found the growth relationship to this variable to be negative across all models. Conversely, the identified positive relationship between growth and the tree stand index predictor can be rationalized given larger trees within the stand generally have a greater competitive edge, greater exposure to environmental resources, and by extension greater opportunity for growth. Ecologic principals such as the law of self-thinning are unlikely to change with shifts in climate normals (Westoby 1984).
Interestingly, we found the deciduous competition variable to have a negative relationship to the growth dependent. Typically, the live crown in a deciduous forest is greater than that of a conifer forest at comparable stem densities. Live crown directly influences competitive interactions within a stand by impacting the availability of solar resources in the understory (Zeide 1987). The proportion of deciduous trees within a stand likely correlates with stand-level competitive intensity and, therefore, negatively correlates with the growth rate of individual trees within the stand.
DBH was found to have a significant and positive relationship to growth, exceeding a 7% ΔBA increase per centimetre DBH in all models. Tree growth has been found to accelerate or otherwise remain constant with tree age (Johnson and Abrams 2009; Sillett et al. 2010; Stephenson et al. 2014). We observed that larger trees experience greater increment growth. As this study utilized an explicit ΔBA term, it is somewhat intuitive that larger trees lend to larger area ΔBA. Further study using a proportional growth term (i.e., a percentage area ΔBA) may provide further insight into the significance of the DBH–ΔBA relationship.
Limitations and caveats
We acknowledge that there are limitations to the work conducted here. Indeed, the relationships between tree growth and climatic terms identified through this work were limited by the data available and could not address a full spectrum of potential sources of tree growth variation. Certain environmental variables such as disturbance (e.g., moose browsing, windthrow, silviculture, disease, pests, fire, etc.), soil characteristics, permafrost depth, extreme weather events, and snowpack can influence species’ growth trends. Any growth variation originating from variables not captured in the models herein is treated as random error in this study.
Age and the successional stage of trees was also not considered. Non-uniform growth response to climate change might be expected from longer-lived species, such as white spruce, where climate-driven changes in forest physiology (e.g., stand age and establishment) have been found to have greater impact their ability to adapt under a changing climate (Hogg et al. 2017). While plots in NL inventory are assigned a generalized stand age, use of this information in a study that sought to conduct analysis on the growth response of individual trees would not be appropriate.
Almost all PSPs used in this analysis are observed within natural forests and commonly contain notable intrastand variability in tree age. Further, less than 6% of the trees in the PSP inventory have had their age recorded (i.e., sample trees). The remaining tree observations within the inventory are assigned an age class by the field crew based on those ages sampled within the plot. Including individual tree age in the model not only introduces biases associated with tree age but also has the effect of limiting the model application. Simply put, we did not have access to robust age information to consider tree age or successional status as a variable. Similarly, tree height was only recorded for sample trees in NL’s inventory program. This inhibited the use of height information for any of the variables prepared in this work. It was for this reason the TSI term was based on measurements of DBH rather than height.
Though the LME models are not intended to be used to inform estimations of forest G&Y, the generalized approach could be used as one input into a forest G&Y modelling approach for the province of NL. Our models determine annual basal area ΔBA (cm2·year−1) through an empirical approach that includes explicit climatic variables as part of its formulation. Such an approach disassociates the relationship between historic growth patterns and future growth trends, where a coupled relationship is characteristic of most current empirical G&Y models, as discussed. This disassociation has the potential to allow key climatic variables to vary within the model horizon, as would be the case in any application in which a model is applied under scenarios of climate change. Further, the models prepared in this work are based on variables sourced from or derivative of data sets presently available in the province. Still, the LME approach prepared in this study only represents a component of a fulsome G&Y modelling method. Any G&Y framework to be applied in an operational context would generally require other key inputs, such as calculations for mortality and ingrowth.
We observed the growth response of species to cumulative annual GDD to contrast between sites in LD and NF for the major commercial species in the province. Further, we found that a small proportion of species responded positively to the cumulative annual precipitation variable, indicating that precipitation is not generally a positive contributor to species productivity within the northeastern Atlantic study area. Where regional climate projections entail both a warmer and wetter future climate across the study area, the commercially significant boreal conifer forests of the province may be anticipated to experience a reduction in forest productivity in warmer settings such as NF but a boost in productivity in cooler settings such as LD.
While the species included in our models supported a negative growth relationship to precipitation across the province and particularly in LD, this is likely more a consequence of limits on evapotranspiration potential associated with historic GDD normals. With warmer climate normals and increased evapotranspiration potential, the productivity relationships fitted through the models prepared in this work may be subject to change, and indeed likely will. Beyond evapotranspiration, the relationship between environmental growing variables and growth are complex and often highly interrelated. With climate change expected to bring forth long-term changes in numerous climatic growing variables beyond simply cumulative annual precipitation and GDD, the inferences put forth herein face some degree of uncertainty in practice, as they are rooted on the present-day interactions between growing variables explored through the empirical methods. Still, the empirical study of existing climate–growth relationships established through this work offers key insights that can better inform forest management practice in northeastern boreal forests. While we did not derive explicit predictions of forest G&Y though this work, the relationships identified herein can contribute to the basis of such efforts.
This research was supported in part by funding provided by the Newfoundland and Labrador Department of Fisheries and Land Resources (grant No. 221325).
Data availability statement
Restrictions apply to the availability of these data sets. Climate data was obtained from the Canadian Forest Service (Natural Resources Canada), and permanent sample plot data from the Province of Newfoundland and Labrador. These data sets are available with the permission of the respective organization.
Declaration on conflicts of interest
The authors declare that they have no conflict of interest.
We thank the Province of Newfoundland and Labrador for the provision of their PSP database and the Canadian Forest Service for the provision of their climate data sets. We would also like to express personal gratitude to John Pedlar as well as Kevin Searls and Tessa Williams for their varied contributions to the review and development of this manuscript.
Adler P.B., Dalgleish H.J., and Ellner S.P. 2012. Forecasting plant community impacts of climate variability and change: When do competitive interactions matter? Climate and species interactions. J. Ecol. 100(2): 478–487.
Ashraf M.I., Bourque C.P.-A., MacLean D.A., Erdle T., and Meng F.-R. 2015. Estimation of potential impacts of climate change on growth and yield of temperate tree species. Mitig. Adapt. Strateg. Glob. Change. 20(1): 159–178.
Boulanger Y., Taylor A.R., Price D.T., Cyr D., McGarrigle E., Rammer W., et al. 2017. Climate change impacts on forest landscapes along the Canadian southern boreal forest transition zone. Landscape Ecol. 32(7): 1415–1431.
Hogg E.H., Michaelian M., Hook T.I., and Undershultz M.E. 2017. Recent climatic drying leads to age-independent growth reductions of white spruce stands in western Canada. Glob. Change Biol. 23(12): 5297–5308.
Lapointe-Garant M.-P., Huang J.-G., Gea-Izquierdo G., Raulier F., Bernier P., and Berninger F. 2010. Use of tree rings to study the effect of climate change on trembling aspen in Québec. Glob. Change Biol. 16(7): 2039–2051.
Michaelian M., Hogg E.H., Hall R.J., and Arsenault E. 2011. Massive mortality of aspen following severe drought along the southern edge of the Canadian boreal forest. Glob. Change Biol. 17(6): 2084–2094.
NLDNR. 2014. Growing our Renewable and Sustainable Forest Economy: Provincial Sustainable Forest Management Strategy (2014–2024). Government of Newfoundland and Labrador, Department of Natural Resources, Newfoundland and Labrador, Canada. Available from https://www.gov.nl.ca/ffa/files/publications-pdf-psfms-14-24.pdf [accessed 9 March 2018].
Taylor A.R., Boulanger Y., Price D.T., Cyr D., McGarrigle E., Rammer W., and Kershaw J.A. 2017. Rapid 21st century climate change projected to shift composition and growth of Canada’s. Acadian Forest Region. For. Ecol. Manag. 405: 284–294.
Vose, J.M., Peterson, D.L., and Patel-Weynand, T. 2012. Effects of climatic variability and change on forest ecosystems: a comprehensive science synthesis for the U.S. Gen. Tech. Rep. PNW-GTR-870. US Department of Agriculture, Forest Service, Pacific Northwest Research Station, Portland, Oreg., USA.
Yaussy D.A. 2000. Comparison of an empirical forest growth and yield simulator and a forest gap simulator using actual 30-year growth from two even-aged forests in Kentucky. For. Ecol. Manag. 126(3): 385–398.
T. Searls, X. Zhu, D.W. McKenney, R. Mazumder, J. Steenberg, G. Yan, and F.-R. Meng. 2021. Assessing the influence of climate on the growth rate of boreal tree species in northeastern Canada through long-term permanent sample plot data sets. Canadian Journal of Forest Research.
51(7): 1039-1049. https://doi.org/10.1139/cjfr-2020-0257
If you have the appropriate software installed, you can download article citation data to the citation manager of your choice. Simply select your manager software from the list below and click Download.
1. 495-Year Wood Anatomical Record of Siberian Stone Pine (Pinus sibirica Du Tour) as Climatic Proxy on the Timberline