Implications of Liebig's law of the minimum for tree-ring reconstructions of climate

The principle that tree growth is determined by the most limiting factor, known as Liebig's law of the minimum (Sprengel 1828, Liebig 1840), underlies efforts at dendroclimatological reconstruction (Fritts 1976, Fritts and Swetnam 1989, Speer 2010). This principle guides sampling efforts towards trees that are particularly stressed by a common growth factor (CGF) experienced by all trees at a site, such as growing-season temperature or moisture availability (LaMarche 1982, Pilcher et al 1990). As one approaches elevational treeline, for example, trees become smaller, a clear visual indication of an upper-limit on growth rates that are often set by cold ambient temperatures (Körner 2012). It follows that limitation on growth at treeline will generally be relaxed in warmer years and permit for greater growth, provided that other, stronger growth limitations are absent.

In environments where growth is limited by temperature, the characteristics of annual growth rings can serve as a proxy of temperature variation. Dendroclimatic reconstructions of climate rely on sampling this growth variability, typically by removing a thin radial core of the wood, and measuring interannual variability in either width (Douglass 1920), density (Parker and Henoch 1971) or color (McCarroll et al 2002). After accounting for the effects of age-growth relationships, disturbance persistence, and the presence of outliers (Cook 1985, Cook 1992), order ten to a hundred records are typically averaged together to form a single chronology. Such a 'mean chronology' can be optimal if tree growth is controlled by a single CGF with various disturbances superimposed on individual trees (Cook 1985, Cook 1992), and this approach has proved successful in reconstructing annually-resolved variability in temperature and precipitation (e.g. Masson-Delmotte et al 2013).

But there exists a tension between averaging across individual records and a local interpretation of Liebig's law of the minimum. If Liebig's law applies at the level of the individual tree, not all records will always respond to the CGF. Tree growth could instead reflect either a CGF or local growth factors (LGFs), depending on which imposes the strongest limitation on a given tree in a given year (Chapin et al 2011). It is common for tree-ring models to invoke Liebig's law of the minimum (Fritts et al 1991, Shashkin and Vaganov 1993, Vaganov et al 2006, Tolwinski-Ward et al 2011). Such models are typically formulated in terms of rates of growth, and are thus expressed in terms of Blackman's law of limiting factors (Blackman 1905). The generalization proposed here posits that unobserved limitations that differ between trees across a site causes within-stand variability in limiting factors. Plausibility of such variability is supported by previous demonstration of environmental sensitivity varying within stands according to details of topographic setting (Bunn et al 2011, Salzer et al 2014, Tran et al 2017). In the following, we present a test for distinguishing between models of tree-ring growth for which all trees at a site experience identical forcing, from those for which trees experience multiple limitations upon growth. We then explore the implications of the latter scenario for where climate information resides in tree-ring series.

Two models are considered for purposes of distinguishing whether unobserved limiting factors modulating growth through Liebig's law of the minimum are important for variability observed across a tree stand. Following Cook (1985), we first adopt a standard model where tree growth is proportional to the CGF plus a contribution from LGFs,

$$\begin{equation} {{\rm{H}}_0}\,(y,i) = {\rm{CGF}}(y) + {\rm{LGF}}(y,i). \end{equation} \tag{ 1 }$$

The CGF is represented as linearly related to temperature anomalies plus a noise term representing non-temperature processes affecting all trees in a stand as well as inaccuracies in our instrumentally-based estimates of temperature. The LGFs represent any stress on growth that is experienced at a scale smaller than the stand-level. Such local stresses could be induced, for example, by disease and insect infestation; wind, ice and fire damage; soil movement or loss; root damage; variability in the soil microbial community; or changes in the competitive environment experienced by individual trees (White 1979, Sousa 1984). The LGFs are assumed to follow a normal distribution with independent draws in each year and at each tree. If noise contributions were instead assumed to be multiplicative, the logarithm of growth could be used to transform the relationships into a similarly additive formulation.

An alternate model invokes Liebig's Law of the Minimum to postulate competing influences between a CGF and LGFs,

$$\begin{equation} {{\rm{H}}_1}\,(y,i) = \min{\kern-4pt}\left\{ {{\rm{CGF}}(y),\quad {\rm{LGF}}(y,i)} \right\}. \end{equation} \tag{ 2 }$$

Both the CGF and LGFs are parameterized as in H₀, but where a non-zero mean is permitted for the LGFs to allow for the possibility that it is more, or less, dominant in determining growth. More complex representations of tree growth and noise contributions are obviously possible, and potentially worth further exploring, but H₀ and H₁ are adopted to simply illustrate the implications of treating unobserved stresses as additive noise as opposed to Law of the Minimum stresses in the climate reconstruction context.

In order to test whether H₀ or H₁ better represent growth variability, we examine the Schweingruber archive of tree-ring density records (Schweingruber and Briffa 1996). These records are generally interpreted as reflecting large-scale variations in temperature (Briffa et al 2004), and are used instead of width because they more closely track temperature variations (Briffa et al 2002) and therefore better facilitate distinguishing between H₀ and H₁. The association between density and temperature is an empirical one, but experimental evidence suggests that cell wall thickness, the primary determinant of density, is strongly influenced by the amount of time available for the cell wall to thicken, which will on average be extended by warmer ambient conditions (Vaganov et al 2006, Vaganov et al 2011). All tree-ring density chronologies north of 50°N are considered that have a significant positive correlation with local summer temperature (P < 0.05) and at least 20 individual cores contributing to the chronology, giving 207 out of a total 496 initial records. Significance is determined at the P = 0.05 level and is assessed taking into account the full auto-covariance structure of the series by building a null distribution by randomizing the phases of the test series (Schreiber and Schmitz 2000).

Tree-ring series are detrended using a linear age-growth curve, and then normalized before taking percentiles. Similar results are recovered using negative exponential, Hugershoff, and constant detrending, though the linear growth curve appears more appropriate from visual examination of density records.

Temperature is taken as the JJA average computed from 1° × 1° gridded monthly temperatures (Rohde et al 2013), with temperature sampled from the grid closest to the location of each record. Results are not appreciably altered when instead using six-month April–September averages, or using another temperature compilation (Jones et al 1999). Both compilations are based upon averages of daily maximum and minimum temperatures from instrumental weather stations and involve detailed corrections for changes in measurement practices.

In order to simulate tree growth, it is necessary to specify the distributions associated with the CGF and LGFs. The CGF is represented as local temperature plus Gaussian noise having a signal-to-noise ratio P₁. The CGF is normalized to zero mean and unit variance without loss of generality because it is only the relative scaling with respect to the LGFs that determines model behavior. The LGF is assumed normal with a variance P₂ and mean P₃ Note that P₃ has no effect on the interannual variability of the growth resulting under H₀ (equation (1)).

Three observables are used to determine these three parameters. The first observable is the correlation between mean chronologies and instrumental temperature (equal to R = 0.5, when averaged across all sites). This correlation primarily constrains P₁ because noise variance decreases the correlation between temperature and the various percentile chronologies in roughly equal proportion. The second observable is the average correlation between mean chronologies and the individual standardized series at each site (R = 0.65). Only these first two observables are used to constrain H₀ (equation (1)). The third observable, used to constrain growth under H₁ (equation (2)), is the variance across standardized series, calculated separately for each year and then averaged across years and normalized by the interannual variance of the mean chronology (0.56). The second and third observables together constrain P₂ and P₃, which together control how frequently an LGF versus a CGF is expressed.

Optimization in a least-squares sense gives model parameters of P₁ = 0.6 and P₂ = 1.4 under H₀, and P₁ = 0.7, P₂ = 1 and P₃ = 0.3 under H₁. These parameters are well-constrained for our purposes (see supplementary data available at stacks.iop.org/ERL/12/114018/mmedia), and we adopt them for all sites in the following analyses, although we discuss sensitivity to variations further below.

Under H₀, the optimal approach for estimating temperature is to form a chronology as the simple average of all standardized tree-ring records. We are unaware of a proof for optimally inferring the CGF under H₁, but a number of studies have found that quantile methods are appropriate when analyzing observations controlled by a Law of the Minimum (Schröder et al 2005, Austin 2007, Huston 2002, Cade et al 1999). Conceptually, it is preferable to exclude contributions from records that are undergoing relatively low growth, as they are more likely to be responding to an LGF than the CGF, suggesting a preference for higher percentiles when attempting to reconstruct the CGF.

Percentile chronologies are formed by taking a percentile of the collection of growth indicators across each given year. For each site we begin with raw density measurement series for each individual core (Schweingruber and Briffa 1996). By convention, density is taken as the maximum density measured in each ring, though results are essentially identical if the mean latewood density is instead used. After detrending each series, we truncate the tree-ring measurements to the longest period of overlap with local gridded instrumental temperature (Rohde et al 2013), which averages 99 years. Each individual core series is normalized to unit variance. Normalization ensures that all trees, in expectation, are equally represented in the percentile series, and prevents trees with a higher average density from dominating the result. The Nth percentile series is then calculated from these normalized records by choosing the Nth percentile value from amongst all measurements in the given year.

As an illustrative example, we consider 21 cores from the Franchise Wood site in Great Britain (figure 1). In the case of 21 cores, the 11th, 19th and 21st largest value in each year (marked with colored dots, figure 1(a)) correspond with the 50th, 90th and 100th percentile values. (Here we have adopted the convention of calling the maximum value the 100th percentile.) The percentile timeseries are then formed from these individual values (figure 1(b)). The correspondence of percentiles with individual measurements is useful for illustrative purposes, but not necessary. In most cases, percentiles will not correspond exactly to values from an individual core, but are instead interpolated between two cores.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

H₀ and H₁ make different predictions as to the signal ratio of percentile series. The signal ratio (figure 2) for each percentile is calculated as the ratio of the squared pearson's correlation (R²) between local temperature and the percentile time series relative to the R² between local temperature and the mean chronology. It thus represents the ratio of the temperature signal explained by the percentile method relative to the usual averaging method. The signal ratio is used because correlations between tree-ring growth and temperature vary widely (see figure S1). Such variability between stands can be understood in the context of H₁ as a consequence of variations in the relative importance of CGF and LGF expression for controlling growth as well as how dominant temperature is in the CGF, as opposed to other macro-scale factors such as moisture availability (Salzer et al 2014, Bunn et al 2005, Vaganov et al 2006). Synthetic records are generated in order to assess signal ratio as a function of percentile under each model. The number of years and trees simulated are specified in accordance with the data, and all values missing from the data are masked from the synthetic records. For statistical stability, we simulate each site 100 times and report the average results.

Under H₀, the best-performing percentile is the 50th percentile because it most nearly approximates the optimal method of computing the mean (figure 2(a)). In this scenario the expected signal ratio is always less than one and symmetrically tapers toward lower values away from the 50th percentile. Under H₁, in contrast, the signal ratio systematically increases toward higher percentiles on account of better isolating variations caused by the CGF (figure 2(b)). This pattern of monotonic increases in reconstruction skill with increasing percentile holds for H₁ for any set of parameters for which both the CGF and LGFs are expressed. For our selected parameters, percentiles outperform the mean at the 60th percentile and above, with the highest percentiles capturing 15% more temperature variance than the mean (figure 2(b)).

To test whether the data better accords with H₀ or H₁, we next assess signal ratios as a function of percentiles in tree-ring density data. Signal ratios decisively correspond with H₁ by systematically increasing toward higher percentiles and outperforming the average chronology at upper percentiles (figure 2(c)). The upper quartile of percentiles has a higher signal ratio than the lower quartile in 96% of the 207 stands that we examine, and half of the stands have a peak signal ratio above the 87th percentile (figure 3). Chronologies computed using the 91st percentile capture, on average, 20% more temperature variance than the average chronology (figure 2(c)). We thus conclude that unobserved stresses, expressed intermittently as modulated through Liebig's law of the minimum, exert a first-order control on tree-ring density variations between individual trees within a site. Results are essentially identical if the average chronology is instead computed using the bi-weighted mean or if autoregressive standardization (Cook 1985) is applied before analysis, regardless of whether the ARSTAN series or the residual series are used. Results for each of the 207 sites are separately presented in the supplementary data.

A pure application of the law of the minimum predicts that skill will improve monotonically up to the highest percentiles (figure 2(b)), but our empirical analysis indicates that the signal ratio falls off sharply above the 91st percentile (figure 2(c)). This drop off can be understood as resulting from outliers. Outliers introduce a tradeoff whereby higher percentiles better exclude LGFs but make chronologies more susceptible to positive outliers.

We append a stochastic noise term to (equation (2)) to give an extended version of our model, referred to as H₁^outlier, that reproduces the drop off in signal ratio (figure 2(b)). Specifically, outliers are simulated under H₁^outlier by adding an additional noise process n_outlier that operates after the min function. The parameters for H₁^outlier are otherwise identical to H₁. Any reasonable method of specifying outliers produces a roll-off in the signal ratio at the highest percentiles (figure 2(b), blue line), although simply adding gaussian noise after the min function produces a roll-off in skill at high percentiles that is more gradual than what is observed. Here we choose n_outlier to have a value of zero in 99.7% of draws, and to be drawn from a log-normal distribution with μ = 0 and σ = 1 in the other 0.3% of cases.

Other sources of noise are almost certainly also present—due, for example, to measurement error or radial variability within individual rings—but our simulations indicate that a small number of positive outliers well explain the sharp roll-off in skill at the highest percentiles. Potential sources of positive outliers include compression wood, incomplete resin extraction, or moisture influencing density.

There are 16 stands out of 207 where our percentile approach qualitatively does not accord with H₁ (figure 3). These instances comprise 2 stands where the mean chronology outperforms any percentile and another 14 where percentiles below the 50th outperform higher percentiles. These anomalies may represent random chance. However, H₁^outlier predicts such anomalous relations should only occur at 4% of sites. Furthermore, when we run H₁^outlier in a mode where the LGF distribution has a higher mean (changing P₃ from 0.3 to 0.65), decreasing the frequency with which the LGF determines tree growth, the frequency with which the mean or lower-percentile outperform upper percentiles increases to 8%. Two lines of evidence support this revised parameterization for these 16 outlier sites. First, these stands coincide with locations estimated to have a temperature limitation (Nemani et al 2003) that is stronger than the rest of the population (p < 0.05, 2-sample t-test). Second, these 16 stands have an average interseries correlation that is significantly greater than that found for the overall population (${\bar R}=0.48$, p < 0.01, 2-sample t-test). The interseries correlation is the average correlation between each pair of standardized tree-ring series contributing to the record and, in the context of H₁, a higher value indicates more frequent expression of the CGF and thus less influence of unobserved LGFs.

A percentile-based approach also permits recovery of climate signal at sites for which the mean chronology contained no significant climate signal. In the initial analysis, 35 of the sites in the network were excluded solely because no significant correlation was found between the mean chronology and local temperature at the P < 0.05 threshold. If instead of the mean chronology, the 91st percentile chronology is taken as a proxy for the large-scale signal at these sites, the correlation with local temperature increases at 29 of the 35 sites and corresponds to an average increase of 0.09. 11 of these sites become significantly correlated with local temperature (P < 0.05).

Local operation of the law of the minimum implies that optimal reconstruction of climate can depend on the climate state itself. The importance of LGFs will vary according to the time-variable values of the CGF. In the case of temperature limitation, warmer temperatures will decrease the expected fraction of trees that express the CGF. Simulations using H₁ indicate that, on average, 16% of density values represent samples of LGFs, while 84% represent samples of the CGF, but this latter value varies from 91% for years with temperatures in the lower quartile, to 61% for years with temperatures in the upper quartile. It is not possible to distinguish whether an individual density measurement samples LGFs or the CGF, but in aggregate we expect larger variance across records during warm years, when LGFs are preferentially sampled. Conversely, in cool years the CGF is more likely to be sampled.

Zoom In Zoom Out Reset image size

Fritts (1976) noted that chronology standard errors were directly proportional to index values at some sites, noting that this was 'undoubtedly a consequence of the law of limiting factors'. Such a warming-induced shift toward LGFs presents a physiological explanation for the apparent increase in spatial variability during the Medieval Warm Period relative to the Little Ice Age found across tree ring temperature reconstructions (Tingley and Huybers 2015). Under H₁, reconstructions incur greater noise during warm intervals through increased contributions from LGFs, where this greater noise variance is expected to translate into greater variability across reconstructions.

Examination of the coherence between chronologies and local temperature permits identifying the frequencies at which skill most improves (figure 4). Coherence is calculated using the multitaper method (Thomson 1982) with 8 windows. Simulations using H₀ predict that percentile chronologies will have lower coherence with instrumental temperatures relative to that calculated using the mean chronology for all frequencies and for all percentiles (figure 4(a)). In contrast, simulations using H₁ predict a decrease in coherence at all frequencies for low percentiles, and an increase in coherence at all frequencies for high percentiles (figure 4(b)). The observations indicate a frequency dependent pattern similar to that predicted by H₁ and unlike that predicated by H₀. Chronologies derived using between the 65th and 86th percentiles show improved coherence across all frequencies relative to an average chronology (figure 4(c)), with a decrease for some frequencies at the highest percentiles, likely associated with the effects of outliers.

The greatest improvement in coherence occurs at the lowest frequencies when using the highest percentiles. H₁ reproduces this pattern, albeit more weakly, when instrumental temperature observations are prescribed as the CGF (figure 4(b)). More generally, the greatest improvement occurs at those frequencies at which the CGF is most energetic (figure S2) because large-amplitude forcing tends to be susceptible to intermittent loss of signal. Temperature estimates are biased low during warmer intervals on account of increased sampling of the LGFs depressing the overall estimate, but this bias is reduced when using high-percentile chronologies relative to an average chronology.

The operation of the law of the minimum is well-established and comes as no particular surprise from an ecological perspective. Indeed, models of tree-ring growth that explicitly parametrize biological laws invoke the law of the minimum to select between multiple explicitly resolved processes controlling growth (Fritts et al 1991, Shashkin and Vaganov 1993, Vaganov et al 2006, 2011, Evans et al 2006, Anchukaitis et al 2006, Tolwinski-Ward et al 2011). Forest dynamics models (Botkin 1993), plankton ecology models (De Baar 1994, Legović and Cruzado 1997) and carbon cycle models (Farquhar et al 1980, Collatz et al 1991) invoke a similar formulation to select between multiple explicitly-resolved environmental limitations on growth. However, standard formulations of tree-ring growth models assume that the factor that is limiting at a given time is common across all trees in a site. Our results establish that the Law of the Minimum operates differentially across trees within a given stand, and that accounting for these differential growth controls allows for more accurate reconstruction of past CGFs. Rather than a noise-plus-signal model, our results favor a model wherein the CGF is only intermittently present in a given tree. In this case, the CGF is better reconstructed through selecting high percentiles of growth across all trees in a stand.

Most dedroclimatological studies rely on measurements of tree ring width, rather than density, because the width measurement is more widely available than are measurements of density. Tree ring width is, however, a noisier climate proxy. For example, in the dataset considered here, the average squared correlation between the mean chronology and local JJA temperature is >4 times larger for density (${\overline{R^2}} = 0.28$) than for width (${\overline{R^2}} = 0.06$). Generally, the addition of any kind of noise after the law of the minimum operator tends to draw down the expected correlation between high percentiles and the CGF, drawing the observed signal ratio behavior towards the prediction under H₀ (as in figure 2(b), blue line) and making detection of any law-of-the-minimum behavior more difficult. In addition, because Schweingruber was explicitly focussed on density, the site selection method used in the present dataset differs from the typical site selection method for dendroclimatological studies using tree ring width, making this a poor dataset for testing hypotheses about controls on width-based climate reconstructions. A proper test for the implications of law of the minimum variability for climate reconstructions based on tree-ring width will require further consideration of additional datasets.

Although an optimal formulation for purposes of reconstructing climate variability is beyond the scope of this paper, some features important for accurate reconstruction can be highlighted. If we repeat our analysis with a fixed network of sites while varying the number of cores included in the analysis from 1 to 20, we find a monotonic increase in the reconstruction skill with the number of cores at all percentile for any reasonable reconstruction method. The trend in correlation between the 91st percentile-based reconstruction and regional temperature variations (0.008 per number of cores included) is somewhat stronger than the analogous trend for the averaging-based reconstruction (0.006 per number of cores included). The different trends arrise because more reconstruction skill is available with the 91st percentile than with the mean chronology whenever more than about 10 cores are available, but the skill of this 91st percentile method drops to be essentially identical to that of the mean chronology once the number of cores drops below 10. This indicates that the number of records available is an even more important determinant of skill when attempting to leverage the additional reconstruction skill available by accounting for action of the law of the minimum. Analysis of the smaller number of sites with more than 20 cores shows that the trend towards increased signal recovery with increase sample depth continues past 30 cores, albeit with smaller marginal gains per core.

It is expected on the basis of order statistics that the maximum value of a random variable will increase with the number of samples. This will be an important consideration for attempts to reconstruct climate further back in time due to the diminishing number of records in earlier periods. Once the number of samples is sufficiently small such that the high percentiles becomes indistinguishable from the maximum value, loss of records will lead to systematic suppression of the expected value. Furthermore, simulations done by varying the number of cores included in the analysis indicate a depression in the percentile most strongly correlated with temperature at sample depths below 10 cores. Thus there will be a need to re-estimate the optimal percentile as well as the mean and uncertainty associated with the reconstruction as a function of the available records if reconstruction is to be extended into periods with very low sample depth.

Despite the complexity implied by local operation of the law of the minimum and utilization of a percentile-based approach to reconstruction, there are grounds to expect systematic improvement in skill based on the higher correlations and coherence. The indication that high-percentiles better reconstruct low-frequency temperature variability is particularly promising. Difficulty in reconstructing low-frequency temperature variations from tree-ring records is most often ascribed to the need to remove growth trends (Cook et al 1995, Briffa et al 1996), but these results suggest that the greater variability in temperature at low-frequencies also contributes to greater noise and biases at these frequencies, and that reconstruction using high percentiles will improve reconstruction of low-frequency changes in temperature. In future work it will be important to more fully characterize the properties of the unobserved stresses that operate through the Law of the Minimum across different regions, species, and time-intervals. Further developing such ecologically-based models of local tree growth should prove useful for better reconstructing past climates.

ARS was supported by NFS Award AGS-1405053. PH was supported by NSF Award AGS-1304309. We thank Fritz Schweingruber for making his density network data available. We thank Pax Lu Huybers for help with field work that motivated this study.