FIG. A1. Gap model schematic. Parameters are
shown in red (see Tables A2 and A3 for values, source and sensitivity
of the parameters). Green shows allometric calculations based on
Table A1. Model driving data is shown in blue (see Fig. A2).
Stochastic processes are shown with a dashed line. The model produces
annual output, but the light and C assimilation functions run on
hourly time steps (see text).

The model is a 1D patch model, explicitly representing each stem
and its light environment, carbon assimilation, growth and mortality.
Each model patch has twenty-five 1 m vertical layers and is not
spatially explicit. The model is written in FORTRAN 90.

Light and seasonal and diurnal model drivers

Within each patch, the model simulates light interception by the
canopy and the resulting photosynthetic assimilation. The model is
driven with a 3.5 year radiation climatology from a weather station
(Skye Instruments, UK) at Chitengo, 25 km from the study site. We
generated a mean diurnal cycle of downwelling PAR, at hourly
resolution, for each calendar month. The model is also driven with an
observed monthly fraction of peak leaf area index (LAI), determined
from 3.5 years of monthly hemispherical photos at two nearby
permanent sample plots.

FIG. A2. Data used to drive the model. A
seasonal cycle of leaf area index was derived from replicated monthly
hemispherical photos on two 1 ha permanent sample plots. PAR
measurements are from a weather station 25 km away. 3.5 years of data
were averaged by hour and month to generate a representative diurnal
cycle for each month.

Plant Area Index (PAI ) was estimated from hemispherical photos
(Nikon Coolpix 4500 with a FC-E8 fisheye converter) collected each
month at two 1 ha permanent sample plots (PSP). Nine photos were
acquired on a 20 m grid on each plot, each month. The photos were
thresholded to separate plant and sky pixels using the algorithm of
Ridler and Calvard (1978), implemented in MATLAB to determine gap
fraction, which was converted to an PAI estimate using the method of
Licor (1989), with the images segmented in a similar manner to van
Gardingen et al. (1999). As the trees are fully deciduous, PAI was
converted to LAI by subtracting the lowest observed PAI from each
observation point. The monthly LAI fraction represents the
deciduousness of the vegetation and is an explicit representation of
limits of the growing season of trees in this environment (Fig. A2).
Each year the vertical leaf area distribution is calculated based on
functions which relate each stem’s DBH to its total height,
canopy depth and leaf area (Shinozaki et al. 1964) using the
allometric relationships shown in Table A1. Leaves were assumed to be
spread evenly throughout the canopy depth. A simple radiative
transfer scheme passes light down through the canopy layers. The
amount of PAR absorbed is calculated according to leaf area density
using the Beer-Lambert law (Jones 1992), with an assumed spherical
leaf angle distribution (k). All radiation was assumed diffuse and
the leaves had no albedo or transmittance.

Carbon assimilation and allocation

The absorbed PAR was converted to assimilated C using
photosynthetic light response curves parameterised from two studies
in Zimbabwean miombo (Tuohy and Choinski 1990; Tuohy et al. 1991).
Two parameters describe the light response curve: the maximum rate of
assimilation, Pmax and the amount of light needed to achieve half
this rate, kp. For each stem, mass of photosynthate was summed for
each canopy layer, for each hour of the 12 diurnal cycles
representative of each month, and scaled up to a yearly total.

Carbon allocation routines are resolved each year for each stem.
We assumed that 50% of assimilate is used for respiration (Waring et
al. 1998, Ra in Table A2). Of the remainder, the required amount is
allocated to leaves to produce the observed phytomass of peak LAI.
Leaf carbon specific area (LCA gC m2) is derived from field
measurements (Nottingham 2004). We assumed that allocation to fine
roots, Afr, is equal to allocation to leaf mass. Any remaining C is
allocated to woody growth, which is partitioned above and belowground
according to a local DBH-specific root:shoot allometric. The increase
in aboveground stem C is converted to a change in DBH using
allometric equations (Table A1), which in turn allows the increase in
height and leaf area to be calculated for the next year.

Mortality and the effect of fire

Once growth is completed, each stem is exposed to a random chance
of mortality, which is either an intrinsic (non-fire) mortality or a
fire-induced mortality derived from our fire experiments. These
fire-induced mortality rates are DBH- and fire intensity
class-specific and are based on piecewise functions of the field
results (Table A3). The intrinsic mortality rate, Mi, was unknown and
is very poorly constrained in this ecosystem. We used a nominal value
of 2% regardless of size (similar to Desanker and Prentice 1994;
Desanker 1996), and undertook sensitivity analysis (Table A2). In the
model, the chance of a fire occurring was determined from a user
supplied control variable, the fire return interval (FRI), and a
random number generator. The chance of a fire occurring was
determined as 1/FRI.

Resprouting

Because of the importance of resprouting in fire-prone ecosystems
(Bond and Midgley 2001; Chidumayo 2004; Mlambo and Mapaure 2006),
aboveground stem mortality (top-kill) was decoupled from belowground
rootstock mortality. Stems > 2 cm DBH (parameter Sresprout) which
are killed, have a probability or resprouting (1 - Smort), which is
based on data from the fire experiments, and we assume they achieve a
2 cm DBH in their first year. In addition, every year there is a
nominal 3% chance (parameter Precruit) of a recruitment event
occurring and 200 new seedlings (parameter Snew) being established in
the patch. We justify this low level of recruitment from studies
which show that most regeneration is from sprouts and that seedling
survival is low in these woodlands (Chidumayo 1992, 1997).

The model was initiated from bare ground with 200 seedlings and
tended to reach equilibrium after ~200 years. All model runs in this
study used 50 patches of 0.02 ha (based on an 8 m crown radius for a
typical large stem) and were 1000 years long. Output for the final
500 years is reported as the 50-patch 500-year mean with the standard
deviation between the 50 patches.

Model sensitivity analysis

Method

We ran a sensitivity test of all model parameters. Each parameter
was varied in turn by a factor of 0.5, 0.75, 1.5, and 2 from its
nominal value (Table A2). The sensitivity, S, is defined as Sx =
([Radj - Rn]/ Rn) / ([Padj - Pn]/ Pn), where x is the factor by which
the nominal value is changed, Radj is the response for the model run
with the adjusted value, Rn is the response with the nominal value,
and Pn and Padj are the parameter values for the nominal and adjusted
cases respectively. Patch basal area was used as the response
variable. The parameters used in the model and their source,
magnitude and sensitivity, are shown in Table A2.

Results

The model output were much more sensitive to the growth parameters
compared to those controlling regeneration (Table A2). In order of
sensitivity, they were: the maximum rate of photo synthesis (Pmax);
the fraction of assimilate allocated to respiration (Ra); the
half-saturation rate of photosynthesis (kp); carbon mass per unit
leaf area (LCA); allocation to fine roots (Afr); the leaf angle
distribution (k); and the leaf area:basal area ratio (LA:BA). These
are basic physiological parameters, which are well constrained by
local field observations (Pmax, and kp, by Tuohy et al. 1991; LCA and
LA:BA by our data) or to a lesser degree, globally (Ra, Waring et al.
1998 and k, Norman and Campbell 1989). Afr is badly constrained,
particularly in savanna environments. Of the mortality and
regeneration parameters, the intrinsic mortality rate (Mi) was the
only one that significantly affected the final patch basal area.

TABLE A1. Allometric equations used in the
model and derived from local field data. Equation forms are defined
as follows: power, Y = aXb; linear, Y = aX ± b. n is the number of
samples of field data used to fit the function.

Dependent (Y) variable

Independent (X) variable

Form

a

b

n

r2

min X

max X

Notes

C stock of stem (kg C)

DBH, m

power

4222

2.6

29

0.93

0.05

0.72

C stem : C root+stem

DBH, m

linear

0.32

0.6

23

0.26

0.05

0.72

Height of top of tree (m)

DBH, m

linear, with saturation

42.6

80

0.67

0.07

0.86

where DBH > 60 cm, ht = 25 m

Height of base of canopy (m)

DBH, m

linear, with saturation

22.3

80

0.58

0.07

0.86

where DBH > 67 cm, canopy base = 15 m

Leaf area, m2

Basal area, m2

linear

1330

10

0.29

4.4

13.9

plot level data

TABLE A2. Parameters in the model, their
source, nominal value and sensitivity, Sx, see text for symbol
definitions. For the sensitivity analysis, we use the 50-patch mean
basal area as the response variable. Each variable was changed in
turn.

Parameter description

Parameter name

Nominal parameter value,
Pn

S0.5

S0.75

S1.5

S2

Source of nominal parameter value

Growth parameters

Fraction of GPP used for autotrophic respiration

Ra

0.5

-1.8

-1.9

-2.4

-1.0

Waring et al. (1998)

Extinction coefficient for Beer-Lambert law

k

0.5

0.0

-0.6

-0.2

0.1

Norman and Campbell (1989)

Amount of C allocated to fine roots, as a fraction of
allocation to leaves

Afr

1

-0.4

-0.2

-1.0

-0.6

Consistent with
Castellanos et al. (2001). Similar to Table 5 of
Hendricks
et al. (2006), which is for temperate systems.

Maximum rate of photosynthesis (µmol C·m-2·s-1)

Pmax

10

2.0

2.0

2.5

2.3

Tuohy and Choinski (1990); Tuohy
et al. (1991)

PAR intensity at which 0.5 Pmax is obtained (µmol·s-1·m-2)

kp

250

-1.7

-1.7

-1.1

-0.7

Tuohy and Choinski (1990); Tuohy et al. (1991)

Leaf carbon per leaf area (g C/m2)

LCA

50

-1.5

-1.1

-1.7

-0.9

Nottingham (2004);
Chidumayo (1997)

Leaf area per unit basal area (m2/m2)

LA:BA

1330

0.0

-0.6

-0.1

-0.1

See Table A1

Mortality and Regeneration parameters

Intrinsic mortality rate

Mi

0.02

-0.8

-1.3

-1.3

-0.5

Estimated, similar to
Desanker and Prentice (1994)

Diameter, m, at which a seedling is considered to have
developed a rootstock and the ability to resprout

Sresprout

0.02

0.0

-0.2

-0.1

0.0

estimated

Number of seedlings per ha established in a recruitment
year

Snew

100

0.0

-0.5

0.2

0.1

estimated

Probability of a recruitment year occurring

Precruit

0.03

0.1

0.1

-0.2

0.0

estimated

Probability of a root stock failing to resprout after
fire

Smort

0.04

-0.2

-0.3

0.1

0.0

this study

TABLE A3. Piecewise parameterisation of
fire-induced mortality rates shown in Fig. A3. These are used to
calculate stem mortality in the model in years when a fire
occurs.

Fire intensity (TA tercile)

High

Med

Low

ax

-70

-70

-60

bx

5

4

2

bx_sat

-2

-3

-3.8

Where:

DBH < 10 cm

Log odds of top-kill = -ax DBH+ bx

DBH > 10 cm

Log odds of top-kill = bx_sat

TABLE A4. Metrological conditions during each burn. The
values are the mean for the period of the fire.

Plot No.

1

2

3

4

5

6

7

8

Temp (°C)

no burn

31.4

29.9

35.4

25.0

33.3

32.8

33.3

Relative Humidity (%)

33

32

23

55

15

28

21

Wind speed (m/s)
-1

2.4

1.4

0.9

1.1

1.7

1.6

0.9

Where:

DBH < 10 cm

Log odds of top-kill = -ax DBH+ bx

DBH > 10 cm

Log odds of top-kill = bx_sat

TABLE A5. Metrological conditions during each
burn. The values are the mean for the period of the fire.

FIG. A3. Model output of patch characteristics
for a low-fire scenario (FRI = 100, low-intensity fires) from
initiation with bare ground at year zero to 1000 years. The mean of
the 50 patches is shown with a black line. One standard deviation of
the 50 patches is shown above and below the mean with a gray shade.
The white triangles on the right Y axis indicate the field data for
the Marrameu forest plot (Site 2), the gray triangles indicate the
50-patch mean of the final 500 model years.