J. theor. Biol. (1982) 94,909922
Competition
and Skewness in Plantations D.J. GATES
CSIRO Division of Mathematics and Statistics, P.O. Box 1965, Canberra City, A.C.T. 2601, Australia (Received 30 December 1980, and in revisedform 31 August 1981)
Data published recently on the stem diameters in experimental Pinus plantations, show a skewnesswhich is initially zero, first becomes negative and later reverses direction, becomes positive and increases indefinitely. This and other behaviour are explained using a zoneofinfluence model based entirely upon competition between neighbouring trees. Negative skewnesscan be identified with the early stagesof competition when only the largest trees compete. The model also generates bimodal distributions when competition is intense, as observed experimentally in annual plants. Further modes are generated as competition is increased further.
Radiata
1. Introduction
In a recent paper (Borough, Gates & McMurtie, 1980) a study was made of the skewness of DBH (stem diameter at breast height) in experimental Pinus Radiata plantations. The experiment was planted in 1940, and measurements taken in 1949, ‘50, ‘52, ‘54, ‘55, ‘58, ‘60, ‘65, ‘70, ‘75 and ‘80. It comprised plots with neighbouring trees at spacings of 6, 7, 8, 9, 10 and 11 feet. Thus it is possible to study the variation of DBH distribution with both age and spacing. A striking feature of the results in Borough et al. (1980), most pronounced in the plots of intermediate spacing, is the variation in skewness of DBH with age, shown in Fig. 1. Positive skewness is widely observed in the literature (Kayama & Kira, 1956; Rabinowitz, 1979; Burkhart & Strub, 1974; Clutter, 1974), though there is no widely accepted theory which explains it. Negative skewness seems not to have been noted previously as a consistent phenomenon. In Gates & Westcott (1980) it was shown mathematically, under rather mild assumptions, that a negative skewness should be expected in the early stages of competition. Our purpose now is to derive in more detail the skewness variation shown in Fig. 1, using a “zoneofinfluence” model which is an extension of the model described in Gates (1978, 1980a, b), Gates & Westcott (1978) and Gates, O’Connor & Westcott (1979). 909
00225193/82/040909+14$02.00/0
@ 1982 Academic Press Inc. (London) Ltd.
910
D.
.I.
GATES
0.1 
o
0.1

0.2

0.3

0.4

05

1
1945
1950
1955
1960
1965
1970
I975
FIG. 1. Variation of skewness with time in three plots with 11 foot spacing, indicated by three different symbols. The lines are subjective best fit curves for each plot. Time of planting was 1940.
Some new equations of growth are derived in section 2, their method of numerical solution described in section 3 and the results described in section 4. 2. The Growth Equations We introduce now a new set of coupled differential equations, one for each tree, which seem to correctly describe the competition process. The emphasis is on a careful accounting of the sharing of resources (light, water, etc.) among competing trees, because oversimplification of this aspect leads to models which fail to predict the observed skewness behaviour. Details of the distribution and metabolism of resources within individual trees are drastically simplified, and are not meant to be quantitatively accurate. As in many growth models (e.g. Hesketh, Baker & Duncan, 1971; Jones & Smerage, 1978; McMurtie, 1981) we assume that the excess absorbed
COMPETITION
AND
SKEWNESS
IN
PLANTATIONS
911
energy of a tree is converted into weight increment +;=
k(EE,w)
(1)
where W is the total weight at time t, E is the rate of absorbed energy, EM is the rate of energy usage for maintenance, and k is the conversion efficiency, assumed constant, of energy into weight. Thus (1) combines, in a naive way, the growth of the many components of a tree and their different conversion efficiencies. Equation (1) has never been experimentally proved or disproved. We assume that the weight is proportional to height h and basal area 7rD2/4 where D is the DBH, i.e. W = phD2
(2)
where p is constant. This formula is clearly good for the part of the weight, @‘ho’ say, contributed by the stem, and tacitly assumes that the root and crown dimensions are proportional to the stem dimensions. It is supported by various data (e.g. Egunjobi, 1975). We further suppose that h=AD’
(3)
where A and { are constants, giving W a D’+‘. On purely dimensional grounds, one might expect 5 = 1, but values as low as l= O3 are sometimes observed (e.g. Egunjobi, 1975; Attiwill & Ovington, 1968). Let us, for the moment, assume that light is the only form of absorbed energy. Then, by analogy with food crops, we assume E = &A
(4)
where q is the light intensity (independent of tree size), S is the area of the projection of the exposed crown on the horizontal plane, (assuming light is radiating vertically downward) and 6, is the fraction of absorbed light per unit area. It is clear that 4 is larger for deeper crowns, and various simple models imply that q5 is proportional to crown depth. From arguments of scale one concludes that crown depth is proportional to tree height, so that c,b=wh
(3
where w is a constant. This will not hold in dense stands when selfpruning of older branches occurs. Equation (4) has not been tested experimentally for trees: it probably fails for large q since photosynthetic conversion is known to be nonlinear.
912
D.
J.
GATES
A further natural assumption is that, for a young tree, the energy required for maintenance is proportional to its weight EM = YW
(6)
where y is independent of tree size. Possibly this is violated in older trees since heart wood should not require maintenance. Again (6) is untested experimentally. Substituting (2) for (6) in (1) yields dB
(7)
x+d=K’S,
where B = 7rD2/4 is the basal area of the stem, and 2ky a=(2+P)y
(8)
rkwq K’=2(2+5)p
(9)
and
are independent of tree size, though they will have some time dependence, on small time scales, through diurnal, seasonal and random variations in q and y. Equivalently, (7) can be written in the form dA ;+uA=KS,
where A = pB is the area of the disc which is the projection on the plane of the crown of an isolated tree with basal area B, and p is a constant while K = FK’. The disc is called the zone of influence (ZOI) of the tree. We note that (7) involves only the areas B and S while (10) involves only A and S because factors proportional to h have been divided out. The process is therefore essentially one of competing discs on the plane. Since S 5 A in (lo), no plant will grow unless K >U.
(11)
For isolated trees (with no competitors), S is the area of the projection of the complete crown, which, by the scaling argument, is proportional to B. Then (7) and (10) reduce to dB
X=(Ku)B
and
$=(KcT)A,
(12)
COMPETITION
AND
SKEWNESS
IN
913
PLANTATIONS
which give an exponential growth (on a time scale of at least one year for which K is constant). This is commonly observed in isolated plants (Kayama & Kira, 1956). For competing trees, S depends on the size of the tree and of its near neighbours, and on the crown shapes, as indicated in Fig. 2a. The projection of the crowns on the horizontal plane (Fig. 2b) shows that S is determined
lb)
FIG. 2. Competing crowns (a), and their projections (b) on the horizontal plane. The shaded area in (b) is S, for this case.
by the intersection of the discs which are projections of the continuations of the crowns. The projections of the space curves through P and Q (the intersection of the crown surfaces) are the plane curves PIP2 and Q1Q2 respectively, and these latter curves are boundaries of S. Their shape depends on the crown shapes: for example, if the crown surfaces are ellipsoids of revolution differing only in scale (i.e. they are all geometrically “similar”) and having centres on a common horizontal plane, then the curves PIP2 and Q1Q2 are common (straight line) chords of pairs of discs. If water and nutrients are in ample supply, while light is limited, then light is the only resource competed for, and (7) is adequate. When water is competed for, the competition within the crowns will follow an equation like (4) for quantity of water rather than energy. Competition for water and nutrients among root systems of different trees can be modelled in the
914
D.
J.
GATES
same way except that the discs and reduced areas S now correspond to the roots areas. Our scaling assumptions imply that the isolated root area of a tree be proportional to its B and its root area under competition to S. Assuming that (1) applies to nutrients and water, with quantity replacing energy on the right side of the equation, leads again to (7) where u and K now contain new factors. Alternatively, (7) can be derived directly from w to represent an abstract “resource” (1) by taking the energies E and E which combines (in an unspecified manner) the contributions from light, water and nutrients. Now however, the ZOI areas A and reduced areas S should be reinterpreted as areas over which the trees effectively use and compete for this combined resource. We need to specify the regions S for this more general situation. In Gates et al. (1979), it was shown, on the basis of several simple assumptions about scaling and growth, that the boundary curves, between two discs such as PIP2 in Fig. 2, have equations of the form IrxI”
Rz
= Iryla
Rr
(13)
where x and y are the centres of two overlapping discs, R, and R, are their radii, r is the locus point on the boundary, and (Y3 1 is a constant measuring the degree of domination of larger discs over smaller ones, and depends on the species of tree. For example, the overlap area is fairly equally shared if LY= 1, but the larger disc takes all of the overlap if (Y= co. Death of a plant due to competition begins when dA,/dt < 0, so that the plant is losing its region of influence. This begins when S,IUA,/K which implies that the plant receives insufficient resources for maintenance. When dA,/dt falls to uA,, at time T say, then S, = 0, so that the plant receives no resources. Thereafter A,(t)
= A,(T)
ea(rT’
(14)
so that A, decays exponentially to zero. The definition of time of death is a matter for convention (e.g. T or T + l/m). The identification A, = pB. does not apply for dA,/dt < 0. When a tree fails to meet its maintenance requirements, DBH remains fixed, or decreases slightly until death is complete. 3. The Simulation
Format
As in the Radiata experiments, we take the trees to be planted on a square lattice with unit spacing. If the tree at a point x has ZOI are A, then its reduced area &({A{,) depends on A, and on the ZOI areas {A}x of those neighbours of x whose discs are large enough to intersect the ZOI
COMPETITION
AND
SKEWNESS
IN
at x. This results in a set of coupled differential
dA. dt+h
PLANTATIONS
915
equations for the Ax's,
= ~sd{A)x).
(15)
We omit the CTA, terms from the simulations: thus the maintenance requirement of trees is ignored, and describes the situastion where energy supply is far in excess of that required for maintenance. This makes little difference for the phenomena studied here where the trees are not seriously stressed. We take 100 trees in a 10 x 10 grid and use reflection boundary conditions, which means that a mirror image of the plantation is present beyond its natural edges: A x,+10,x2 = A loxI.x2 A x,.x~+lo =A
(16)
x,.10x2
This helps to eliminate edge effects caused by untypical competition faced by trees on the perimerter, and imitates, to some extent a much larger plantation. The ZOI areas A.(O) at time 0 are taken to be independent and identically distributed with the beta density function ?(A:)(:A)
f(A) =
i 0
fori
otherwise.
(17)
Distributions with long tails were avoided because inaccuracies in the tails at later times lead to large inaccuracies in higher moments like the skewness. The range of initial A’s roughly imitates observed data for young plantations, that is, before competition begins. It implies an initial ZOI radius 5 5, so that there is no overlap of discs and no competition at our chosen zero of time. The time dependence in (15) was traced in discrete steps using the difference equations Ax0 + 1) = A.(t) + wSx({AO)L)
(18)
where w is K times the time increment: it determines the timescale of the steps in (18). For the degree of domination (Y the values 1, 2, 4 and co were used thereby spanning the possible range of relative dominance. The simulations were run for 10 discrete times values t = 0, 1,2, . . . ,9. The values w = 0.1 and O2 were both used, these values being compromises which give both a reasonable smoothness of time development and sufficiently advanced development at time t = 9. Results for w = 0.2 case, only, are presented.
916
D. J. GATES
The S,‘s were computed by superimposing a fine square grid of points at a spacing & of the tree spacing, and testing whether each point lies on the appropriate side of each boundary of S,. Thus each S, contained typically about 100 points, when the canopy became closed, and was estimated with about I % error. Using this range of parameters, the B,‘s and S,‘s were computed for the 100 trees and the 10 time values. At each time value the mean, variance, skewness and kurtosis of both the Ax's and S,‘s were computed together with histograms of both variables on 20 equal interval classes. Histograms and moments of the diameters dx = (4A,/#*
= p “‘DBH
were also calculated. 4. Results of the Simulations We present a fairly complete description of the results of the simulation with cx= 2 and w = 0.2, followed by an outline of the other cases, showing the influence of the parameters and their appropriateness for real plantations. The density function (17) implies that the 100 ZOI’s are nonintersecting and have a small range of radii 0.354% R 50.5. Actually the random sample used contained only radii in the range 0.3701 RI 0.485. As t increases the discs grow according to (18). At time t = 9 the plantation area is completely covered by the reduced areas S, as shown in Fig. 3. The
FIG. 3. The reduced zones S. of plants in a 10x 10 simulated plantation at time t = 9. Zones of influence for four of the plants are also shown.
COMPETITION
AND
SKEWNESS
IN
PLANTATIONS
917
boundaries of the S, are just common chords of the discs, four of which are shown in Fig. 3. All discs have grown markedly, the largest has a radius of 1.033. Figure 3 enables us to see at a glance the relative sizes (the S,) of the reduced crowns (or reduced ZOI’s), which by (7) (with c = O), are also the relative growth rates of the basal areas of the trees, i.e.
(19) Figure 3 also illustrates a negative correlation between nearest neighbour S,‘s; i.e. a large S, tends to have small neighbours. This is deduced theoretically in the following papers Gates (1980~) and Gates & Westcott (1980). The step by step growth of the plantation is illustrated in Figs 4(a) and (b) where histograms of both the disc areas A, and reduced areas S, are plotted for times t = 0, 1, . . .9. The variation in skewness with t is very clear for both the Ax's and S,‘s. Both have skewness =O at t = 0, which fall to minima of 2.8 at t = 5 and 8.7 at t = 4 respectively, and rise to positive values by the times t = 8 and t = 6 respectively. It is evident that skewness changes are much more pronounced and somewhat earlier for the S,‘s than for the A,%. This is broadly consistent with the basis of our model which sees the S,‘s (the crowns or roots) as the driving force for the stem growth. To make a direct comparison with the radiata data, the skewness of the DBH
is plotted in Fig. 5. To relate the simulation time t = 0, 1,2, . . . ,9 to the age in Fig. 1 we must compare initial conditions. Our t = 0 in Fig. 5 corresponds to the instant when the largest tree reaches ZOI radius $ on a unit grid (see equation (17)). This is 5$ feet on the llft grid. Seedling size at time of planting seems not to have been recorded so that ZOI size at planting is not known. A reasonable orderofmagnitude estimate is maximum radius 1 foot. In this precompetition stage, the growth equation ( 12) (U = 0) have solution
A(t)=A(O)e"'
(20)
so that Kt = log, {A(t)/A(O)}. Substituting the above values A(0) = 7r ft2 and A(t) = %r ft’ gives Kf
= log,
y
=
3.4.
(21)
918
D.
J.
GATES
..:L1.Ll 12
3
12
34
1234
1234
FIG. 4. Histograms of (a) disc (ZOI) areas A, and (b) reduced areas S. in the simulation attimesf=O,l,... 9. The horizontal axis is graduated in multiples of r/8.
Thus time of planting in Fig. 5 is t = 3.4, Choosing a time scale in the model so that t = 3.4
compared to 1940 in Fig. 1.
corresponds to 10 years
(22)
gives the fitted dates in Fig. 5. Comparing this with experimental plots of skewness in Fig. 1 shows an agreement which is surprisingly good. The numerical agreement is largely a coincidence: The qualitative agreement, on the other hand, is a confirmation of the theory.
COMPETITION
AND
SKEWNESS
IN
PLANTATIONS
919
0.1 . . .
0.2
l
.
0.3

0.4

.
. . .
1
l
1
I
1970 1950 1960 I980 FIG. 5. Variations of skewness of DBH with time t in the simulation with (2 = 2. This may be compared with the data in Fig. 1. Corresponding time of planting is t = 3.4. Fitted dates follow from (22).
The t = 9 histogram of the S.‘s shows the beginning of two distinct peaks. This “bimodality” has been observed in experiments with annual plants (Rabinowitz, 1979; Ford, 1973, and has been derived from various simplified theoretical models Gates (1978), Diggle (1976), Ford & Diggle (1980) and Hesketh et al. (1971). The effect becomes more pronounced when the dominance of larger plants over smaller ones is greater. This is illustrated in Fig. 6, which gives histograms of disc areas for the (Y= 4,
2
4 (bl
6
2
4
6
(Cl
FIG. 6. Histograms of (a) disc areas at time r = 9 in the (I = 4 case (b) disc areas at time t = 7 in the a = COcase and (c) reduced areas at time t = 7 in the (Y= m case. The graduations are again multiples of rr/8.
920
D.
J.
GATES
w = 0.2 case at t = 9 and of disc areas and reduced disc areas for CY= 03, w = 0.2 at t = 7 (the simulation was not run further because this would necessitate taking account of next neighbour interactions). Not only are the two principal modes quite distinct here, but a third mode becomes evident also. Trimodal structure appears not to have been observed: data usually exhibits too much variation to resolve three modes. Possibly also, the degrees of dominance implied by (Y= 4 and (Y= 00 is not usually realized. One can give a plausible explanation of the modal structure. Plants with an initial advantage in size gradually increase this advantage through competition. Later on the separation into dominated and dominating plants is so great that they appear as two distinct modes. Later still the plants in the dominant mode become large enough to compete among themselves, thereby splitting this mode into two distinct modes. In theory this process could continue to produce new modes indefinitely. Bimodality does not occur in the Pinus Radiatu data. It is a phenomenon characteristic of the latter stages of competition when mortality has reached about 75% (Ford, 1975). 5. Conclusion
The model studied here predicts, for the first time, skewness variation that agrees qualitatively with data on Pinus Radiuta. It also predicts bimodal distributions like those observed in annual plant experiments. The model is undoubtedly a gross simplification of real growth and competition, but it does have a biological basis which is at least plausible. No serious attempt at model validation has been attempted. At the modelling level this would require physiological data relevant to (l)(6), and this, for the most part, is not available. At the prediction level, we have only one parameter (Y (describing crown shape) which has a qualitative effect: this was varied through values 1, 2, 4 and 00, the value 2 appearing to give the best qualitative fit. The modelling assumptions (l)(6) could be improved in various ways. In more mature plantations where selfpruning of older branches occurs, C$ in (5) would probably lose its h dependence, that is, all trees might absorb about the same fraction of light per unit crown area. This gives
d Wx= cS,  kE,, dt
(23)
where c = kq4 is constant. Then W, is related to B. as before while EM, is some more complex function. When the canopy is closed, the S, cover
COMPETITION
AND
SKEWNESS
IN
PLANTATIONS
921
all the available area CS,+Na’, x
where N is the number of trees and a is the stem spacing. Defining mean weight by
and mean maintenance
requirement [email protected] rca
,!?M similarly,
the
gives
’  k,!?,w.
If maintenance requirement is a small part of available energy then this predicts that weight increment asymptotes to a constant value, as commonly observed by foresters. The equations (23) probably give predictions qualitatively similar to the studied equations (15). The model can be made biologically more realistic by accounting for different tree components such as crown, stem and root. However, the data considered here appear not to warrant such a complex model. I am indebted to R. McMurtrie and M. Westcott for helpful discussions, to K. Malafant for some preliminary computer programming and to the referees for helpful comments and corrections. REFERENCES A~IWILL, P. M. & OVINGTON, J. D. (1968). Forest Sci. 14, 13, Table 1. BOROUGH, C., GATES, D. J. & MCMURTIE, R. (1980). Skewness reversal of DBH in Pinus Radiuta
plantations (preprint).
BURKHART, H. E. & STRUB, M. R. (1974). A model for simulation of planted loblolly pine stands. In: Growth Models for Tree and Stand Simulation. (Fries, J. ed.) Royal College of Forestry, Stockholm, Res Notes 30, 379 pp. CLUTITR, J. L. (1974). A growth and yield model for Pinus Radiata in New Zealand, in Fries, J. lot. cit. DIGGLE, P. J. (1976). J. appl. Prob. 13, 662. EGUNJOBI, J. K. (1975). Ecol. Plant. 11, 109, Table II. FORD, E. D. (1975). J. Ecol. 63,311. FORD, E. D. & DIGGLE, P. (1980). Ann. Borany. GATES, D. G. (1978). J, rheor. Biol. 72, 525. GATES, D. J. (1980a). Mnrh. Biosci. 48, 157. GATES, D. J. (19806). Math. Biosci. 48, 195. GATES, D. J. (1980~). Math. Biosci. 57. GATES,
D. J. & WESTCOTT,
M.
(1978).
Adu.
Appl.
Prob. 10,499.
922
D.
J.
GATES
GATES, D. J. & WESTCOTT, M. (1980). J. m&h. Bid. GATES, D. J., O’CONNOR, A. J. & WESTCOTT, M. (1979). Proc. Roy. SOC. Lond. HESKETH, J. D., BAKER, D. N. & DUNCAN, W. G. (1971). Crop. Sci. 11,394. JONES, J. W. & SMERAGE, G. H. (1978). ASAE Tech. Paper No. 784024. KAYAMA, H. & KIRA, T. (19.56). J. Itwt. Polytech. Osaka Cy. Ser. B33, 326. MCMURTRIE, R. (1981). J. theor. Biol. 89, 151.
A367,59