Geostatistically estimation and mapping of forest stock in a natural unmanaged forest in the Caspian region of Iran

Geostatistically estimation and mapping of forest stock in a natural unmanaged forest in the Caspian region of Iran

 R. Akhavan^1*, H. Kia-Daliri², V. Etemad³

 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 E-mail: akhavan@rifr-ac.ir

(Received: April. 07.2014   Accepted: July. 30.2014)

INTRODUCTION

The implementation of close-to-nature silvi-culture that aimed at conserving and restor-ing unique, important and ecologically and economically endangered (due to human dis-turbances and illegal exploitation)  Caspian forests of Iran (Abdollahpour & Assadi Atui, 2005) that resemble natural forests, requires a thorough understanding of ecological pro-cesses that have created the structures and composition of these natural archetypes. For-tunately, natural mixed uneven-aged hard-wood deciduous forests that encompass all three major development stages (i.e., the de-cay, initial, and optimal stages; Leibundgut, 1993, Korpel, 1995) still exist along the southern part of the Caspian Sea and the north-facing aspects of Elborz Mountain in northern Iran (Sagheb-Talebi et al., 2004). These Caspian forests are characterized by a high plant species diversity of more than 80 woody species (Sagheb-Talebi et al., 2004) and afford a unique opportunity to study the spatial structures found in intact natural old-growth forests. Although proper maps of spatial patterns of trees are a prerequisite for forest planning and management, planners rarely have wall-to-wall 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 at-tributes between sample points. Further, non-spatial 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 un-biased 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 pro-vide geo-referenced estimates of forest attributes. Although the spatial distribution of trees in a particular stand represents a point pat-tern of discrete objects (Dale, 2000), stand at-tributes such as basal area, volume stock and density can be thought to be directly influenced by different spatially continuous varia-bles such as solar radiation, soil characteris-tics and water nutrient availability, thus al-lowing to be considered as spatially continuous (Kint et al., 2003). Geostatistics, thus, pro-vides 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 forest-scale 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 geo-referenced data maps for basal area, density, standing volume and tree species groups at scales where these variables often reveal spatial auto-correlation. Following two studies in which geostatistics was applied to estimate and map forest stock at-tributes in a natural managed forest(Akhavan et al., 2010) and a plantation forest (Akhavan & Kia-Daliri, 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 north-facing aspects of Elborz Mountain in northern Iran, occur in five main vegetation types (i.e., Querco-Buxetum, Querco-Carpinetum, Parrotio-Carpinetum, Fagetum hyrcanum and Carpinetum orientale) and are remnant of the Tertiary era that extend for 800 km in an east to west direction (Sagheb-Talebi 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 lon-gitude and 36˚ 32' N latitude (Fig. 1). Eleva-tion 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 investi-gation is a typical naturally seeded, uneven-aged, mixed hardwood old-growth stand in the Fagetum hyrcanum association that experi-enced no silvicultural activity and has no known management history. The inventoried forest area is a mixture of broad-leaf decidu-ous tree species. The dominant tree species are beech (Fagus orientalis Lipsky) and horn-beam (Carpinus betulus L.), alongside with maple (Acer velutinum Boiss. and A. cappadoci-cum 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 (N-S) × 200 m (W-E) 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. Stand-level 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 Sagheb-Talebi 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 auto-correlation analyses

To identify and describe the spatial depend-ency (i.e., spatial auto-correlation) of BA, V and N of all trees and the diameter size clas-ses of the two dominant tree species, we used variogram (semi-variance) 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 auto-correlation. A low relative nugget effect (≤ 25%) is a sign of strong spatial auto-correlation where the application of geostatistical techniques is particularly beneficial. Nugget effects between 25% and 75% indicate moderate spatial auto-correlation 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 experi-mental variogram surfaces showed no vario-gram 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 un-sampled points based on the spatial structure defined by the experimental semi-variogram. 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 as-sumed to be stationary and unknown and because no large-scale 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 ob-tained using the nearest 16 data plots within the maximum effective range of variograms that corresponded to the scale of auto-correlation. A cross-validation (leave-one-out) was performed for selected datasets to evaluate the kriging results. Cross-validation was evaluated by calculating the Mean Bias Error (MBE), which should ideally be equal zero, because kriging is unbiased (Webster & Oli-ver, 2000):

 (4)

                                   (5)

The accuracy of kriging was measured using the Root Mean Square Error (RMSE):

              (6)

                              (7)

 Where  is the estimated value of re-gionalized variable at the location of ; N is

The number of sample plots and  is the mean value of measured samples of interest-ed 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 dependen-cies of BA, V and N, we computed several indices for quadrat counts (e.g., variance-to-mean ratio, Morisita’s index of dispersion, and Morisita’s standardized index of disper-sion) to analyze separately the spatial point patterns for all trees, trees in different diame-ter size classes, and for beech and hornbeam. The variance-to-mean ratio, attributed to Fisher et al. (1922) is one of the oldest and simplest measures of dispersion. The ratio        ( ) usually called the index of disper-sion (I) (Bailey & Gatrell, 1995) and is based on the observation in a random pattern, de-scribed by the Poisson distribution, the vari-ance equals the mean, so I = 1 for a random pattern. Ratios larger than 1 indicate clump-ing, while smaller ratios indicate regular or uniform pattern. The other frequently used quadrat-based dispersion index is the Stand-ardized 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 chi-squared with (n-1) 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% con-fidence limits at +0.5 and -0.5. Random pat-terns 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, conse-quently, 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 uneven-aged 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 aniso-tropies, we fitted only omni-directional vario-grams using a spherical model to which a nugget effects was added (Fig. 3, Table 3). No spatial auto-correlation was revealed in the experimental variograms for BA and V. Con-sequently, these attributes represented a pure nugget effect. However, N showed a moderate level of spatial auto-correlation (50%).

Table 1. Summary statistics of 309 sample plots

SD: Standard Deviation; CV: Coefficient of Variation

Table 2. Frequency and proportion of trees in the four diameter classes in the study area

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)

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 auto-correlation, 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 point-estimates. The minimum number of pairs for each lag. distance is 274

Table 4. Summary statistics of kriging results for forest stem density (n ha-1)

.

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 cross-validation 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 re-sults were quite different when we computed the experimental variograms by tree size clas-

ses and species, individually. When consider-ing 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 horn-beam individually exhibited moderate spatial auto-correlation (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 ha-1)

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. Cross-validation 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 point-estimates. 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

variance-to-mean 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 point-estimates. Minimum number of pairs for each lag distance is 274.

,

DISCUSSION

The spatial structure of the three forest attrib-utes (i.e., basal area, volume and stem density), expressed by their spatial autocorrelation, differed in this unmanaged natural uneven-aged 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 auto-correlation 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 old-growth 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 medium-sized 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 close-to-nature 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 man-aged 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).