Bayesian thermometer

Salekhard (66.53 °N, 66.67 °E) is located in Yamalia at the Northern end of the Ural mountain range. 95% of the daily average temperatures observed at Salekhard weather station since 1882 fall within the shaded area of the above graph.  The red line at 5°C is an indicative lower threshold for plant growth. Clearly, the growing season in Yamalia is temperature limited, variable and short.

Tree ring data have been collected from specimens of living and sub-fossil larch trees (Larix Sibirica) found in the vicinity. This is known as the Polar Ural dataset. Specialists have created growth index chronologies extending back over more than a thousand years, allowing inferences to be made about local climatic history. The data confirm that annual tree growth is correlated with “summer temperature” (usually taken as the average of June-July-August).  The Polar Ural dataset includes X-ray measurements of  annual maximum latewood density (MXD). It turns out that MXD is more strongly correlated with temperature than are noisier tree ring widths, making the Polar Ural data an especially useful tree ring thermometer.

Plants have a more complex response to temperature than does, say, a mercury thermometer. In horticulture daily growth is usually taken to be zero below a base temperature threshold T_B and linear above. However this response is unlikely to be uniform over the growing season. Indeed for most plant species, warmth is more beneficial in the early “vegetative” part of the growing season than in late summer or autumn. Putting these facts together, I assume that annual tree growth is sensitive to an effective growing season temperature S defined by

(1)   \begin{equation*}S(\{T_1,...T_{D}\}) = \frac{1}{D}\sum_{i=1}^{D} w_i \begin{cases} 0 & T_i < T_B \\ T_i-T_B & T_i > T_B \end{cases}\end{equation*}

where T_i is temperature on day i and D is an interval of days which spans the growing season (e.g. for Polar Urals I took D=200 from 11 Apr to 27 Oct).  Daily weights w_i (\sum_i w_i = 1) can be modelled flexibly using a two parameter \mathrm{Beta} function (defined on the interval [0,1]):

(2)   \begin{equation*}w_i  = \mathrm{B}(\frac{i}{D}, \alpha,\alpha + \epsilon) \end{equation*}

\alpha > 1 and \epsilon = 0 produces a centrally peaked distribution, while \epsilon controls the degree of asymmetry. Positive values of \epsilon give a higher weight to earlier parts of the growing season and visa versa.

Although wood density is expected to increase with S, it cannot exceed some physical limit, implying an additional non-linearity in the relationship between S and MXD growth index. A simple model of expected MXD is:

(3)   \begin{equation*} \overline{\mathrm{MXD}} = A(1-e^{-\lambda S})\end{equation*}

For the same reason, errors about expected MXD are skewed to the downside. I assume a skew-Normal distribution:

(4)   \begin{equation*} \mathrm{MXD} \sim \mathcal{N}_{skew} \left(\overline{\mathrm{MXD}},\omega,\mu \right)\end{equation*}

For zero skew (\mu=0) this reduces to a standard normal distribution with mean \overline{\mathrm{MXD}} and standard deviation \omega.

With suitable priors, joint posterior probability distributions for the 7 parameters of Equation 1-4 were obtained using the remarkable Markov Chain Hamiltonian Monte Carlo package stan.

Rendered by

For stan aficionados, the calculation used 4 chains each with 1000 iterations of which 500 were warmup. All \hat{\mathrm{R}} values of <1.01 suggesting good convergence. The fit to MXD data is illustrated below using mean values of the parameters (showing 5% and 50% confidence intervals of the skew Normal error distribution).

A few more things.

  1. The posterior distribution of \alpha and \epsilon can be used to estimate the point in the growing season where growth is most sensitive to temperature. This turns out to be 27 June \pm 17 days. As expected, this is in the earlier part of the growing season.
  2. Prior and posterior distributions of T_B are shown in the graph below. T_B is quite well constrained by the data at biologically reasonable values i.e. within the range expected for a cold weather plant species.
  3. The posterior “summer temperature” S determined by Bayesian analysis is similar but not identical to the standard JJA summer temperature often used by specialists to reconstruct of past climates close to the tree-line in the Northern Hemisphere. The relationship is illustrated below using mean posterior parameters values.


In a follow up post I will look at the calibration problem of reconstructing S from observed \mathrm{MXD}.

Briffa et al 2013


weather GHCN