Document Type: Research Paper
Authors
^{1} Research Institute of Forests and Rangelands Tehran
^{2} Science and Research Branch Tehran
^{3} University of Tehran
Abstract
Keywords
Geostatistically estimation and mapping of forest stock in a natural unmanaged forest in the Caspian region of Iran
R. Akhavan^{1*}, H. KiaDaliri^{2}, V. Etemad^{3}
1 Forest Research Division, Research Institute of Forests and Rangelands, Tehran, Iran
2 Faculty of natural Resources, Islamic Azad University, Science and Research Branch, Tehran, Iran
3 Dept. of Forestry, University of Tehran, Karaj, Iran
* Corresponding author’s Email: akhavan@rifrac.ir
(Received: April. 07.2014 Accepted: July. 30.2014)
ABSTRACT
Estimation and mapping of forest resources are preconditions for management, planning and research. In this study, we applied kriging interpolation of geostatistics for estimation and mapping of forest stock attributes in a natural, unevenaged, unmanaged forest in the Caspian region of northern Iran. The site of the study has an area of 516 ha and an elevation that ranges from 1100 to 1450 m a.s.l. Field sampling was performed on a 75m × 200m systematic grid using 309 georeferenced circular sample plots of 1000 m2 area. Experimental variograms were calculated and plotted for basal area (BA), volume (V) and stem density (N). Whereas the calculated variograms of BA and V exhibited spatial autocorrelation only after data stratification based on diameter size classes and tree species, the variogram of stem density displayed a moderate spatial structure that was fitted by a spherical model. Stem density was estimated by ordinary block kriging and the accuracy of estimation was validated by crossvalidation result. We conclude that geostatistical approaches have the potential to more accurately capture and describe the spatial variability of forest stock, and thus reduce the uncertainty in estimates of stem density as well as produce more accurate stem density maps of forests in comparison with the spatially uninformed classic method. Geostatistical methods provide a very suitable tool to derive more accurate estimates of growing stock, particularly in structurally complex, unmanaged, unevenaged forest such as this one from the Caspian region of northern Iran.
Keywords: Geostatistics, mapping, forest stock, unmanaged forest, Caspian region, Iran.
INTRODUCTION
The implementation of closetonature silviculture that aimed at conserving and restoring unique, important and ecologically and economically endangered (due to human disturbances and illegal exploitation) Caspian forests of Iran (Abdollahpour & Assadi Atui, 2005) that resemble natural forests, requires a thorough understanding of ecological processes that have created the structures and composition of these natural archetypes. Fortunately, natural mixed unevenaged hardwood deciduous forests that encompass all three major development stages (i.e., the decay, initial, and optimal stages; Leibundgut, 1993, Korpel, 1995) still exist along the southern part of the Caspian Sea and the northfacing aspects of Elborz Mountain in northern Iran (SaghebTalebi et al., 2004). These Caspian forests are characterized by a high plant species diversity of more than 80 woody species (SaghebTalebi et al., 2004) and afford a unique opportunity to study the spatial structures found in intact natural oldgrowth forests. Although proper maps of spatial patterns of trees are a prerequisite for forest planning and management, planners rarely have walltowall coverage of spatially explicit tree stem maps that depict the spatial distribution of forest stock attributes (e.g., stem density, basal area, and volume) over large areas. Instead, forest stock attributes are typically inventoried using a sampling process that lacks spatial resolution of stock attributes between sample points. Further, nonspatial independence suggests that many statistical tools and inferences may not be appropriate and a spatial extrapolation of results to any given point in the whole sample area may not be warranted. Thus, an unbiased interpolation method has to be employed for mapping forest stock attributes after sampling that permits forest managers to more closely delineate the major development stages based on natural spatial structures of tree species and enable researchers to more realistically model stand structures and dynamics. Geostatistics, which is concerned with the detection, modeling and estimation of spatial dependence of continuously distributed or regionalized variables (Isaak & Srivastava, 1989, Goovaerts, 1997), is a useful tool to describe the spatial structure and provide georeferenced estimates of forest attributes. Although the spatial distribution of trees in a particular stand represents a point pattern of discrete objects (Dale, 2000), stand attributes such as basal area, volume stock and density can be thought to be directly influenced by different spatially continuous variables such as solar radiation, soil characteristics and water nutrient availability, thus allowing to be considered as spatially continuous (Kint et al., 2003). Geostatistics, thus, provides a natural framework for estimation techniques in forest inventory sampling (Mandallaz, 1991) and has frequently been used to estimate and map forest resources based on forestscale surveys (Samra et al., 1989, Biondi et al., 1994, Gunnarsson et al., 1998, Tuominen et al., 2003, Montes et al., 2005, Freeman & Moisen, 2007, Tröltzsch et al., 2009, Brus et al., 2012), enabling the production of georeferenced data maps for basal area, density, standing volume and tree species groups at scales where these variables often reveal spatial autocorrelation. Following two studies in which geostatistics was applied to estimate and map forest stock attributes in a natural managed forest(Akhavan et al., 2010) and a plantation forest (Akhavan & KiaDaliri, 2010) in the Caspian region of Iran, this study aimed at extending this approach to a natural, unmanaged forest in this region of Iran. This study specifically aimed at (1) investigating the spatial variability and spatial structure of forest stock attributes (stem density, basal area, and volume) using variogram analysis in geostatistics, (2) estimating and mapping forest stock attributes using the geostatistical kriging interpolation method, and (3) examining whether kriging improves the estimation accuracy of these forest attributes in comparison with a spatially uninformed classic approach.
MATERIALS AND METHODS
Study Area
Caspian forests along the southern part of the Caspian Sea, span the elevational gradient of northfacing aspects of Elborz Mountain in northern Iran, occur in five main vegetation types (i.e., QuercoBuxetum, QuercoCarpinetum, ParrotioCarpinetum, Fagetum hyrcanum and Carpinetum orientale) and are remnant of the Tertiary era that extend for 800 km in an east to west direction (SaghebTalebi et al., 2004). The study area (516 ha) is located in the northern part of the fourth district (Chelir) of the educational and research forest station of Tehran University (Kheiroud) at 51˚ 40' E longitude and 36˚ 32' N latitude (Fig. 1). Elevation varies from 1100 m to 1450 m above sea level and slopes range from 5% to 85%. The mean annual temperature, precipitation and relative humidity are 15.3 ˚C, 1458 mm, and 83%, respectively. The climate is cold and wet in the winter and temperate in the summer without any dry season. The growing season is 270 days per year. The forest under investigation is a typical naturally seeded, unevenaged, mixed hardwood oldgrowth stand in the Fagetum hyrcanum association that experienced no silvicultural activity and has no known management history. The inventoried forest area is a mixture of broadleaf deciduous tree species. The dominant tree species are beech (Fagus orientalis Lipsky) and hornbeam (Carpinus betulus L.), alongside with maple (Acer velutinum Boiss. and A. cappadocicum Gled.), alder (Alnus subcordata C.A.M.)and oak (Quercus castaneifolia C.A.M.). Large topographic variation contributes to the heterogeneous nature of the stand (Fig. 1).
Fig. 1. Study area and sampling grid (75 m × 200 m).
Data collection
A network with 75 m (NS) × 200 m (WE) systematic rectangular grid was used for sapling. At each grid point, we established a circular sample plot of 1000 m2 surface area where we recorded the UTM coordinates of the plot center (Fig.1) and the diameter at breast height (d.b.h. in 1.3 m above the ground) of each tree with a d.b.h. that exceeded 7.5 cm. Further, the height of the closest tree to its respective plot center point and the largest tree in each plot were measured for volume estimation. The inventory was accomplished in the summer of 2011. Standlevel attributes of interest, i.e., basal area (BA), volume stock (V) and stem density (N) were computed for all trees as well as for trees in different size classes (i.e., small size [d.b.h. ≤ 32.5 cm, S], medium size [32.5< d.b.h. ≤ 52.5, M], large size [52.5 < d.b.h. ≤72.5 , L], and extra large size [d.b.h. > 72.5 cm, EL]; after SaghebTalebi et al., 2005), and separately for the two most dominant tree species (i.e., beech and hornbeam) that accounted for about 75% of the total species frequency in the study area.
Geostatistical approach for spatial autocorrelation analyses
To identify and describe the spatial dependency (i.e., spatial autocorrelation) of BA, V and N of all trees and the diameter size classes of the two dominant tree species, we used variogram (semivariance) analysis as our geostatistical approach. Three parameters are commonly used to describe and model the behavior of variogram: range, sill and nugget effect. The range is the distance where the spatial correlation disappears and the variogram levels off, the sill corresponds to the height of the variogram after leveling off, and the nugget effect is represented by the intercept of variogram on the ordinate axis. The ratio of the nugget effect to the sill is known as the relative nugget effect. This is a measure of the percentage of the variability in the data from sources other than spatial autocorrelation. A low relative nugget effect (≤ 25%) is a sign of strong spatial autocorrelation where the application of geostatistical techniques is particularly beneficial. Nugget effects between 25% and 75% indicate moderate spatial autocorrelation and high relative nugget effect (≥75%) can indicate either weak spatial autocorrelation in the population or spatial patterns at scales smaller than the sampling distance (Cambardella et al., 1994, Ganawa et al., 2003, Freeman & Moisen, 2007). In the current study, all the experimental variograms obtained were modeled using either a pure nugget effect or an additional spherical model. The selection criterion used was the minimal residual sum of squares. The spherical model is given by (Webster & Oliver, 2000):
2000): for
for (1)
Where h, c0, c, and a represent a particular lag vector, nugget effect, structural variance, and range, respectively. After normalization of the data, only omnidirectional (isotropic) variograms were calculated, because experimental variogram surfaces showed no variogram anisotropy. By common convention, the analysis was restricted to distances of half of the study area dimension (i.e., 2200 m). We further used the Kriging interpolation method as our geostatistical approach for possible interpolating values among sample plots and mapping the distribution of BA, V and N. Kriging computes surfaces of the best linear unbiased estimation of regionalized variables at unsampled points based on the spatial structure defined by the experimental semivariogram. Ordinary kriging (the most common type of kriging in practice, particularly in environmental sciences) of the regionalized variable at point is given by (Webster & Oliver, 2000):
(2)
Where, is the weight associated with the value of at the sampled point with the nonbiased condition:
(3) Kriging may be used for estimation at a single point (point kriging) or over an area (block kriging). In this study, since the mean is assumed to be stationary and unknown and because no largescale trend was observed, ordinary block kriging without trend was used. A 32 m × 32 m mesh (approximately the same area as a sample plot to emphasize the local variation around the sampling plots) was used to discretize the study area for block kriging interpolation. The estimates were obtained using the nearest 16 data plots within the maximum effective range of variograms that corresponded to the scale of autocorrelation. A crossvalidation (leaveoneout) was performed for selected datasets to evaluate the kriging results. Crossvalidation was evaluated by calculating the Mean Bias Error (MBE), which should ideally be equal zero, because kriging is unbiased (Webster & Oliver, 2000):
(4)
(5)
The accuracy of kriging was measured using the Root Mean Square Error (RMSE):
(6)
(7)
Where is the estimated value of regionalized variable at the location of ; N is
The number of sample plots and is the mean value of measured samples of interested attribute. The software package used for geostatistical analysis was GS+ version 9 (Gamma Design Software, LLC, Plain Well, MI).
Point pattern analyses
To further understand the spatial dependencies of BA, V and N, we computed several indices for quadrat counts (e.g., variancetomean ratio, Morisita’s index of dispersion, and Morisita’s standardized index of dispersion) to analyze separately the spatial point patterns for all trees, trees in different diameter size classes, and for beech and hornbeam. The variancetomean ratio, attributed to Fisher et al. (1922) is one of the oldest and simplest measures of dispersion. The ratio ( ) usually called the index of dispersion (I) (Bailey & Gatrell, 1995) and is based on the observation in a random pattern, described by the Poisson distribution, the variance equals the mean, so I = 1 for a random pattern. Ratios larger than 1 indicate clumping, while smaller ratios indicate regular or uniform pattern. The other frequently used quadratbased dispersion index is the Standardized Morisita index of dispersion. Smith
Gill (1975) set out to improve Morisita's index (Morisita, 1962) just described by putting it on an absolute scale from 1 to +1. To calculate this index, Morisita's index of dispersion (Id) was first calculated, along with two critical values of the uniform index (Mu) and the clumped index (Mc):
(8)
(9)
(10)
Where n is the sample size, x is the number of individuals, and are the values
of chisquared with (n1) degrees of freedom that have 2.5% or 97.5% of the area to the
right. Morisita's Id is 1 for a random distribution, >1 for a clumped distribution, and
a regular/uniform distribution. The Standardized Morisita index (Ip) is then calculated by one of the four following formulae:
When
(11)
When (12)
When (13)
When
(14)
The Ip ranges from 1.0 to +1.0, with 95% confidence limits at +0.5 and 0.5. Random patterns give an Ip of zero, clumped patterns give an index value of above zero and uniform patterns below zero.
RESULTS
Sampling
A total of 309 plots were sampled in the study area. Normalization tests showed that the data were not normally distributed, consequently, the data were transformed using the square root and log transforms to more closely approximate normal distributions. Table 1 shows the summary statistics of the sample plots. The distribution of tree sizes in the study area reveals an unevenaged stand diameter distribution with the majority of trees in the small and medium tree size classes (Table 2) and a few very large trees up to 300 cm d.b.h. (Fig. 2).
Geostatistical analyses
Because we did not find any variogram anisotropies, we fitted only omnidirectional variograms using a spherical model to which a nugget effects was added (Fig. 3, Table 3). No spatial autocorrelation was revealed in the experimental variograms for BA and V. Consequently, these attributes represented a pure nugget effect. However, N showed a moderate level of spatial autocorrelation (50%).
Table 1. Summary statistics of 309 sample plots
Attribute 
Mean 
Min 
Max 
SD 
CV % 
Skewness 
Kurtosis 
BA (m2 ha1) 
35.8 
8.6 
99.6 
11.5 
32.2 
1.34 
5.13 
V (m3 ha1) 
459.7 
95.1 
1452.4 
191.7 
41.7 
1.37 
4.15 
N (n ha1) 
279 
20 
1140 
166.2 
59.6 
1.61 
3.86 
SD: Standard Deviation; CV: Coefficient of Variation
Table 2. Frequency and proportion of trees in the four diameter classes in the study area
Relative frequency 
Absolute frequency 
Diameter class 
64.6% 
5561 
S 
17.4% 
1493 
M 
9.8% 
850 
L 
8.2% 
708 
EL 
100% 
8612 
Total 
S, Small size (d.b.h. ≤ 32.5 cm); M, Medium size (32.5
L, Large size (52.5 72.5 cm)
Table 3. Characteristics of fitted variograms for the three attributes of basal area (BA), volume (V), and stem density (N)
Attribute 
Lag (m) 
Fitted model 
Nugget effect 
Sill

Range (m) 
R2 % 
SpD % 
BA (m2 ha1) 
130 
Pure nugget effect 
0.86 
0.89 
 
7 
96.6 
V (m3 ha1) 
130 
Pure nugget effect 
0.11 
0.12 
 
37 
91.6 
N (n ha1) 
130 
Spherical 
0.11 
0.22 
533 
93 
50 
SpD (Spatial Dependence) = (SpD ≥75%, Weak; 25%< SpD <75%, Moderate; SpD ≤ 25%, Strong)
Fig. 2.Diameter distribution of trees in the study area.
Because BA and V did not show any spatial autocorrelation, kriging interpolation was used only for N (Table 4). A comparison of tables 1 and 4 reveals that the estimated mean stem densities are not very different. However,
compared to the spatially uninformed classic method, a variance reduction in the estimate of N of approximately 72% was achieved using kriging interpolation. Figure 4 shows a kriging map and an error map for stem density over the study area.
Fig. 3. Graphs of isotropic variograms and fitted models for the three attributes of basal area (BA), volume (V), and stem density (N) in the study area. Solid lines represent the models, filled circles represent pointestimates. The minimum number of pairs for each lag. distance is 274
Table 4. Summary statistics of kriging results for forest stem density (n ha1)
Mean 
Min 
Max 
SD 
CV% 
260 
67 
597.32 
82.5 
31.7% 
.
Evaluating the kriging results of N with the Mean Bias Error and Root Mean Square Error and their respective relative values (Table 5)
revealed a low relative Mean Bias Error of less than 10% (MBEr = 6.9%). A crossvalidation graph of stem density confirms the accuracy of our estimation (Fig. 5).
Variography after data classification based on tree size classes and species
Although BA and V did not show any spatial variability when considering all trees, the results were quite different when we computed the experimental variograms by tree size clas
ses and species, individually. When considering trees of different diameter size classes, the spatial variability of BA (Fig. 6) and V (Fig. 7) clearly weakens from the small diameter size class towards the extra large class, with an approximately moderate spatial dependence observed in the small and medium size classes and only a pure nugget effect in the large and extra large classes. Similarly, when examining the spatial variability of BA by tree species, both beech and hornbeam individually exhibited moderate spatial autocorrelation (Fig. 8) that was completely masked when examining the spatial variability of all species (Fig. 6)
Table 5. Validation results of kriging interpolation for stem density (n ha1)
MBE 
RMSE 
MBEr 
RMSEr 
19.4 
147.7 
6.9% 
52.9% 
MBE: Mean Bias Error; RMSE: Root Mean Square Error; MBEr: relative MBE; RMSEr: relative RMSE
Fig. 4. Kriging map (a) and error (standard deviation) map (b) of stem density in the study area
Fig. 5. Crossvalidation graph for stem density (N). The solid line represents the regression and the dotted line is the 45 degree line; the best unbiased estimation occurs when the solid line coincides with the dotted line
Point pattern analyses
Point pattern analyses revealed that the spatial distributions of all trees, the two dominant tree
species and all diameter size classes except the extra large class were clumped in the study area (Table 6).
Fig. 6. Graphs of isotropic variograms and fitted models for basal area (BA) based on the four diameter size classes in the study area. Solid lines represent the models, filled circles represent pointestimates. Minimum number of pairs for each lag distance is 258. Note: S= Small size, M= Medium size, L= Large size, and EL= Extra Large size.
Fig. 7. Graphs of isotropic variograms and fitted models for volume (V) based on the four diameter size classes in the study area. The abbreviations and explanations are the same as for Fig. 6.
Table 6. Spatial patterns for all trees, size classes and the two dominant tree species in the study area

Index 



Id 
Ip 
Spatial patterns 

All trees 
9.91 
1.32 
0.500 
Clumped* 
S 
15.91 
1.83 
0.501 
Clumped* 
M 
2.6 
1.33 
0.500 
Clumped* 
L 
1.58 
1.21 
0.500 
Clumped* 
EL 
0.958 
0.98 
0.139 
Uniform* 
Beech 
6.47 
2.01 
0.501 
Clumped* 
Hornbeam 
14.86 
1.87 
0.501 
Clumped* 
variancetomean ratio; Id, the Morisita’s index; Ip, the Standardized Morisita’s index; *, significant at p
Fig. 8. Graphs of isotropic variograms and fitted models for basal area (BA) for the two most abundant tree species in the study area. Solid lines represent the models, filled circles represent pointestimates. Minimum number of pairs for each lag distance is 274.
,
DISCUSSION
The spatial structure of the three forest attributes (i.e., basal area, volume and stem density), expressed by their spatial autocorrelation, differed in this unmanaged natural unevenaged deciduous forest. Whereas stem density behaved as a regionalized variable and exhibited a moderate spatial structure (Table 3), overall stand basal area and volume did not show any spatial autocorrelation and exhibited a pure nugget effect (Fig. 3). The reason for these apparent differences can be found in the spatial distribution of trees of different diameter size classes. While smaller trees exhibited a clumped spatial distribution, larger trees increasingly tended toward a random and regular spatial distribution. This development toward spatial randomness or regularity with increasing tree size has been demonstrated in many forest types (e.g., Szwagrzyk & Czerwczak, 1993, Zenner & Peck, 2009) and is the main reason why the kriging interpolation method embedded in a geostatistical approach was able to estimate stem density more accurately to reduce the variance for stem density estimation by approximately 70%, and to produce a smaller coefficient of variation with acceptable estimation accuracy (MBEr ≈ 7%; Table 5) in comparison with the spatially uninformed classic approach (Table 1). Because smaller trees were more numerous than larger trees in this oldgrowth forest (Fig. 2) and each individual tree contributes equally to stem density, the moderate spatial structure of the stem density of the stand is largely driven by the
clumped spatial pattern of the smaller trees (Table 6). This is in stark contrast to the observed lack of spatial structure of the stand basal area and volume (Fig. 3) that are principally driven by the spatial randomness of the larger trees that, despite their low density, account for a large proportion of the growing stock. Thus, the spatial distribution of larger trees has more influence on the spatial structure of the growing stock (basal area and volume) in forests than that of small and mediumsized trees that typically exhibit a moderate spatial structure for basal area, and volume (Figs. 6 and 7). This lack of spatial structure of basal area and volume appears to be independent of changes induced by management, as was shown for a stand in the nearby Namkhane district (Fig. 1) that was subject to a management plan that prescribed a closetonature silviculture (Akhavan et al., 2010) and resulted in a diameter distribution that was similar to that of the unmanaged Chelir forest investigated in this study. Because larger trees were retained in the Namkhane forest following two harvest entries (Fig. 9), no spatial structures were detected for stand basal area and volume in Namkhane as well. The spatial structure of stem density differed between the unmanaged and managed stands, however, and was most likely influenced by management interventions such as road construction, plantation, and harvesting that may have masked a clear signal of the spatial structure of stem density in the Namkhane district (Akhavan et al., 2010).
Fig. 9. Diameter distribution of trees in the Namkhane district; after Akhavan et al. (2010).
The spatial structures of stem density and basal area in the current study are quite different from those reported by Akhavan and KiaDaliri (2010), who applied the same geostatistical approach in an eighteen yearold maple (Acer velutinum Boiss.) plantation in the Caspian region of Iran. In that study, stem density did not exhibit any spatial structure whereas basal area did. This was likely due to the regular spacing of the trees and the homogeneity of tree sizes in the plantation that was still in the early stage of stand development. Neither windstorms and droughts nor intraspecific competition had been strong enough to remove stems and break up the initial regular spatial distributions (initial planting space of 3 m × 3 m) after 18 years. Consequently, the regularity of the spatial distributions of the plants remained intact and no spatial autocorrelation was found for stem density. In contrast, the homogeneity of tree sizes in the maple plantation had a low coefficient of variation for d.b.h (35%), which was less than half of the coefficient of variation for d.b.h. observed in this study (79%), was sufficient to induce spatial autocorrelation for stand basal area. However, it is expected that as stem numbers decline until harvesting time (at approximately 80 years), stem density will exhibit spatial autocorrelation and the spatial structure of basal area will disappear as the homogeneity of tree sizes decreases. The low, medium and high density areas that are typically associated with mature, middleaged, and young stands, respectively, have become clearly visible in the kriged stem density maps (Fig. 4a). These kriged stem density maps, thus, permit an indirect estimation of forest stock for any point in the area. For example, the point I on the digital map of figure 4a has a stem density of 170 n ha1. Based on the diameter size classification in table 2, there are 64.6% (110 trees), 17.4% (30 trees), 9.8% (16 trees), and 8.2% (14 trees) of trees in the small, medium, large and extra large size classes at the point I, respectively. Therefore, if we use the midpoint value of each size class (namely, 20, 42.5 and 62.5 cm for the first three size classes, respectively; for the largest size class this depends on the maximum diameter size in the studied area), we can indirectly estimate the growing stock at point I. To obtain a more precise estimate, more and narrower size classes in the upper range of the size distribution would of course be beneficial due to their relatively larger impact on forest volume stock. Nonetheless, the advantage of a geostatistical approach that enables the precise quantification of the statistical error in a kriged error map (Fig. 4b) is that it identifies areas with higher errors that can then be covered with extra sample plots to reduce this estimation error. Each tree species has a specific physiological age and longevity that is reflected in a specific autocorrelation diagram (i.e., variogram) that may be masked when investigating the overall stand structure. While we failed to detect any spatial autocorrelation for basal area when analyzing all tree species that were present at the same time (Fig. 3), we did detect a moderate spatial structure for basal area for the two dominant species (i.e., beech and hornbeam with clumped spatial distribution; Table 6), individually (Fig. 8). Hence, it appears that the spatial structures of growing stock attributes in mixed stands, which are often unevenaged in this region, may be weaker than those observed in pure and evenaged stands, particularly if not analyzed separately by species. For this reason, Akhavan and Kia Daliri (2010) may have been able to detect a spatial structure for basal area in their plantation forest. In general, we hypothesize that if the diameter distribution observed in a natural forest has a narrower range (e.g., no larger trees present) and has no missing size classes (i.e., is without any interruption), the growing stock appears to have a stronger spatial structure with a higher autocorrelation, but this tentative hypothesis needs to be put to more rigorous testing. Geostatistical approaches were instrumental for the identification of the spatial structures of stem density in this natural unmanaged forest. In closetonature forestry, which is increasingly being applied throughout the world, kriged stem density maps can provide greater assistance for forest management planning and for the identification of forest development stages (Akhavan et al., 2012) than traditional tree size class maps. Whereas the identification of development stages is straightforward and fairly easy in monospecific and evenaged forests, the detection of these stages is often difficult and imprecise in mixed and unevenaged forests due to different habitat requirements of tree species. Here, kriged stem density maps could make an effective contribution toward a more accurate identification and delineation of these stages. Spatial interpolation of forest tree attributes using geostatistical approaches (kriging) can provide useful outputs at the regional scale as well, facilitating rapid spatial analyses that enable the user to test the influence of different factors on forest stock attributes (Merganič et al., 2004). These maps can be very useful for forest owners and managers as a suitable tool to gain more insight into the spatial distribution of the growing stock. This information is also beneficial as a guide map for writing forest management plans, i.e., to identify harvesting areas and to plan road networks, interplanting, and silvicultural interventions (i.e., thinning and light thinning) that are based on forest density and stock distribution, to ensure that forestry is done more sustainably.
ACKNOWLEDGMENTS
This research was financially supported by Iran National Science Foundation (INSF) with grant No. 88001186 and partially supported by the Research Institute of Forests and Rangelands (RIFR), Iran. The authors thank M. Hasani for assistance with fieldwork and data collection and Kh. Mirakhorlou for assistance with the preparation of maps.