Simulation of stand dynamics for miombo woodlands in Kitulangalo Forest Reserve, Morogoro, Tanzania

ABSTRACT

Miombo woodlands constitute the most extensive type of forest cover within South and East Africa, Tanzania in particular, and they are seriously threatened by intensive degradation and deforestation. The dynamics of miombo stands were examined in this study using a simple simulation system applied to experimental stands at the Kitulangalo Forest Reserve (KFR) in Tanzania. Prediction of tree-level diameter increment, mortality, and ingrowth across species groups with empirical models based on data from experimental sites constituted the essential elements of the simulation system. Diameter distribution, stand basal area, volume yield, and above-ground biomass accumulation were assessed in three harvesting alternatives across three species groups. The results indicate a high degree of positive response to careful selective harvesting in lightly degraded, recovering miombo forests. Simulation systems compiled from experimentally based models seem promising as tools for research on stand dynamics and sustainable management alternatives in miombo woodlands.

KEYWORDS:

• Simulation
• stand dynamics
• diameter increment
• harvesting alternatives
• miombo woodlands

Introduction

Miombo woodland is the most extensive tropical seasonal woodland and dry forest formation in Africa, which is dominated by the legume family (Fabaceae, subfamily Caesalpinioideae)Frost (1996). Miombo woodlands cover about 270 million ha, extending from Central, Eastern, to Southern African regions (Campbell et al. 2007). Miombo woodlands are a crucial source of a wide range of ecosystems, and environmental services for over 75 million people in the rural settings of Africa, and about 25 million people in urban and suburban areas in the region (Syampungani 2009). The climatic condition of Miombo woodlands is arid and semi-arid, with highly leached soils with poor nutrients where tree growth is relatively slower compared to forests in wetter areas (Backéus et al. 2006; Malmer 2007). Another general distinguishing characteristic of miombo vegetation is the co-dominance of ecosystems of continuous herbaceous cover and sedges (Osborne et al. 2018) marking evident seasonal variations that are related to water stress (Malmer 2007). Unlike the wetter miombo woodlands in the Angolan, Southern miombo, and Zambezian ecoregions receiving over 1000 mm of annual precipitation, the Tanzanian miombo woodlands fall in drier areas receiving less than 1000 mm of rain annually. The Tanzanian miombo occupies about 92% of the forested land, forming the main vegetation type in the country (Mauya et al. 2015; Ene et al. 2017). Generally, the floristic composition and structural variability are so enormous across the landscape (Godlee et al. 2020); with higher species diversity and dimensions in the southern, and northern highlands and the western wet plateaus (Malmer 2007). The central plateaus are dominated by a mosaic of miombo species and other numerous woodland species such as Combretum, Acacia, and Commiphora species (Rajala et al. 2022).

The capacity of the Tanzanian miombo woodlands to provide various ecosystem services is diminishing rapidly due to the escalating anthropogenic activities (Gordon and Enfors 2008) driven by agriculture, illegal logging, wood extraction for poles, bioenergy (charcoal and firewood), unsustainable livestock grazing, and uncontrolled annual fires (Ribeiro et al. 2015a; Doggart et al. 2020). Although forest disturbances that involve opening up the canopy through selective harvesting have been shown to activate regeneration, Chidumayo (1987), Piiroinen et al. (2008) urged that regeneration in miombo woodlands may not be fast and abundant enough to cover the losses when anthropogenic disturbances grow beyond the threshold levels, thus making sustainability nearly impossible. The limited and fragmented information on regeneration, tree growth and the overall stand dynamics in miombo woodlands as a result of habitat modifications by human disturbances and climate change constrains the process of sustainable management of these woodlands, and thus the ecosystem services for the well-being of the local communities. The need for whole-stand simulation systems that can predict the stand development under different management and disturbance regimes is therefore essential for resource assessment planning, and management of miombo woodlands for improved growth, biodiversity values and ecosystem functioning. It is therefore, so essential to understanding the key ecological processes and how miombo woodlands, which is the largest Savannah type in sub-Saharan Africa (Campbell 1996; Frost 1996) respond to both local and global changes.

Unlike temperate and boreal forests, where tree growth and stand development have been extensively studied (Kellomäki et al. 1992; Bugmann 1996; Ruckstuhl et al. 2008), there has been a scarcity of simulation and matrix models for tropical forests and woodlands (Fearnside 2004; Meir et al. 2015; Stanturf and Mansourian 2020). One of the possible reasons is its immense and complex floristic composition and structural variability (Wright 2002; Draper et al. 2020). However, in recent years, there has been a growing concern about the modelling and simulation of tropical forests and woodlands (Bonal et al. 2016). Currently, the allometric models for the estimation of stock variables for individual, dominant tree species in miombo woodlands are readily available (Mugasha et al. 2013, 2016; Hofstad and Araya 2015; Naesset et al. 2016; Manyanda et al. 2020), yet, the yield estimation and overall population dynamic models are much more uncertain (Hofstad and Araya 2015), constraining the forest planners and managers who would like to estimate allowable cuts or non-decreasing harvest levels in miombo woodlands. It is also problematic for the assessment of stand productivity, deforestation, and degradation applicable for REDD+ Projects, which has become essential for the analysis of climate change mitigation (Kapinga et al. 2018; Ribeiro et al. 2021)

Since the application of growth rings is somehow limited and also irregular in miombo tree species (Trouet et al. 2001; Grundy 2006; Syampungani et al. 2010), the only practical method for studying tree growth and population dynamics essential for growth modelling is through re-measurement of trees at Permanent Sampling Plots (PSPs). The exploratory analyses on tree growth patterns by Valkonen et al. (2008) and Elifuraha et al. (2008) using one-year growth data from Kitulangalo Forest Reserve (KFR) PSPs indicated that an empirical modelling approach for miombo conditions needs a medium to long-term growth period and the appropriate methods capable of accounting for the complex inherent and human-induced variability within the stands and dealing with the variation in characteristics, roles, and mutual interaction of individual trees and regeneration (Mbwambo et al. 2008).

The modern methods of combining distance-independent variables, mortality, area-based recruitment models and other allometric models have enabled the formulation of simulation models (Vanclay 1995) to simulate stand development (Bollandsås et al. 2008; Mugasha et al. 2017) and, the same time, predicting diameter growth (Mugasha et al. 2017). The need for long time monitoring of forest growth with an appropriate representation of climatic conditions, is, however, important for the development of reliable and meaningful simulation models for miombo woodlands.

The current study explores the dynamics of miombo stands using the nine-year interval (2007 to 2016) datasets from Kitulangalo study plots in Morogoro, Tanzania, as a study area. The overall objective was to develop a simple and biologically sound whole stand simulation system that would predict the stand dynamics for Kitulangalo miombo trees over a medium-term period of 99 years. The simulation system was modified from Bollandsås et al. (2008) and Yahya et al. (2012); by integrating a diameter increment model (Njoghomi et al. 2021), and the density-dependent multiplicator into the system for addressing the following research questions.

1. Is it feasible to predict and explore the stand dynamics in miombo woodlands with a simulation system composed of basic models based on a relatively medium-term monitoring period on experimental plots at KFR?

2. Are the miombo stands developing towards recovery from past exploitation in the medium term without human disturbances?

3. How would harvesting affect stand dynamics and development as compared to no disturbance? Does such harvesting seem silviculturally feasible? What kind of differences would emerge in wood yields and above-ground biomass carbon stores between the treatment alternatives?

4. What are the implications of harvesting alternatives for research and forest management?

Material and methods

Experimental setting

The study was based on data from Permanent Sample Plots (PSPs) at Kitulangalo Forest Reserve (KFR) in Morogoro, Tanzania. KFR represents the typical lowland dry miombo woodland in Eastern Tanzania. KFR is currently considered as regenerating forest from the past abuse of various disturbances (Malimbwi et al. 2000; Karenget JA-s 2009) which happened before Participatory Forest Management reforms (Lund and Nielsen 2005; Blomley 2006). The high level of accessibility to the highway made this area a prime charcoal production area for the supply of the nearby Morogoro municipality and Dar es Salaam city (Malimbwi et al. 2000). The KFR suffered from timber extraction through the activities of local pit-sawyers, and from cutting of tree stems for building poles (Zahabu 2006) because of insufficient human resources of the Forest Department to maintain control over the area and to prevent the over use of forest resources. This was a typical de facto open access resource (Mbilinyi et al. 2007; Skutsch and Solis 2011). The empirical basis of this study (experimental setting and setup, stand structure and dynamics including species composition, tree size distribution, tree growth, mortality, and recruitment) is described comprehensively and only a few key features are revisited here (Figure 1) (Njoghomi et al. 2021, 2021).

Figure 1. The study area map: The Kitulangalo Forest Reserve in the Morogoro Region, Tanzania.

Two sets of experimental plots had been established in two subjectively selected stands (approximately 500 m apart) to get a proper representation of the miombo woodlands in KFR. Each stand had two blocks with a total of six adjoining PSPs measuring 30 × 30 m. The first measurement was carried out in 2007 and the final measurement in 2016, making nine (9) years intervals. Measurements involved diameter at breast height (Dbh at 1.3 m from the ground), species name, height (for sample trees), and coordinates (x, y) for each tree. Randomization of treatments using nested plots, blocks and stands to study tree growth, mortality and ingrowth (Njoghomi et al. 2021) aimed to reduce bias in estimating mortality in the plots. Because of the scarcity of empirical data and on long-term studies in the Tanzanian miombo woodlands, the use of the small plots (30 × 30 m) was objectively intended to test the applicability of the design in miombo conditions for research purposes.

Data preparation

All trees with a minimum diameter at breast height (Dbh) of 5 cm were measured for the diameter. A tree forking above the Dbh point (1.3 m from the ground level) was considered a single tree; otherwise, all stems from a multiple-stemmed tree were treated individually. The observed tree increment (i_d9) was calculated as the difference between two measurements with an interval of 9 years. Trees that were present and alive in both measurements were included only. The diameter of the trees for all measured trees was categorized into 11 diameter classes, with class size (w) of 5 cm as follows from <5–9.9 cm, 10–14.9 cm, 15–19.9 cm, 20–24.9 cm, 25–29.9 cm, 30–34.9 cm, 35–49.9 cm, 50–54.9 cm, 55–59.9 cm, 60–64.9 cm & 65–69.9 cm. The designation of species grouping followed that of Njoghomi, et al. (2021) where we assigned all tree species to three groups: 1) trees that grow relatively rapidly to be dominants in top canopy layers such as Julbernardia globiflora, Brachystegia boehmii, and Sclerocaria birrea; 2) trees that stay mainly in the lower and middle canopy levels dominated by Combretum species; and 3) trees that grow slowly but persistently and may eventually rise to dominant and codominant canopy positions such as Senegalia nigrescens (syn. Acacia nigrescens) and Vachelia robusta (syn. Acacia robusta). Species group 1 had 170 observations, while groups 2 and 3 had 157 and 81 observations, respectively. The observed number of trees per species group was used for the calculations of different primary variables which were later converted into per ha basis. Individual tree volume was calculated using Mauya et al. (2014) (1a), while above-ground biomass (AGB), and associated carbon were calculated using a height-free allometric model by Mugasha et al. (2013) (1b and 1c)

(1a) $V_{i,j,t}=0.0049+0.0009db{h}_{i,j,t}$(1a)

(1b) $AG{B}_{i,j,t}={\beta }_{0}db{h}_{i,j,t}^{{\beta }_i}\hspace{0.17em}$(1b)

(1c) $AG{C}_{i,j,t}=0.47AGB$(1c)

‘Where = $V_{i,j,t}$, $AG{B}_{i,j,t},$and $AG{C}_{i,j,t}$ are the individual tree volume, above-ground biomass and carbon, respectively, for a tree in species group i, diameter class j at period t. and $dbh{h}_{i,j,t}=$ stem diameter for species group i, diameter class j at period t.

Overview of the simulation system

A simulation system is a common paradigm for analysing complex processes by creating a simplified representation of the system under study (Altiok and Melamed 2010). Simulations require the use of matrix models. In practical forestry, these systems are built to provide managers with predictions of forest growth and yield response to treatments (Yahya et al. 2012). We developed a simulation system to predict stand dynamics as a response to management treatments (harvesting intensities) across three species groups and size classes that are assumed to exhibit similar growth. The developed simulation system is based on a stand-level Dbh distribution matrix model that reflects every single effect of management option (harvesting intensity) that were theoretically applied in the system. A three-level diameter growth model was developed for estimating growth rates in the stand-level matrix model whereby trees with high growth rates upgrow into a higher Dbh class, the slow-growing trees live and stay in the same diameter class, or a tree dies (mortality), a tree is harvested or retained at the beginning of the period, or a tree emerges into the smallest Dbh class during the respective period as a recruit according to Markov chain rule (Acevedo et al. 1995). The density-dependent matrix model (2), therefore, predicts the number of trees per hectare by diameter classes for each species group as influenced by the treatment alternatives. The matrix is highly sensitive to the number of recruits (vector R) added to the smallest diameter class (Dbh 5–10 cm) and the number of trees dying in each diameter class in all species groups.

(2) $y_{t+p}=G\left(y_t-h_t\right)+R$(2)

(Bollandsås et al. 2008)

where $y_{t+p}$= is the state vector, in which each entry is the number of live trees per ha for species group i and diameter class j at time t, G is a transition matrix with a submatrix for each species group,$h_t=$is the harvesting vector in which each entry is the number of trees harvested, and p is the projection interval whereas R is the in-growth (recruitment) vector showing the number of recruits for tree species i entering the lowest diameter class during time t+p.

Other parameters such as stand basal area, volume, above-ground biomass and carbon were simultaneously estimated using the projected number of stems for each harvesting alternative across each species group.

In addition, the systems applied three harvesting alternatives, 1) “No harvesting” whereby trees develop and die with no disturbance, 2) “Uniform intensity” whereby all stems with DBH ≥40.0 cm for all species groups were harvested and finally, 3) “Varying intensity” whereby all stems with DBH ≥24 cm (species group 2) and Dbh ≥40 cm (group 1 and 3) were harvested, respectively. Note that, the simulation system applied a harvesting cycle of 36 years (with four, 9-year time intervals) implemented at the 36^th and 72^nd years under the assumption that potentially harvestable timber volume is calculated only from trees already above the legal minimum harvestable diameter at breast height (Dbh). For example, a prescribed legal minimum harvestable diameter for Julbernardia globiflora is 40 cm, and 45 cm for other dominant tree species such as Brachystegia and Pterocarpus species. However, for other species in the middle and lower canopy levels such as Dalbergia melanoxylon, a prescribed minimum harvestable diameter is 24 cm (Lovett 2003). Although the individual trees in species group 2 are limited in diameter growth and rarely exceed a Dbh of 30 cm, yet they have been harvested even at smaller diameters because of their wide range of uses.

Similarly, illegal timber and pole harvesting are usually targeted at minimum tree sizes and species which are of economic importance, which does not practically differ from what is legally harvested.

Stand-level Dbh distribution matrix model

A matrix model framework modified from Bollandsås et al. (2008) and (Yahya et al. 2012) was applied for the simulation (3): The basic stand table was composed of the number of live trees in each dbh class (5 cm class width) by species groups.

(3) $\left[\begin{array}{c}{y}_{i,1,t+p}\\ y_{i,2,t+p}\\ y_{i,3,t+p}\\ ⋮\\ y_{i,j,t+p}\\ ⋮\\ y_{i,n,t+p}\end{array}\right]=\left[\begin{array}{ccccccc}{a}_{i,1,t}& \hspace{0.17em}& \hspace{0.17em}& \hspace{0.17em}& \hspace{0.17em}& \hspace{0.17em}& \hspace{0.17em}\\ b_{i,1,t}& a_{i,2,t}& \hspace{0.17em}& \hspace{0.17em}& \hspace{0.17em}& \hspace{0.17em}& \hspace{0.17em}\\ \hspace{0.17em}& b_{i,2,t}& a_{i,3,t}& \hspace{0.17em}& \hspace{0.17em}& \hspace{0.17em}& \hspace{0.17em}\\ \hspace{0.17em}& \hspace{0.17em}& \ddots & \ddots & \hspace{0.17em}& \hspace{0.17em}& \hspace{0.17em}\\ \hspace{0.17em}& \hspace{0.17em}& \hspace{0.17em}& b_{i,j-1,t}& a_{i,j,t}& \hspace{0.17em}& \hspace{0.17em}\\ \hspace{0.17em}& \hspace{0.17em}& \hspace{0.17em}& \hspace{0.17em}& \ddots & \ddots & \hspace{0.17em}\\ \hspace{0.17em}& \hspace{0.17em}& \hspace{0.17em}& \hspace{0.17em}& \hspace{0.17em}& b_{i,n-1}& a_{i,n,t}\end{array}\right]\left[\begin{array}{c}{y}_{i,1,t}-h_{i,1,t}\\ y_{i,2,t}-h_{i,2,t}\\ y_{i,3,t}-h_{i,3,t}\\ ⋮\\ y_{i,j,t}-h_{i,j,t}\\ ⋮\\ y_{i,n,t}-y_{i,n,t}\end{array}\right]+\left[\begin{array}{c}{r}_{i,t}\\ 0\\ 0\\ ⋮\\ 0\\ ⋮\\ 0\end{array}\right]$(3)

$y_{i,j,t}$ is the number of alive trees of species i, Dbh class j at time t, n is the number of Dbh classes, and p is the projection interval (9 years).

$h_{i,j,t}$ is a number of harvested trees of species i, Dbh class j at time t.

$a_{i,j,t}$is a proportion of trees staying the same Dbh class of species i, Dbh class j at time t.

$b_{ij,t}$is a proportion of trees upgrowing to higher Dbh class of species i, Dbh class j at time t.

$r_{i,t}$is a number of recruitment trees of species i at time t.

Elements of the Dbh matrix model

The proportions for tree movement into upper diameter classes varied with individual species growth rate, stand basal area and density, mortality, the recruitment and harvesting alternatives. The percentage of trees in species group i that live and stay in the same diameter class j from period t to t+p was denoted by $a_{ij}.$ The predicted number of stems upgrowing to higher diameter classes ($b_{ij})$and the mortality rates ($m_{ij})$ were used to calculate the proportion of trees staying in the same diameter class, expressed as

(4) $a_{i,j,t}=\left(1-b_{i,j,t}-m_{i,j,t}\right)\ast 100$(4)

where

$m_{i,j,t}$is a mortality rate (proportion of the number of stems) in species group i and Dbh class j at the final measurement.

The likelihood that a tree grows from diameter class j to j + 1 between time t and t +9 years at a given harvesting alternative and species group i was estimated as;

(5) $b_{i,j,t}=\frac{i{d}_{i,j,t}}{w}$(5)

Where

$b_{i,j,t}$ = the number of stems for species group i upgrowing to higher diameter classes

The diameter growth model predicted individual increments for species group i and dbh class j in each harvesting alternative after every 9 years. The levels of the diameter increment curves were, however, regulated by the stand basal area in each species group as shown by the coefficient ${\beta }_2$in the model (3). Generally, the anticipated higher basal area values from “no harvesting” and “uniform harvesting” were expected to reduce diameter increments irrespective of species’ groups and diameter class when compared to the “varying alternative”.

Mortality and recruitment sub-models

Different candidate models to predict the probability of tree mortality to happen for species groups 1, 2 and 3 were developed and tested. Unlike mortality and growth models, we tried the formulation and testing of the two-step, area-based recruitment models in which we first predicted the probability of recruitment to happen (in plots with recruits) and then predict the conditional number of recruits using our plot data.

However, because of the insufficient number of observations for both mortality and recruitment, models were discarded because they did not perform. Therefore, since it was not feasible to develop mortality and recruitment models, we applied a simple linear multiplicator developed from the probabilities observed in the data (in each Dbh class across species groups) to make mortality data density-dependent. The mortality rate, therefore, involved the proportion of all dead stems (at the final measurement) to the initial total number of livingstems as shown in Table 1.

Table 1. The observed number of stem mortality and their probability rates.

A. Observed mortality				B. Calculated probabilities
Dbh_class	Species Group1	Species Group2	Species Group3	Species Group1	Species Group2	Species Group3
5.0–10.0	2.0	9.0	5.0	0.033	0.148	0.082
10–15.0	4.0	6.0	5.0	0.066	0.098	0.082
15–20.0	3.0	5.0	1.0	0.049	0.082	0.016
20–25.0	2.0	3.0	2.0	0.033	0.049	0.033
25–30.0	2.0	1.0	1.0	0.033	0.016	0.016
30–35.0	2.0	1.0	0.0	0.033	0.016	0.000
35–40.0	1.0	0.0	0.0	0.016	0.000	0.000
40–45.0	2.0	0.0	0.0	0.033	0.000	0.000
45–50.0	1.0	0.0	0.0	0.016	0.000	0.000
50–55.0	1.0	0.0	0.0	0.016	0.000	0.000
55–60.0	0.0	0.0	0.0	0.000	0.000	0.000
60–65.0	1.0	0.0	0.0	0.016	0.000	0.000
65–70.0	0.0	0.0	0.0	0.000	0.000	0.000

The density-dependent mortality rate was therefore estimated using the equation form:

(6) $m_{i,j,t}=\frac{O{d}_{i,j}}{{\sum }^O{y}_{i,j}}\ast G_{prop,t}$(6)

$O{d}_{i,j}$= the number of dead stems for species group i and diameter class j during the final census,

$O{y}_{i,j}$ = initial total number of living stems for species group i and diameter class j at the end of the study

$Gpro{p}_t=$ proportion of the predicted stand basal area at time t to the stand basal area at time t = 0

A similar, linear multiplicator was used to convert the observed number of recruits into a density-dependent recruitment. The number of recruits was therefore estimated as the observed newly emerged stems into the smallest diameter class (Dbh = 5–10 cm) for each species group i at the start of the experiment conditioned by predicted stand basal area at time t (7). The recruitment vector ($r_i$), therefore, consisted of one sub-vector for each species group with the number of recruits (in-growth) of a specific species i that entered the smallest diameter class from time t to t+p:

(7) $r_{i,t}=n{r}_i\ast G_{prop,t}$(7)

$n_r=$ the number of observed recruit trees of species i at the final census.

The tree-level diameter growth model

This study used a tree-level, nonlinear mixed effect diameter growth model to predict a 9-year diameter increment between 2007 and 2016. The model was adjusted from that of Njoghomi et al. (2021) who used a similar dataset with a measurement interval of 8 years (2008 to 2016). It was worth using the 9-year interval because the 2007 dataset had baseline information for other parameters such as mortality and recruitment which were crucial for simulation purposes.

The change in diameter growth was considered as a response variable while initial stem diameter, species groups and the density-independent stand basal area were considered explanatory variables (fixed effects of the model). The experimental hierarchy (plots/stands) provided random effects of the model. Upgrowth movement of diameters was therefore predicted with diameter increment models (8) constructed with our measurement data with a 9-year interval. An aggregated version of the model (Njoghomi et al. 2021) (8) was developed and applied, where the distance-dependent local basal area was replaced with a distance-independent stand basal area (G). The main point was to eliminate a possible substantial source of bias in the recursive calculation framework where a small systematic error may lead to a large prediction error after several simulation steps.

The diameter growth model was:

(8) $i{d}_{ijk}={\beta }_0+{\mathrm{d}\mathrm{b}\mathrm{h}}_{ij}^{{\beta }_1}\ast e^{\left({\beta }_{2\mathrm{o}}+{\beta }_{2i}\ast \mathrm{S}{\mathrm{P}}_{\mathrm{i}\mathrm{j}}\right)}\ast \mathrm{d}\mathrm{b}{\mathrm{h}}_{ij}{\ast \mathrm{G}}_{ky}^{{\beta }_3}+{\mathrm{u}}_{ij+}{\mathrm{e}}_{ij}$(8)

Where.

$i{d}_{i,j,k}$= 9-year diameter increment of tree k in species group i and diameter class j,

$db{h}_{i,j}$ = initial stem diameter of a tree in species group i and diameter class j,

SP = species groups

$G_{ky}$= Stand basal area

${\beta }_2=$Coefficient for the factor of species groups with three levels (three species

groups) entered through a sub-model $({\beta }_{2o}+{\beta }_{2i}\ast S{P}_{i,j})$.

${\beta }_0\dots .{\beta }_3$= Other coefficients for intercept, initial stems diameter, and plot basal area (the fixed variables)

${\mathrm{u}}_{i,j}$ = random effect of a tree in plot i and diameter class j

${\mathrm{e}}_{i,j}$ = random residual for a tree in plot i and diameter class j

Results

The tree-level diameter growth model

Except for the model intercept, all other model parameters were found to be significant (p < 0.05). The model explained the sources of variation in diameter increment by 38% (R² = 0.38) (Table 2). Residual and density plots (Supplementary 2), however, showed a normal distribution, which indicated the accuracy and reliability of the model in predicting the diameter increments across species groups and diameter classes as influenced by harvesting alternatives. The lower SD values of the model showed that most of the diameter increments values were very close to the mean. The negative model coefficients$,{\beta }_2$indicated that the diameter increment was negatively correlated to the initial diameter of the stems for all species groups (Groups 1, 2 & 3).

Table 2. Parameters for diameter increment model (p < 0.05).

Model expression with predictors	Parameter	Coefficient	STD	P-value	R^2	RSME
$i{d}_{ijk}={\beta }_0+{\mathrm{d}\mathrm{b}\mathrm{h}}_{ij}^{{\beta }_1}\ast e^{\left({\beta }_{2\mathrm{o}}+{\beta }_{2i}\ast \mathrm{S}{\mathrm{P}}_{\mathrm{i}\mathrm{j}}\right)}\ast \mathrm{d}\mathrm{b}{\mathrm{h}}_{ij}{\ast \mathrm{G}}_{ky}^{{\beta }_3}+{\mathrm{u}}_{ij+}{\mathrm{e}}_{ij}$	${\beta }_0$	0.111	0.084	0.073	0.38	0.1541
	${\beta }_1$(dbh)	2.145	0.330	0.0000
	${\beta }_{2\hspace{0.17em}}$(SP1)	−0.125	0.022	0.0000
	${\beta }_2($SP2)	−0.254	0.017	0.0000
	${\beta }_2\hspace{0.17em}$(SP3)	−0.264	0.026	0.0000
	${\beta }_3\left(\mathrm{G}\right)$	−0.217	0.184	0.003

In addition, the increase in plot basal area (${\beta }_3)$, however, decreased with the diameter increments. The coefficient ${\beta }_2$ adjusted the variation in growth levels between species groups. The higher coefficient value (${\beta }_2)$for species group indicated higher growth potentials for species group 1 than other species groups with lower values of coefficient ${\beta }_2$. The model output demonstrated in Figure 2 indicates a single line of prediction for each species group for the initial plot basal area. A detailed model output with multiple prediction lines (t + 0 to t + 99) (Supplementary 1) showed how stand developed through basal area growth, which was a prime factor to regulate diameter increment and hence stand development. The top line in each curve represented a diameter growth under lower stand basal, while the lowest curves, species group and harvesting alternative, however, represented diameter increment and higher stand basal area values.

Figure 2. Diameter increment curves by three species groups.

Diameter distribution and species proportion

A great majority (92%) of the initial density was concentrated in the lower diameter classes (Dbh 5–25 cm) at the beginning and the end of the 99 years in all treatment alternatives. Species group 1 dominated the stands across all harvesting alternatives by >45% of the number of stems (N/ha) throughout. The number of stems for the middle and lower canopy (species group 2) decreased in the “no harvesting” alternative compared to all other harvesting alternatives (Figure 3). The top canopy species (group 1) contributed over 47% of all ingrowths during the simulation, while species groups 2 and 3 contributed about 33% and 18%, respectively. The proportion of trees upgrowing initially increased and then decreased with greater Dbh. The highest proportion of up-growing trees was in species group 1. The highest proportion of stems growing into higher classes was attained in the “varying intensity” harvesting alternatives under species group 1 between Dbh class 20–25 cm, but at around 40–45 cm in the other alternatives as indicated in supplementary materials.

Figure 3. The number of stem densities across species groups and harvesting options.

The “no harvesting” alternative resulted in an overall lowest number of stems in Dbh class 5–25 cm but a relatively higher number of bigger diameter stems in species group 1 compared to other groups. There was a relatively small drop in the number of stems for smaller stems (Dbh 5–25 cm) for all species groups but an increased number of harvestable stems (Dbh ≥25–50 cm) for species group 1 for both “Varying” and “Uniform intensity” harvesting alternatives. Species groups 2 and 3, on the other hand, showed a relatively small increase in the number of stems with Dbh 20–30 cm for both “Varying” and “Uniform intensity” harvesting alternatives during the simulation process. Generally, species group 1 had a higher basal area proportion (%) compared to other species groups at the end of the study (Figure 4).

Figure 4. The change in the number of stems across species groups.

Stand-level growth trajectory

In the “no harvesting” alternative, the stand basal area increased from 8.9 m²ha⁻¹ to 18.7 m²ha⁻¹, where it tended to level off. In the “uniform intensity” and “varying harvesting intensity” the basal area increased less, to 13.8 m² ha⁻¹ and 13.9 m²ha⁻¹, respectively (Figure 5). The two harvesting alternatives had a sharper rise in basal area recovery following the first harvesting cycle compared to the second cycle.

Figure 5. Basal area growth across species groups.

The observed initial stand volume was 87.5 m³ha⁻¹. At the end of the simulation period; the “no harvesting” alternative had the highest stand volume of 175 m³ha⁻¹ with a net growth of 83 m³ha⁻¹ and annual growth of 0.8 m³ha⁻¹yr⁻¹. The “uniform intensity” and “varying intensity”, however, attained a relatively similar stand volume of 130 m³ ha⁻¹ and 134 m³ ha⁻¹ with a net growth of 9.7 and 5.5 m³ha⁻¹yr⁻¹ which included the harvested volume. The total harvested volume for “uniform intensity” was 33 m³ha⁻¹ and that for “varying intensity was 41 m³ha⁻¹. In all harvesting alternatives, the stand volume was dominated by species group 1 with a proportion of 71–73%.

The initial above-ground biomass carbon (AGB) for the stand was 24.1tC ha⁻¹. The tree harvest and mortality resulted in the loss of AGB carbon. Nevertheless, the “No harvesting” alternative had the highest net AGB carbon, valuing 42.6tC as compared to the “intensity (39.6tC) and “Uniform intensity” (33.2tC), respectively, at the end of the simulation period (Table 2). On the other hand, the stand basal area growth was assumed to have a direct linear relationship with in-growth (from saplings into trees with Dbh ≥5 cm) and mortality (Figure 6).

Figure 6. Relationship between stand basal area growth and mortality.

In the “no harvesting” alternative, the stand basal area increased from 8.9 m²ha⁻¹ to 18.7 m²ha⁻¹, tending to level off. In the “uniform intensity” and “varying intensity”, however, the basal area increased less, to 13.8 m² ha⁻¹ and 13.9 m²ha⁻¹, respectively. There was a sharper rise in basal area recovery following the first harvesting cycle compared to the second cycle in the two harvesting alternatives.

Discussion

Reliability of the growth model

The formulated diameter growth model was the basis for the development of the simulation system for Kitulangalo miombo stands. Formulating an individual species growth model could be necessary to facilitate diameter growth modelling and thus simulation system; however, the vast nature of miombo species necessitated the species grouping system using canopy structure and growth characteristics of miombo species (Njoghomi et al. 2021), which, however, seemed to improve the model fit. Several other models, such as for Julbernnadia globiflora, and Combretum molle were formulated during the process but were dropped due to unexpected biological behaviour resulting from interactions among the independent variables. Because the aim was to get a relatively simple model with robust diameter growth to facilitate stand growth simulation, the selected, simple model with Dbh and basal area (G) which introduced competition seemed statistically fit. The applied harvesting alternatives across the species groups also regulated the stand density, and basal area, which consequently controlled stand growth. Harvesting a few bigger size tree species (for groups 1 and 2) had a higher impact on the stand basal area, and thus stand growth. Generally, the model performed well statistically. Despite low regression coefficient of the model, its alignment to the biological behaviours of miombo growth accredited the selected model. We considered this feature as one of the strengths of the developed model which depicted the actual patterns of miombo growth and development (Mugasha et al. 2017). On the other hand, the limited number of observed mortality and recruits limited the formulation of mortality and recruitment sub-models, which could have disavailed our simulation matrix system. We consider these as crucial areas to improve simulation studies on stand growth and management in miombo woodlands in the future. In addition, we considered the observed steady-state conditions and quick stabilization in stand growth and development some years after second harvesting, as the effect of improved growth when the competition is reduced. This phenomenon is also parallel to the effect of the application of the density-dependent in-growth, and mortality. Similar to the findings of Valkonen et al. (2008), Mbwambo et al. (2008), and (Njoghomi et al. 2021), suggesting that the stand mortality rate must be lower than the recruitment rate for the sustainability of the miombo woodlands. For research purposes, our experiment did not involve fire (as it was excluded in the setting of the experiment); however, frequent annual fires would reduce the fuel load on the forest floor, consequently affecting regeneration, and thus recruitment.

The increment model was, however, applied only to the data it was constructed with which had a small range (5–60 cm). We, therefore, do not know how it will behave in other miombo with a wide range of diameters and other climatic conditions.

Stand dynamics without management

Without harvesting, the growth model parameter (${\beta }_2)$ continued to predict higher diameter increments, particularly to individuals in the top canopy species group (group 1). This increased diameter growth results in increased basal area growth for this group, ultimately inducing strong growth competition against the middle and lower canopy species. Further, the resulting basal area growth from species group 1, therefore, accelerates tree mortality for other species groups, while suppressing an overall stand ingrowth (Szwagrzyk and Szewczyk 2001). The dominance of species group 1 will therefore increase, thus compromising the future species composition, stand structure and diversity (Mangora 2005; Kalaba et al. 2012). Sustaining stands at basal area threshold levels shown by (Njoghomi et al. 2021) are appropriate for sustainable growth and the perpetual supply of ecosystem services and biodiversity functions in miombo woodlands. The “No harvesting” alternative implicitly explains the complete enclosed management system in miombo woodlands which is more common in many state-owned miombo forests in Tanzania, where harvesting is restricted (Luoga et al. 2005; Kajembe et al. 2015). This simulation result is important for discussion on sustainable management of miombo woodlands, thus threatening the future stand structure and composition.

Stand dynamics with management

Generally, harvesting trees create canopy gaps that essentially reduce growth competition among stand trees; thus, higher diameter increments (predicted by the model parameter ${\beta }_2)$, more ingrowth, and therefore, more stems up-growing into higher diameter classes for species group 1 (Acevedo et al. 1995; Valkonen et al. 2008). This phenomenon was demonstrated by harvesting treatments with wood removals from the large diameter classes. This was followed by higher diameter increments from the remaining stand trees, and thus, quick recovery in basal area, volume and other parameters’ growth for species group 1(the top canopy trees) after every cutting cycle. In addition, the created canopy gaps through harvesting promoted higher number of small stems recruited into the smallest diameter class (5–10 cm), thus adding to the overall stand basal area growth. Similar findings were reported by Muvengwi et al. (2020) in Zimbabwe and Syampungani (2009) in Zambia whereby extraction of fuel wood induced higher recruitments, and thus, overall stand growth. This situation, however, confirmed a general notion of a reversed J-shaped diameter distribution demonstrated across all species groups, throughout the simulation period in the treatment alternatives (Isango et al. 2007; Zheng and Zhou 2010). We consider the two harvesting options as examples of treatments to think of for sustainable forest development in Miombo woodlands. The extent of wood removal in the two harvesting alternatives and the resulting stand basal area recovery at KFR demonstrated the capacity of miombo woodlands to immediately recover from the slight to medium disturbances (Ribeiro et al. 2015b; Gonçalves et al. 2017) under sustainable management. The sustainability of the stand structure, especially the two-dimensional diameter distribution-species composition matrix will have to be involved in the consideration.

Suggestions on stand management from aspects of stand-level growth

Greater estimates of stand basal area, volume, above-ground biomass (AGB) and carbon under the “no harvesting” alternative point toward the presence of bigger diameter trees which were subject to harvesting under other alternatives. While the “no harvesting” management option promotes the accumulation of the AGB and carbon storage in large trees (Bertsch-Hoermann et al. 2021; Mikolas et al. 2021), the simulation results through harvesting treatments diminished AGB storage through harvesting and mortality, but subsequently increased carbon sequestration because of the increased growth rate and the resulted number of ingrowths (Malimbwi et al. 1994; Mbwambo et al. 2008; Lyimo and Shaaban 2015; Njoghomi et al. 2021). The Government of Tanzania, conservationists and the communities are currently in dilemma whether to store AGB carbon in our forests and woodlands (which benefits biodiversity conservation, ecosystem functioning and less economic gains) or, to embark on selective timber harvesting (demonstrated in this study) resulting into enhanced growth, carbon sequestration and timber trades as a response towards the global climate change mitigation initiatives and the needs to conserve of biological diversity. Management options that can balance the AGB storage and carbon sequestration dynamics (Mikolas et al. 2021) are a new area of research that will enhance sustainable forest management, increase community benefits and climate change mitigations through the implementation of future REDD+ programs and other forest landscape restoration initiatives in the country such as AFR100 and the Paris Challenge 2016.

Future research works

At the initial steps of simulation, the model would give rather unbiased increment estimates by default, as more recurrent steps were added on, deviations from the initial conditions could have emerged in background factors, especially related to standing density in the unmanaged alternatives and the species composition within the aggregate species groups. It was, however, impossible to gather sufficient observations for modelling all the species (61 in this data) in such complex woodlands in the context of limited results and time. It would not be justifiable to conclude very much from the results on the ecological and social sustainability of the management of miombo woodlands. Stand density, diameter distribution, and species composition are just the first step in that direction. For instance, the tree species have different ecological niches (Menon et al. 2010; Baldeck et al. 2013), regeneration, survival, and competition strategies (Vieira and Scariot 2006), and a multitude of timber and nontimber uses (Chirwa et al. 2008).

To guarantee the sustainability of miombo stands constitutes a hugely complex undertaking requiring major steps in research and the accumulation of knowledge (Valkonen et al. 2008). Further, longer observation periods on extended networks of permanent experimental plots and the involvement of other potential parameters such as fire, animal disturbance, and the influence of changing climates in miombo woodlands are crucial for a better understanding of stand dynamics and the resulting carbon sequestration dynamics in the structure and composition of miombo woodlands for sustainable planning and management and perpetual supply of ecosystem services the community.

Conclusions

The study showed that it is possible to build relevant simulation systems for typical miombo forests composed of rather basic empirical models based on datasets gained from relatively short observation periods on sample plots. The development of such a research approach and infrastructure for miombo is still at the very beginning, along with the research infrastructure as a whole. However, the results point towards a high degree of positive response to careful selection harvesting in lightly degraded, recovering miombo woodlands. The role of Kitulangalo Forest Reserve as a prime example of such a semi-intact forest and a key research platform is confirmed.

Acknowledgements

I acknowledge partial financial support from Luke’s INFORES-Tanzania project, funded by the Ministry of Foreign Affairs of Finland.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Supplementary material

Supplemental data for this article can be accessed online at https://doi.org/10.1080/13416979.2023.2239136