Environmental Research Letters
ACCEPTED MANUSCRIPT • OPEN ACCESS
Implications of Liebig's law of the minimum for tree-ring reconstructions
of climate
To cite this article before publication: Alexander (Zan) Stine et al 2017 Environ. Res. Lett. in press https://doi.org/10.1088/1748-9326/aa8cd6
Manuscript version: Accepted Manuscript
Accepted Manuscript is “the version of the article accepted for publication including all changes made as a result of the peer review process,
and which may also include the addition to the article by IOP Publishing of a header, an article ID, a cover sheet and/or an ‘Accepted
Manuscript’ watermark, but excluding any other editing, typesetting or other changes made by IOP Publishing and/or its licensors”
This Accepted Manuscript is © 2017 The Author(s). Published by IOP Publishing Ltd.
 
As the Version of Record of this article is going to be / has been published on a gold open access basis under a CC BY 3.0 licence, this Accepted
Manuscript is available for reuse under a CC BY 3.0 licence immediately.
Everyone is permitted to use all or part of the original content in this article, provided that they adhere to all the terms of the licence
https://creativecommons.org/licences/by/3.0
Although reasonable endeavours have been taken to obtain all necessary permissions from third parties to include their copyrighted content
within this article, their full citation and copyright line may not be present in this Accepted Manuscript version. Before using any content from this
article, please refer to the Version of Record on IOPscience once published for full citation and copyright details, as permissions may be required.
All third party content is fully copyright protected and is not published on a gold open access basis under a CC BY licence, unless that is
specifically stated in the figure caption in the Version of Record.
View the article online for updates and enhancements.
This content was downloaded from IP address 96.241.23.20 on 04/11/2017 at 03:50
Page 1 of 28 AUTHOR SUBMITTED MANUSCRIPT - ERL-104091.R1
1
2
3
4 Implications of Liebig’s law of the minimum for tree-ring
5 reconstructions of climate
6
7
8
9 A.R. Stine1 & P. Huybers2
10
11
12
13 1Department of Earth & Climate Sciences, San Francisco State University, San Francisco, Califor-
14
15
16 nia, USA
17
18 2Department of Earth and Planetary Sciences, Harvard University, Cambridge, Massachusetts,
19
20
21 USA
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58 1
59
60
Accepted Manuscript
AUTHOR SUBMITTED MANUSCRIPT - ERL-104091.R1 Page 2 of 28
1
2
3
4 A basic principle of ecology, known as Liebig’s Law of the Minimum, is that plant growth
5
6 reflects the strongest limiting environmental factor. This principle implies that a limiting en-
7
8
9 vironmental factor can be inferred from historical growth and, in dendrochronology, such re-
10
11 construction is generally achieved by averaging collections of standardized tree-ring records.
12
13 Averaging is optimal if growth reflects a single limiting factor and noise but not if growth
14
15
16 also reflects locally variable stresses that intermittently limit growth. In this study a collec-
17
18 tion of Arctic tree ring records is shown to follow scaling relationships that are inconsistent
19
20
21 with the signal-plus-noise model of tree growth but consistent with Liebig’s Law acting at the
22
23 local level. Also consistent with law-of-the-minimum behavior is that reconstructions based
24
25
26 on the least-stressed trees in a given year better-follow variations in temperature than typ-
27
28 ical approaches where all tree-ring records are averaged. Improvements in reconstruction
29
30
31 skill occur across all frequencies, with the greatest increase at the lowest frequencies. More
32
33 comprehensive statistical-ecological models of tree growth may offer further improvement in
34
35
36 reconstruction skill.
37
38
39
40 Introduction
41
42
43
44 The principle that tree growth is determined by the most limiting factor, known as Liebig’s law of
45
46 the minimum (Sprengel 1828, Liebig 1840), underlies efforts at dendroclimatological reconstruc-
47
48
49 tion (Fritts 1976, Fritts & Swetnam 1989, Speer 2010). This principle guides sampling efforts to-
50
51 wards trees that are particularly stressed by a common growth factor (CGF) experienced by all trees
52
53
54 at a site, such as growing-season temperature or moisture availability (LaMarche 1982, Pilcher,
55
56
57
58
59
60
Accepted M
nu cript
Page 3 of 28 AUTHOR SUBMITTED MANUSCRIPT - ERL-104091.R1
1
2
3
4 et al. 1990). As one approaches elevational treeline, for example, trees become smaller, a clear
5
6 visual indication of an upper-limit on growth rates that are often set by cold ambient temperatures
7
8
9 (Körner 2012). It follows that limitation on growth at treeline will generally be relaxed in warmer
10
11 years and permit for greater growth, provided that other, stronger growth limitations are absent.
12
13
14
15 In environments where growth is limited by temperature, the characteristics of annual growth
16
17
rings can serve as a proxy of temperature variation. Dendroclimatic reconstructions of climate rely
18
19
20 on sampling this growth variability, typically by removing a thin radial core of the wood, and mea-
21
22 suring interannual variability in either width (Douglass 1920), density (Parker & Henoch 1971) or
23
24
25 color (McCarroll, et al. 2002). After accounting for the effects of age-growth relationships, dis-
26
27 turbance persistence, and the presence of outliers (Cook 1985, Cook 1992), order ten to a hundred
28
29
30 records are typically averaged together to form a single chronology. Such a mean chronology can
31
32 be optimal if tree growth is controlled by a single CGF with various disturbances superimposed
33
34
35 on individual trees (Cook 1985, Cook 1992), and this approach has proved successful in recon-
36
37 structing annually-resolved variability in temperature and precipitation (e.g. Masson-Delmotte, et
38
39
40 al. 2013).
41
42
43
44 But there exists a tension between averaging across individual records and a local interpre-
45
46 tation of Liebig’s law of the minimum. If Liebig’s law applies at the level of the individual tree,
47
48
49 not all records will always respond to the CGF. Tree growth could instead reflect either a CGF or
50
51 local growth factors (LGFs), depending on which imposes the strongest limitation on a given tree
52
53
54 in a given year (Chapin, et al. 2011). It is common for tree-ring models to invoke Liebig’s law
55
56
57
58
59
60
Accepted Manuscript
AUTHOR SUBMITTED MANUSCRIPT - ERL-104091.R1 Page 4 of 28
1
2
3
4 of the minimum (Fritts, et al. 1991, Shashkin & Vaganov 1993, Vaganov, et al. 2006, Tolwinski-
5
6 Ward, et al. 2011). Such models are typically formulated in terms of rates of growth, and are
7
8
9 thus expressed in terms of Blackman’s law of limiting factors (Blackman 1905). The general-
10
11 ization proposed here posits that unobserved limitations that differ between trees across a site
12
13 causes within-stand variability in limiting factors. Plausibility of such variability is supported by
14
15
16 previous demonstration of environmental sensitivity varying within stands according to details of
17
18 topographic setting (Bunn, et al. 2011, Salzer, et al. 2014, Tran, et al. 2017). In the following,
19
20
21 we present a test for distinguishing between models of tree-ring growth for which all trees at a
22
23 site experience identical forcing, from those for which trees experience multiple limitations upon
24
25
26 growth. We then explore the implications of the latter scenario for where climate information
27
28 resides in tree-ring series.
29
30
31
32
33 Two models of growth factor response
34
35
36
37 Two models are considered for purposes of distinguishing whether unobserved limiting factors
38
39 modulating growth through Liebig’s law of the minimum are important for variability observed
40
41
42 across a tree stand. Following Cook (1985), we first adopt a standard model where tree growth is
43
44 proportional to the CGF plus a contribution from LGFs,
45
46
47 H0(y, i) = CGF(y) + LGFs(y, i). (1)48
49
50
51 The CGF is represented as linearly related to temperature anomalies plus a noise term repre-
52
53 senting non-temperature processes affecting all trees in a stand as well as inaccuracies in our
54
55
56 instrumentally-based estimates of temperature. The LGFs represent any stress on growth that
57
58
59
60
Acc pted Ma us ript
Page 5 of 28 AUTHOR SUBMITTED MANUSCRIPT - ERL-104091.R1
1
2
3
4 is experienced at a scale smaller than the stand-level. Such local stresses could be induced, for
5
6 example, by disease and insect infestation; wind, ice and fire damage; soil movement or loss; root
7
8
9 damage; variability in the soil microbial community; or changes in the competitive environment
10
11 experienced by individual trees (White 1979, Sousa 1984). The LGFs are assumed to follow a
12
13 normal distribution with independent draws in each year and at each tree. If noise contributions
14
15
16 were instead assumed to be multiplicative, the logarithm of growth could be used to transform the
17
18 relationships into a similarly additive formulation.
19
20
21
22 An alternate model invokes Liebig’s Law of the Minimum to postulate competing influences
23
24
25 between a CGF and LGFs,
26
27
28 H1(y, i) = min{CGF(y) , LGFs(y, i)}. (2)29
30
31
32 Both the CGF and LGFs are parameterized as in H0, but where a non-zero mean is permitted for
33
34 the LGFs to allow for the possibility that it is more, or less, dominant in determining growth.
35
36
37 More complex representations of tree growth and noise contributions are obviously possible, and
38
39 potentially worth further exploring, but H0 and H1 are adopted to simply illustrate the implications
40
41
42 of treating unobserved stresses as additive noise as opposed to Law of the Minimum stresses in the
43
44 climate reconstruction context.
45
46
47
48
49 Data and model parameterization
50
51
52
53 In order to test whetherH0 orH1 better represent growth variability, we examine the Schweingruber
54
55 archive of tree-ring density records (Schweingruber & Briffa 1996). These records are generally
56
57
58
59
60
Accep ed Manuscript
AUTHOR SUBMITTED MANUSCRIPT - ERL-104091.R1 Page 6 of 28
1
2
3
4 interpreted as reflecting large-scale variations in temperature (Briffa, et al. 2004), and are used
5
6 instead of width because they more closely track temperature variations (Briffa, et al. 2002) and
7
8
9 therefore better facilitate distinguishing between H0 and H1. The association between density and
10
11 temperature is an empirical one, but experimental evidence suggests that cell wall thickness, the
12
13 primary determinant of density, is strongly influenced by the amount of time available for the
14
15
16 cell wall to thicken, which will on average be extended by warmer ambient conditions (Vaganov
17
18 et al. 2006, Vaganov, et al. 2011). All tree-ring density chronologies north of 50◦N are considered
19
20
21 that have a significant positive correlation with local summer temperature (P <0.05) and at least
22
23 20 individual cores contributing to the chronology, giving 207 out of a total 496 initial records.
24
25
26 Significance is determined at the P=0.05 level and is assessed taking into account the full auto-
27
28 covariance structure of the series by building a null distribution by randomizing the phases of the
29
30
31 test series (Schreiber & Schmitz 2000).
32
33
34
35 Tree-ring series are detrended using a linear age-growth curve, and then normalized before
36
37 taking percentiles. Similar results are recovered using negative exponential, Hugershoff, and con-
38
39
40 stant detrending, though the linear growth curve appears more appropriate from visual examination
41
42 of density records.
43
44
45
46 Temperature is taken as the JJA average computed from 1°x1°gridded monthly temperatures
47
48
49 (Rohde, et al. 2013), with temperature sampled from the grid closest to the location of each record.
50
51 Results are not appreciably altered when instead using six-month Apr-Sep averages, or using an-
52
53
54 other temperature compilation (Jones, et al. 1999). Both compilations are based upon averages
55
56
57
58
59
60
Accepted Manus r pt
Page 7 of 28 AUTHOR SUBMITTED MANUSCRIPT - ERL-104091.R1
1
2
3
4 of daily maximum and minimum temperatures from instrumental weather stations and involve
5
6 detailed corrections for changes in measurement practices.
7
8
9
10 In order to simulate tree growth, it is necessary to specify the distributions associated with
11
12
13 the CGF and LGFs. The CGF is represented as local temperature plus Gaussian noise having a
14
15 signal-to-noise ratio P1. The CGF is normalized to zero mean and unit variance without loss of
16
17
generality because it is only the relative scaling with respect to the LGFs that determines model
18
19
20 behavior. The LGF is assumed normal with a variance P2 and mean P3 Note that P3 has no effect
21
22 on the interannual variability of the growth resulting under H (Eqn. 1).
23 0
24
25
26 Three observables are used to determine these three parameters. The first observable is the
27
28
29 correlation between mean chronologies and instrumental temperature (equal to R=0.5, when av-
30
31 eraged across all sites). This correlation primarily constrains P1 because noise variance decreases32
33
34 the correlation between temperature and the various percentile chronologies in roughly equal pro-
35
36 portion. The second observable is the average correlation between mean chronologies and the
37
38
39 individual standardized series at each site (R=0.65). Only these first two observables are used to
40
41 constrain H0 (Eqn. 1). The third observable, used to constrain growth under H1 (Eqn. 2), is the
42
43
44 variance across standardized series, calculated separately for each year and then averaged across
45
46 years and normalized by the interannual variance of the mean chronology (0.56). The second and
47
48
49 third observables together constrain P2 and P3, which together control how frequently an LGF
50
51 versus a CGF is expressed.
52
53
54
55 Optimization in a least-squares sense gives model parameters of P1=0.6 and P2=1.4 underH0,
56
57
58
59
60
Accepted M
nuscript
AUTHOR SUBMITTED MANUSCRIPT - ERL-104091.R1 Page 8 of 28
1
2
3
4 and P1=0.7, P2=1 and P3=0.3 under H1. These parameters are well-constrained for our purposes
5
6 (see Supplementary Data), and we adopt them for all sites in the following analyses, although we
7
8
9 discuss sensitivity to variations further below.
10
11
12
13 Distinguishing between models
14
15
16
17 Under H0, the optimal approach for estimating temperature is to form a chronology as the simple
18
19
20 average of all standardized tree-ring records. We are unaware of a proof for optimally inferring
21
22 the CGF under H1, but a number of studies have found that quantile methods are appropriate
23
24
25 when analyzing observations controlled by a Law of the Minimum (Schröder, et al. 2005, Austin
26
27 2007, Huston 2002, Cade, et al. 1999). Conceptually, it is preferable to exclude contributions from
28
29
30 records that are undergoing relatively low growth, as they are more likely to be responding to an
31
32 LGF than the CGF, suggesting a preference for higher percentiles when attempting to reconstruct
33
34 the CGF.
35
36
37
38
Percentile chronologies are formed by taking a percentile of the collection of growth indica-
39
40
41 tors across each given year. For each site we begin with raw density measurement series for each
42
43 individual core (Schweingruber & Briffa 1996). By convention, density is taken as the maximum
44
45
46 density measured in each ring, though results are essentially identical if the mean latewood den-
47
48 sity is instead used. After detrending each series, we truncate the tree-ring measurements to the
49
50
51 longest period of overlap with local gridded instrumental temperature (Rohde et al. 2013), which
52
53 averages 99 years. Each individual core series is normalized to unit variance. Normalization en-
54
55
56
57
58
59
60
Acc pted Manuscript
Page 9 of 28 AUTHOR SUBMITTED MANUSCRIPT - ERL-104091.R1
1
2
3
4 sures that all trees, in expectation, are equally represented in the percentile series, and prevents
5
6 trees with a higher average density from dominating the result. The Nth percentile series is then
7
8
9 calculated from these normalized records by choosing the Nth percentile value from amongst all
10
11 measurements in the given year.
12
13
14
15 As an illustrative example, we consider 21 cores from the Franchise Wood site in Great
16
17
Britain (figure 1). In the case of 21 cores, the 11th, 19th and 21st largest value in each year (marked
18
19
20 with colored dots, figure 1(a)) correspond with the 50th, 90th and 100th percentile values. (Here we
21
22 have adopted the convention of calling the maximum value the 100th percentile.) The percentile
23
24
25 timeseries are then formed from these individual values (figure 1(b)). The correspondence of
26
27 percentiles with individual series is useful for illustrative purposes, but not necessary. In most
28
29
30 cases, percentiles will not correspond exactly to values from an individual core, but are instead
31
32 interpolated between two cores.
33
34
35
36 H0 and H1 make different predictions as to the signal ratio of percentile series. The signal
37
38
39 ratio (figure 2) for each percentile is calculated as the ratio of the squared pearson’s correlation
40
41 (R2) between local temperature and the percentile time series relative to the R2 between local tem-
42
43
44 perature and the mean chronology. It thus represents the ratio of the temperature signal explained
45
46 by the percentile method relative to the usual averaging method. The signal ratio is used because
47
48
49 correlations between tree-ring growth and temperature vary widely (see figure S1). Such variabil-
50
51 ity between stands can be understood in the context of H1 as a consequence of variations in the
52
53
54 relative importance of CGF and LGF expression for controlling growth as well as how dominant
55
56
57
58
59
60
Accepted Manuscript
AUTHOR SUBMITTED MANUSCRIPT - ERL-104091.R1 Page 10 of 28
1
2
3
4 temperature is in the CGF, as opposed to other macro-scale factors such as moisture availability
5
6 (Salzer et al. 2014, Bunn, et al. 2005, Vaganov et al. 2006). Synthetic records are generated in
7
8
9 order to assess signal ratio as a function of percentile under each model. The number of years and
10
11 trees simulated are specified in accordance with the data, and all values missing from the data are
12
13 masked from the synthetic records. For statistical stability, we simulate each site 100 times and
14
15
16 report the average results.
17
18
19
20 Under H0, the best-performing percentile is the 50
th percentile because it most nearly ap-
21
22 proximates the optimal method of computing the mean (figure 2(a)). In this scenario the expected
23
24
25 signal ratio is always less than one and symmetrically tapers toward lower values away from the
26
27 50th percentile. Under H1, in contrast, the signal ratio systematically increases toward higher per-28
29
30 centiles on account of better isolating variations caused by the CGF (figure 2(b)). This pattern of
31
32 monotonic increases in reconstruction skill with increasing percentile holds forH1 for any set of pa-
33
34
35 rameters for which both the CGF and LGFs are expressed. For our selected parameters, percentiles
36
37 outperform the mean at the 60th percentile and above, with the highest percentiles capturing 15%
38
39
40 more temperature variance than the mean (figure 2(b)).
41
42
43
44 To test whether the data better accords with H0 or H1, we next assess signal ratios as a
45
46 function of percentiles in tree-ring density data. Signal ratios decisively correspond with H1 by
47
48
49 systematically increasing toward higher percentiles and outperforming the average chronology at
50
51 upper percentiles (figure 2(c)). The upper quartile of percentiles has a higher signal ratio than the
52
53
54 lower quartile in 96% of the 207 stands that we examine, and half of the stands have a peak signal
55
56
57
58
59
60
Accep ed Manuscript
Page 11 of 28 AUTHOR SUBMITTED MANUSCRIPT - ERL-104091.R1
1
2
3
4 ratio above the 87
th percentile (figure 3). Chronologies computed using the 91st percentile capture,
5
6 on average, 20% more temperature variance than the average chronology (figure 2(c)). We thus
7
8
9 conclude that unobserved stresses, expressed intermittently as modulated through Liebig’s law of
10
11 the minimum, exert a first-order control on tree-ring density variations between individual trees
12
13 within a site. Results are essentially identical if the average chronology is instead computed using
14
15
16 the bi-weighted mean or if autoregressive standardization (Cook 1985) is applied before analysis,
17
18 regardless of whether the ARSTAN series or the residual series are used. Results for each of the
19
20
21 207 sites are separately presented in the Supplementary Data.
22
23
24
25 A pure application of the law of the minimum predicts that skill will improve monotonically
26
27 up to the highest percentiles (figure 2(b)), but our empirical analysis indicates that the signal ratio
28
29
30 falls off sharply above the 91
st percentile (figure 2(c)). This drop off can be understood as resulting
31
32 from outliers. Outliers introduce a tradeoff whereby higher percentiles better exclude LGFs but
33
34
35 make chronologies more susceptible to positive outliers.
36
37
38
39 We append a stochastic noise term to Eq. 2 to give an extended version of our model, referred
40
41 to as Houtlier, that reproduces the drop off in signal ratio (figure 2(b)). Specifically, outliers are
42 1
43
44 simulated under H
outlier
1 by adding an additional noise process noutlier that operates after the min
45
46 function. The parameters for H outlier1 are otherwise identical to H1. Any reasonable method of
47
48
49 specifying outliers produces a roll-off in the signal ratio at the highest percentiles (figure 2(b), blue
50
51 line), although simply adding gaussian noise after the min function produces a roll-off in skill at
52
53
54 high percentiles that is more gradual than what is observed. Here we choose noutlier to have a value
55
56
57
58
59
60
Accepted Manuscr pt
AUTHOR SUBMITTED MANUSCRIPT - ERL-104091.R1 Page 12 of 28
1
2
3
4 of zero in 99.7% of draws, and to be drawn from a log-normal distribution with µ = 0 and σ = 1 in
5
6 the other 0.3% of cases.
7
8
9
10 Other sources of noise are almost certainly also present—due, for example, to measurement
11
12
13 error or radial variability within individual rings—but our simulations indicate that a small number
14
15 of positive outliers well explain the sharp roll-off in skill at the highest percentiles. Potential
16
17
sources of positive outliers include compression wood, incomplete resin extraction, or moisture
18
19
20 influencing density.
21
22
23
24 There are 16 stands out of 207 where our percentile approach qualitatively does not accord
25
26 with H1 (figure 3). These instances comprise 2 stands where the mean chronology outperforms27
28
29 any percentile and another 14 where percentiles below the 50th outperform higher percentiles.
30
31 These anomalies may represent random chance. However, Houtlier predicts such anomalous rela-
32 1
33
34 tions should only occur at 4% of sites. Furthermore, when we run Houtlier in a mode where the1
35
36 LGF distribution has a higher mean (changing P3 from 0.3 to 0.65), decreasing the frequency with
37
38
39 which the LGF determines tree growth, the frequency with which the mean or lower-percentile
40
41 outperform upper percentiles increases to 8%. Two lines of evidence support this revised pa-
42
43
44 rameterization for these 16 outlier sites. First, these stands coincide with locations estimated to
45
46 have a temperature limitation (Nemani, et al. 2003) that is stronger than the rest of the population
47
48
49 (p<0.05, 2-sample t-test). Second, these 16 stands have an average interseries correlation that is
50
51 significantly greater than that found for the overall population (R̄ = 0.48, p<0.01, 2-sample t-test).
52
53
54 The interseries correlation is the average correlation between each pair of standardized tree-ring
55
56
57
58
59
60
Ac ep ed Manuscript
Page 13 of 28 AUTHOR SUBMITTED MANUSCRIPT - ERL-104091.R1
1
2
3
4 series contributing to the record and, in the context of H1, a higher value indicates more frequent
5
6 expression of the CGF and thus less influence of unobserved LGFs.
7
8
9
10 A percentile-based approach also permits recovery of climate signal at sites for which the
11
12
13 mean chronology contained no significant climate signal. In the initial analysis, 35 of the sites in
14
15 the network were excluded solely because no significant correlation was found between the mean
16
17
chronology and local temperature at the P<0.05 threshold. If instead of the mean chronology, the
18
19
20 91st percentile chronology is taken as a proxy for the large-scale signal at these sites, the correlation
21
22 with local temperature increases at 29 of the 35 sites and corresponds to an average increase of
23
24
25 0.09. 11 of these sites become significantly correlated with local temperature (P<0.05).
26
27
28
29 Climate influence on reconstruction
30
31
32
33
34 Local operation of the law of the minimum implies that optimal reconstruction of climate can
35
36 depend on the climate state itself. The importance of LGFs will vary according to the time-variable
37
38
values of the CGF. In the case of temperature limitation, warmer temperatures will decrease the
39
40
41 expected fraction of trees that express the CGF. Simulations using H1 indicate that, on average,
42
43 16% of density values represent samples of LGFs, while 84% represent samples of the CGF, but
44
45
46 this latter value varies from 91% for years with temperatures in the lower quartile, to 61% for
47
48 years with temperatures in the upper quartile. It is not possible to distinguish whether an individual
49
50
51 density measurement samples LGFs or the CGF, but in aggregate we expect larger variance across
52
53 records during warm years, when LGFs are preferentially sampled. Conversely, in cool years the
54
55
56
57
58
59
60
Acc pted Manuscript
AUTHOR SUBMITTED MANUSCRIPT - ERL-104091.R1 Page 14 of 28
1
2
3
4 CGF is more likely to be sampled.
5
6
7
8 Fritts (1976) noted that chronology standard errors were directly proportional to index values
9
10 at some sites, noting that this was undoubtedly a consequence of the law of limiting factors. Such a
11
12
13 warming-induced shift toward LGFs presents a physiological explanation for the apparent increase
14
15 in spatial variability during the Medieval Warm Period relative to the Little Ice Age found across
16
17
tree ring temperature reconstructions (Tingley & Huybers 2015). Under H1, reconstructions incur18
19
20 greater noise during warm intervals through increased contributions from LGFs, where this greater
21
22 noise variance is expected to translate into greater variability across reconstructions.
23
24
25
26 Examination of the coherence between chronologies and local temperature permits identi-
27
28
29 fying the frequencies at which skill most improves (figure 4). Coherence is calculated using the
30
31 multitaper method (Thomson 1982) with 8 windows. Simulations using H0 predict that percentile32
33
34 chronologies will have lower coherence with instrumental temperatures relative to that calculated
35
36 using the mean chronology for all frequencies and for all percentiles (figure 4(a)). In contrast,
37
38
39 simulations using H1 predict a decrease in coherence at all frequencies for low percentiles, and an
40
41 increase in coherence at all frequencies for high percentiles (figure 4(b)). The observations indi-
42
43
44 cate a frequency dependent pattern similar to that predicted by H1 and unlike that predicated by
45
46 H0. Chronologies derived using between the 65
th and 86th percentiles show improved coherence
47
48
49 across all frequencies relative to an average chronology (figure 4(c)), with a decrease for some
50
51 frequencies at the highest percentiles, likely associated with the effects of outliers.
52
53
54
55 The greatest improvement in coherence occurs at the lowest frequencies when using the
56
57
58
59
60
Accep ed Manuscript
Page 15 of 28 AUTHOR SUBMITTED MANUSCRIPT - ERL-104091.R1
1
2
3
4 highest percentiles. H1 reproduces this pattern, albeit more weakly, when instrumental temperature
5
6 observations are prescribed as the CGF (figure 4(b)). More generally, the greatest improvement
7
8
9 occurs at those frequencies at which the CGF is most energetic (figure S2) because large-amplitude
10
11 forcing tends to be susceptible to intermittent loss of signal. Temperature estimates are biased low
12
13 during warmer intervals on account of increased sampling of the LGFs depressing the overall
14
15
16 estimate, but this bias is reduced when using high-percentile chronologies relative to an average
17
18 chronology.
19
20
21
22
23 Further discussion and conclusions
24
25
26
27 The operation of the law of the minimum is well-established and comes as no particular surprise
28
29
30 from an ecological perspective. Indeed, models of tree-ring growth that explicitly parametrize bio-
31
32 logical laws invoke the law of the minimum to select between multiple explicitly resolved processes
33
34 controlling growth (Fritts et al. 1991, Shashkin & Vaganov 1993, Vaganov et al. 2006, Vaganov
35
36
37 et al. 2011, Evans, et al. 2006, Anchukaitis, et al. 2006, Tolwinski-Ward et al. 2011). Forest dy-
38
39 namics models (Botkin 1993), plankton ecology models (De Baar 1994, Legović & Cruzado 1997)
40
41
42 and carbon cycle models (Farquhar, et al. 1980, Collatz, et al. 1991) invoke a similar formulation to
43
44 select between multiple explicitly-resolved environmental limitations on growth. However, stan-
45
46
47 dard formulations of tree-ring growth models assume that the factor that is limiting at a given time
48
49 is common across all trees in a site. Our results establish that the Law of the Minimum operates
50
51
52 differentially across trees within a given stand, and that accounting for these differential growth
53
54 controls allows for more accurate reconstruction of past CGFs. Rather than a noise-plus-signal
55
56
57
58
59
60
Accepted Manuscri t
AUTHOR SUBMITTED MANUSCRIPT - ERL-104091.R1 Page 16 of 28
1
2
3
4 model, our results favor a model wherein the CGF is only intermittently present in a given tree. In
5
6 this case, the CGF is better reconstructed through selecting high percentiles of growth across all
7
8
9 trees in a stand.
10
11
12
13 Most dedroclimatological studies rely on measurements of tree ring width, rather than den-
14
15 sity, because the width measurement is more widely available than are measurements of density.
16
17
Tree ring width is, however, a noisier climate proxy. For example, in the dataset considered here,
18
19
20 the average squared correlation between the mean chronology and local JJA temperature is >4
21
22 times larger for density (R2 = 0.28) than for width (R2 = 0.06). Generally, the addition of any
23
24
25 kind of noise after the law of the minimum operator tends to draw down the expected correlation
26
27 between high percentiles and the CGF, drawing the observed signal ratio behavior towards the
28
29
30 prediction under H0 (as in Figure 2b, blue line) and making detection of any law-of-the-minimum
31
32 behavior more difficult. In addition, because Schweingruber was explicitly focussed on density,
33
34
35 the site selection method used in the present dataset differs from the typical site selection method
36
37 for dendroclimatological studies using tree ring width, making this a poor dataset for testing hy-
38
39
40 potheses about controls on width-based climate reconstructions. A proper test for the implications
41
42 of law of the minimum variability for climate reconstructions based on tree-ring width will require
43
44
45 further consideration of additional datasets.
46
47
48
49 Although an optimal formulation for purposes of reconstructing climate variability is beyond
50
51 the scope of this paper, some features important for accurate reconstruction can be highlighted. If
52
53
54 we repeat our analysis with a fixed network of sites while varying the number of cores included
55
56
57
58
59
60
Accepte  Manuscript
Page 17 of 28 AUTHOR SUBMITTED MANUSCRIPT - ERL-104091.R1
1
2
3
4 in the analysis from 1 to 20, we find a monotonic increase in the reconstruction skill with the
5
6 number of cores at all percentile for any reasonable reconstruction method. The trend in correlation
7
8
9 between the 91st percentile-based reconstruction and regional temperature variations (0.008 per
10
11 number of cores included) is somewhat stronger than the analogous trend for the averaging-based
12
13 reconstruction (0.006 per number of cores included). The different trends arrise because more
14
15
16 reconstruction skill is available with the 91st percentile than with the mean chronology whenever
17
18 more than about 10 cores are available, but the skill of this 91st percentile method drops to be
19
20
21 essentially identical to that of the mean chronology once the number of cores drops below 10. This
22
23 indicates that the number of records available is an even more important determinant of skill when
24
25
26 attempting to leverage the additional reconstruction skill available by accounting for action of the
27
28 law of the minimum. Analysis of the smaller number of sites with more than 20 cores shows that
29
30
31 the trend towards increased signal recovery with increase sample depth continues past 30 cores,
32
33 albeit with smaller marginal gains per core.
34
35
36
37 It is expected on the basis of order statistics that the maximum value of a random variable
38
39
40 will increase with the number of samples. This will be an important consideration for attempts
41
42 to reconstruct climate further back in time due to the diminishing number of records in earlier
43
44
45 periods. Once the number of samples is sufficiently small such that the high percentiles becomes
46
47 indistinguishable from the maximum value, loss of records will lead to systematic suppression of
48
49
50 the expected value. Furthermore, simulations done by varying the number of cores included in the
51
52 analysis indicate a depression in the percentile most strongly correlated with temperature at sample
53
54
55 depths below 10 cores. Thus there will be a need to re-estimate the optimal percentile as well as
56
57
58
59
60
Accepted Manus ript
AUTHOR SUBMITTED MANUSCRIPT - ERL-104091.R1 Page 18 of 28
1
2
3
4 the mean and uncertainty associated with the reconstruction as a function of the available records
5
6 if reconstruction is to be extended into periods with very low sample depth.
7
8
9
10 Despite the complexity implied by local operation of the law of the minimum and utilization
11
12
13 of a percentile-based approach to reconstruction, there are grounds to expect systematic improve-
14
15 ment in skill based on the higher correlations and coherence. The indication that high-percentiles
16
17
better reconstruct low-frequency temperature variability is particularly promising. Difficulty in re-
18
19
20 constructing low-frequency temperature variations from tree-ring records is most often ascribed to
21
22 the need to remove growth trends (Cook, et al. 1995, Briffa, et al. 1996), but these results suggest
23
24
25 that the greater variability in temperature at low-frequencies also contributes to greater noise and
26
27 biases at these frequencies, and that reconstruction using high percentiles will improve reconstruc-
28
29
30 tion of low-frequency changes in temperature. In future work it will be important to more fully
31
32 characterize the properties of the unobserved stresses that operate through the Law of the Minimum
33
34
35 across different regions, species, and time-intervals. Further developing such ecologically-based
36
37 models of local tree growth should prove useful for better reconstructing past climates.
38
39
40
41
42 Supplementary Data accompanies this manuscript.
43
44
45
46 Acknowledgements ARSwas supported by NFS Award AGS-1405053. PH was supported by NSF Award
47
48 AGS-1304309. We thank Fritz Schweingruber for making his density network data available.
49
50
51
52 Author Contributions ARS and PH contributed equally to formulating the project and communicating
53
54 its results. ARS performed the calculations.
55
56
57
58
59
60
Accept d Manuscript
Page 19 of 28 AUTHOR SUBMITTED MANUSCRIPT - ERL-104091.R1
1
2
3
4
5
6
7 Competing Interests The authors declare that they have no competing financial interests.
8
9
10
11 Correspondence Correspondence and requests for materials should be addressed to
12
13 ARS (email: stine@sfsu.edu).
14
15
16
17
18 K. Anchukaitis, et al. (2006). ‘Forward modeling of regional scale tree-ring patterns in the south-
19
20
21 eastern United States and the recent influence of summer drought’. Geophys. Res. Lett 33.
22
23
24 M. Austin (2007). ‘Species distribution models and ecological theory: a critical assessment and
25
26
27 some possible new approaches’. Ecological modelling 200(1):1–19.
28
29
30 F. F. Blackman (1905). ‘Optima and limiting factors’. Annals of Botany 19(74):281–295.
31
32
33 D. B. Botkin (1993). Forest dynamics: an ecological model. Oxford University Press on Demand.
34
35
36
37 K. Briffa, et al. (1996). Climatic variations and forcing mechanisms of the last 2000 years, chap.
38
39 Tree-ring variables as proxy-climate indicators: problems with low-frequency signals, pp. 9–42.
40
41
42 Springer.
43
44
45 K. Briffa, et al. (2004). ‘Large-scale temperature inferences from tree rings: a review’. Global and
46
47 Planetary Change 40(1-2):11–26.
48
49
50
51 K. Briffa, et al. (2002). ‘Tree-ring width and density data around the Northern Hemisphere: Part
52
53 1, local and regional climate signals’. The Holocene 12(6):737–757.
54
55
56
57
58
59
60
Accepted Manuscript
AUTHOR SUBMITTED MANUSCRIPT - ERL-104091.R1 Page 20 of 28
1
2
3
4 A. G. Bunn, et al. (2011). ‘Topographically modified tree-ring chronologies as a potential means
5
6 to improve paleoclimate inference’. Climatic Change 105(3):627–634.
7
8
9 A. G. Bunn, et al. (2005). ‘Topographic mediation of growth in high elevation foxtail pine (Pinus
10
11
12 balfouriana Grev. et Balf.) forests in the Sierra Nevada, USA’. Global Ecology and Biogeography
13
14 14(2):103–114.
15
16
17
18 B. S. Cade, et al. (1999). ‘Estimating effects of limiting factors with regression quantiles’. Ecology
19
20 80(1):311–323.
21
22
23 F. S. Chapin, et al. (2011). Principles of terrestrial ecosystem ecology. Springer Science & Busi-
24
25
26 ness Media.
27
28
29 G. J. Collatz, et al. (1991). ‘Physiological and environmental regulation of stomatal conductance,
30
31
32 photosynthesis and transpiration: a model that includes a laminar boundary layer’. Agricultural
33
34 and Forest Meteorology 54(2):107–136.
35
36
37
38 E. Cook (1985). A time series analysis aprproach to tree ring standardization. Ph.D. thesis,
39
40 University of Arizona.
41
42
43 E. Cook (1992). Methods of dendrochronology: Applications in the environmental sciences, chap.
44
45
46 A conceptual linear aggregate model for tree rings, pp. 98–104 pp. Kluwer Academic Publishers.
47
48
49 E. R. Cook, et al. (1995). ‘The ’segment length curse’ in long tree-ring chronology development
50
51
52 for palaeoclimatic studies curse in long tree-ring chronology development for palaeoclimatic
53
54 studies’. The Holocene 5(2):229–237.
55
56
57
58
59
60
Accepted Manuscript
Page 21 of 28 AUTHOR SUBMITTED MANUSCRIPT - ERL-104091.R1
1
2
3
4 H. De Baar (1994). ‘von Liebig’s Law of the Minimum and Plankton Ecology (1899–1991)’.
5
6 Progress in Oceanography 33(4):347–386.
7
8
9 A. E. Douglass (1920). ‘Evidence of climatic effects in the annual rings of trees’. Ecology 1(1):24–
10
11
12 32.
13
14
15 M. Evans, et al. (2006). ‘A forward modeling approach to paleoclimatic interpretation of tree-ring
16
17
18 data’. J. Geophys. Res 111.
19
20
21 G. v. Farquhar, et al. (1980). ‘A biochemical model of photosynthetic CO2 assimilation in leaves
22
23 of C3 species’. Planta 149(1):78–90.24
25
26
27 H. Fritts (1976). Tree rings and climate. London, New York, San Francisco.: Academic Press.
28
29
30 H. Fritts, et al. (1991). ‘Climate variation and tree-ring structure in conifers. Empirical and mecha-
31
32
33 nistic models of tree-ring width, number of cells, cell size, cell-wall thickness and wood density’.
34
35 Clim. Res 1(2):97–116.
36
37
38 H. C. Fritts & T. W. Swetnam (1989). ‘Dendroecology: a tool for evaluating variations in past and
39
40
41 present forest environments’. Advances in ecological research 19:111–188.
42
43
44 M. A. Huston (2002). ‘Introductory essay: critical issues for improving predictions’. Predicting
45
46
47 species occurrences: issues of accuracy and scale pp. 7–21.
48
49
50 P. D. Jones, et al. (1999). ‘Surface Air Temperature and its Changes Over the Past 150 Years’.
51
52
Reviews of Geophysics 37(2):173–199.
53
54
55
56
57
58
59
60
Accepted Manuscript
AUTHOR SUBMITTED MANUSCRIPT - ERL-104091.R1 Page 22 of 28
1
2
3
4 C. Körner (2012). Alpine treelines: functional ecology of the global high elevation tree limits.
5
6 Springer Science & Business Media.
7
8
9 V. LaMarche (1982). Climate from Tree Rings, chap. Sampling Strategies. Cambridge University
10
11
12 Press.
13
14
15 T. Legović & A. Cruzado (1997). ‘A model of phytoplankton growth on multiple nutrients based on
16
17
18 the Michaelis-Menten-Monod uptake, Droop’s growth and Liebig’s law’. Ecological Modelling
19
20 99(1):19–31.
21
22
23 J. Liebig (1840). Organische Chemie in ihrer Anwendung auf Agricultur und Physiologie. English.
24
25
26 Taylor and Walton.
27
28
29 V. Masson-Delmotte, et al. (2013). Climate Change 2013: The Physical Science Basis. Contri-
30
31
32 bution of Working Group I to the Fifth Assessment Report of the Intergovernmental Panel on
33
34 Climate Change, chap. 5. Information from paleoclimate archives, pp. 383–464. Cambridge
35
36
37 University Press.
38
39
40 D. McCarroll, et al. (2002). ‘Blue reflectance provides a surrogate for latewood density of high-
41
42
latitude pine tree rings’. Arctic, antarctic, and alpine research pp. 450–453.
43
44
45
46 R. Nemani, et al. (2003). ‘Climate-driven increases in global terrestrial net primary production
47
48 from 1982 to 1999’. Science 300(5625):1560–1563.
49
50
51
52 M. Parker &W. Henoch (1971). ‘The use of Engelmann spruce latewood density for dendrochrono-
53
54 logical purposes’. Canadian Journal of Forest Research 1(2):90–98.
55
56
57
58
59
60
Accepted Manuscript
Page 23 of 28 AUTHOR SUBMITTED MANUSCRIPT - ERL-104091.R1
1
2
3
4 J. Pilcher, et al. (1990). ‘Primary data’. In Methods of dendrochronology, pp. 23–96. Springer.
5
6
7 R. Rohde, et al. (2013). ‘Berkeley earth temperature averaging process’. Geoinfor Geostat: An
8
9 Overview 1(2):1–13.
10
11
12
13 M. W. Salzer, et al. (2014). ‘Changing climate response in near-treeline bristlecone pine with
14
15 elevation and aspect’. Environmental Research Letters 9(11):114007.
16
17
18 T. Schreiber & A. Schmitz (2000). ‘Surrogate time series’. Physica D 142:346–382.
19
20
21
22 H. K. Schröder, et al. (2005). ‘Rejecting the mean: Estimating the response of fen plant species
23
24 to environmental factors by non-linear quantile regression’. Journal of Vegetation Science
25
26
27 16(4):373–382.
28
29
30 F. Schweingruber & K. Briffa (1996). ‘Tree-ring density networks for climate reconstruction’.
31
32 Climate Variations and Forcing Mechanisms of the Last 2000 Years pp. 43–66.
33
34
35
36 A. Shashkin & E. Vaganov (1993). ‘Simulation Model of Climatically Determined Variability of
37
38 Conifer’s Annual Increment (on the Example of Common Pine in the Steppe Zone)’. Russian
39
40
41 Journal of Ecology 24(5):275–280.
42
43
44 W. P. Sousa (1984). ‘The Role of Disturbance in Natural Communities’. Annual Review of Ecology
45
46 and Systematics 15:353–391.
47
48
49
50 J. H. Speer (2010). Fundamentals of tree-ring research. University of Arizona Press.
51
52
53 C. Sprengel (1828). ‘Von den Substanzen der Ackerkrume und des Untergrundes’. Journal fur
54
55
56 Tecnische und Okonomische Chemie 2:423–474.
57
58
59
60
Accepted Manuscript
AUTHOR SUBMITTED MANUSCRIPT - ERL-104091.R1 Page 24 of 28
1
2
3
4 D. J. Thomson (1982). ‘Spectrum estimation and harmonic analysis’. Proceedings of the IEEE
5
6 70(9):1055–1096.
7
8
9 M. P. Tingley & P. Huybers (2015). ‘Heterogeneous warming of Northern Hemisphere surface tem-
10
11
12 peratures over the last 1200 years’. Journal of Geophysical Research: Atmospheres 120(9):4040–
13
14 4056.
15
16
17
18 S. E. Tolwinski-Ward, et al. (2011). ‘An efficient forward model of the climate controls on inter-
19
20 annual variation in tree-ring width’. Climate Dynamics 36(11):2419–2439.
21
22
23 T. J. Tran, et al. (2017). ‘Cluster analysis and topoclimate modeling to examine bristlecone pine
24
25
26 tree-ring growth signals in the Great Basin, USA’. Environmental Research Letters 12(1):014007.
27
28
29 E. A. Vaganov, et al. (2011). ‘How well understood are the processes that create dendroclimatic
30
31
32 records? A mechanistic model of the climatic control on conifer tree-ring growth dynamics’. In
33
34 Dendroclimatology, pp. 37–75. Springer.
35
36
37
38 E. A. Vaganov, et al. (2006). Growth dynamics of conifer tree rings: Images of past and future
39
40 environments, vol. 183 of Ecological studies. Springer.
41
42
43 P. S. White (1979). ‘Pattern, process, and natural disturbance in vegetation’. The botanical review
44
45
46 45(3):229–299.
47
48
49
50
51
52
53
54
55
56
57
58
59
60
Accepted Ma uscript
Page 25 of 28 AUTHOR SUBMITTED MANUSCRIPT - ERL-104091.R1
1
2
3
4
5
6 (a)
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
1900 1920 1940 1960
28
(b) year
29
30
31
32
1900 1920 1940 1960
33 year
34 (c)
35
36
37
38 1900 1920 1940 1960 1980
39 year
40 Figure 1 Illustration of method for calculating percentile series: Demonstration of how
41
42
43 percentile series are formed at one of the 207 sites used in this study. (a) Individual density series
44
45 from21 cores from Franchise Wood, Great Britain. Magenta, blue and green dots respectively
46
47
48 indicate the core that represents the 50
th, 90th and 100th percentiles for each year. (b) Percentile
49
50 series formed from the individual core series shown in (a). Magenta, blue and green lines re-
51
52
53 spectively represent the 50
th, 90th, and 100th percentile series. (c) Mean chronology for the site,
54
55 formed by taking the average of all density series for each year.
56
57
58
59
60
Acce t mean chronology percentile series normalized single-core density seriesd Manuscript
AUTHOR SUBMITTED MANUSCRIPT - ERL-104091.R1 Page 26 of 28
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
Accepted Manuscript
Page 27 of 28 AUTHOR SUBMITTED MANUSCRIPT - ERL-104091.R1
1
2
3
4
5  180° E 100
6
° W   157 0 0 °  15 E
8 90
°
9  60  N  
10
11 80
12
13
14  75° N  70
15
16 60
17
18
19  90° N  50
20
21
22 40
23
24
25 30
26
27 20
28
29
30 10
31   30 ° °  E 
32 W   3
0
33 0
   0°
34
35
36
37
38
39
40 Figure 3 Percentile with strongest temperature signal. Circles give locations of tree-ring
41
42
43 sites used in this study. Color indicates the percentile chronology having the highest correlation
44
45 with local temperature for the 207 sites where the highest percentile correlation with tempera-
46
47
48 ture exceeds the correlation between temperature and the mean chronology. Black X’s located
49
50 in central Russia indicate the 2 sites for which the mean chronology has a higher correlation with
51
52
53 temperature than any percentile. The close proximity of the two sites make the two X’s indistin-
54
55 guishable.
56
57
58
59
60
Accep  120 ° Wted Manus
  60 ° E c
  90° E ri
maximum correlation percentile pt
°
0  
E
 12
°  W
  60
  90° W
AUTHOR SUBMITTED MANUSCRIPT - ERL-104091.R1 Page 28 of 28
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
Accepted Manuscript