Tree growth analysis issues, again, part four

Continuing….

Last time, I derived an equation that gives constant tree mass growth rates with time, given an apriori, defined relationship between biomass and diameter. Varying two parameters thereof (q and s) produced continuously increasing or decreasing mass and radial (diameter) growth rates, but only the continuously increasing mass growth rates also had reasonable radial growth rates over the full range of small and large trees. I now want to investigate whether there are radial growth models that can give more complex relationships of mass to time, unimodal responses in particular. But before I get to that I’ll show some graphs illustrating these dynamics, and provide the R code as well (at the end).

I’ve altered two things from the previous post, for convenience. First, I re-defined q so as to make the connection between M = f(D), and D = f(t), easier to follow. Specifically,

$M = de^{-1.864 + 2.608\ln(D)}$ and
$ln(D) = \frac{\ln(st)}{2.608q}$ and therefore
$M = de^{-1.864 + \frac{\ln(st)}{q}}$

where M is tree mass, D is diameter, d is wood density, e is the base of natural logs, t is time in years and s and q are computed parameters that jointly give equal mean growth rates at year t = 200 (as explained in the previous post). In this formulation, any model having q < 1.0 will result in continuously increasing mass growth rates as functions of either time or diameter. Second, I changed the target diameter from 156 cm to 120 cm, at 200 years (i.e. slower mean growth rate). This slower mean growth rate just allows me to extend growth beyond 200 years and not get unrealistically sized trees as quickly, at say 300 years.

But before going on, below are four graphs, each showing the growth dynamics of six models, given what I’ve done so far. All produce increasing mass growth rates as a function of time, result from varying q and s in the above equation, and have q < 1.0. The graphs show: (1) diameter and (2) mass, as functions of time, and then mass growth rates as functions of (3) time and (4) diameter. The values of q and s listed in the legends are for the third equation above $M = de^{-1.864 + \frac{\ln(st)}{q}}$

But there’s nothing to say that mass growth rates have to continually increase or decrease with time of course. In fact, there are outstanding arguments as to why they should not, or cannot do so. In the next post in the series I’ll explore the implications of other radial growth models on mass growth rates, and on the conclusions of Stephenson et al. (2014).

R code:

## Create constant, or constantly increasing or decreasing, mass growth rates as f(time) ##
# Using Chave's simplified "moist forest" equation here (linear model, ln transform), for pedagogy
# Chave et al. (2005) Oecologia 145:87–99 DOI 10.1007/s00442-005-0100-x

density = 0.6						# g/cm^3; value not critical, assumed constant over time
q = seq(from=0.8, to=0.3, by= -0.1)			# 1.0 gives constant mass growth rate, <1.0 gives accelerating rate, as f(t)
t.ref = 200; d.ref = 120				# reference point at which equal mean growth rates across models is true
t = 1:((3/2)*t.ref)					# take to t = 300 years to see diverging dynamics
s = exp(q*2.608*log(d.ref))/t.ref
pars = data.frame(q,s,t.ref,d.ref)
n.mods = nrow(pars)
rownames(pars) = paste("model",c(1:n.mods), sep="")

d = rep(NA,length(t)); r = d; w = d; m = d		# templates
res = array(dim=c(length(t),n.mods,3))

for(i in 1:n.mods){
diameter = exp(log(s[i]*t)/(q[i]*2.608))		# Diameter = f(t), derived analytically
# w = c(r[1],diff(r))+1; d=2*cumsum(w)			# for setting minimum allowed ring widths
mass = density * exp(-1.864 + 2.608*log(diameter)) 	# from Chave et al (2005)
res[,i,1] = diameter
res[,i,2] = width
res[,i,3] = mass
}

## Plots
mods.plot = n.mods
labels = c("time (yr)", "ring width (mm)", "diameter (cm)", "mass (kg)", "mass growth rate (kg/yr)", "diameter growth rate (cm/yr)")
plotcolors = palette()[1:mods.plot]
lgnd.text = paste("q=",pars[,1], ", s=", round(pars[,2],2),sep="");lgnd.text

jpeg("S14_1.jpeg",quality=100,width=800,height=800)
plot(t,res[,1,1], xlab=labels[1],ylab=labels[3],ylim=c(0,max(res[,,1])))
for (j in 1:mods.plot) points(t,res[,j,1], col = plotcolors[j])
legend("topleft", legend=lgnd.text, fill = plotcolors, cex=1.5)
dev.off()

jpeg("S14_2.jpeg",quality=100,width=800,height=800)
plot(t,res[,1,3], xlab=labels[1],ylab=labels[4],ylim=c(0,max(res[,,3])))
for (j in 1:mods.plot) points(t,res[,j,3], col = plotcolors[j])
legend("topleft", legend=lgnd.text, fill = plotcolors, cex=1.5)
dev.off()

jpeg("S14_3.jpeg",quality=100,width=800,height=800)
plot(t[-1],diff(res[,1,3]), xlab=labels[1],ylab=labels[5],ylim=c(0,max(diff(res[,,3]))))
for (j in 1:mods.plot) points(t[-1],diff(res[,j,3]), col = plotcolors[j])
legend("topleft", legend=lgnd.text, fill = plotcolors, cex=1.5)
dev.off()

jpeg("S14_4.jpeg",quality=100,width=800,height=800)
plot(res[-1,1,1], diff(res[,1,3]), xlab=labels[3], ylab=labels[5],
ylim=c(0,max(diff(res[,,3]))), xlim = c(0,max(res[-1,,1])) )
for (j in 1:mods.plot) points(res[-1,j,1],diff(res[,j,3]), col = plotcolors[j])
legend("topleft", legend=lgnd.text, fill = plotcolors, cex=1.5)
dev.off()