Abstract

Aim of study: To examine the relationships between net primary productivity (NPP) and leaf area index (LAI) and modeling these parameters with stand parameters such as stand median diameter (dg), dominant height (htop), number of trees (N), stand basal area (BA) and stand volume (V).

Area of study: Pure Anatolian black pine (Pinus nigra J.F. Arnold) stands in semi-arid climatic conditions in the Black Sea backward region of Türkiye.

Material and methods: In this study, the Carnegie-Ames-Stanford Approach model was used to calculate NPP; LAI, dg, h_top, N, BA and V were calculated in 180 sample plots. The relations between NPP and LAI with stand parameters were modeled using multiple regression analysis, support vector machines (SVM) and deep learning (DL) techniques. Relationships between NPP and LAI were investigated according to stand developmental stages and crown closure classes.

Main results: The highest level of relations was obtained in the stands containing the a-b developmental stages (r=0.84). The most successful model in modeling NPP with stand parameters was obtained by DL method (model R²=0.64, test R²=0.51). Although DL method had higher success in modeling LAI with stand parameters, SVM method was found to be more successful in terms of model-test fit, and modeling success in independent data set.

Research highlights: Grouping parameters affecting NPP and LAI increased the level of correlation between them. In modeling NPP and LAI in relation to stand parameters, machine learning algorithms performed better than linear approach. The overfitting problem can be eliminated substantially by including arguments such as early stopping, network reduction and regularization in the network structure.

Additional key words: machine learning algorithms; support vector machines; deep learning.

Abbreviations used: ANN (artificial neural networks); BA (basal area); CASA (Carnegie-Ames-Stanford Approach); dg (median diameter); DL (deep learning); FE (Forest Enterprise); FPU (Forest Planning Unit); htop (dominant height); LAI (leaf area index); MRA (multiple linear regression); N (number of trees); NPP (net primary productivity); RDF (Regional Directorate of Forestry); SVM (support vector machines); V (stand volume).

Citation: Bulut, S; Günlü, A; Keles, S (2023). Assessment of the interactions among net primary productivity, leaf area index and stand parameters in pure Anatolian black pine stands: A case study from Türkiye. Forest Systems, Volume 32, Issue 1, e003.
https://doi.org/10.5424/fs/2023321-19615

Received: 02 Jun 2022. Accepted: 14 Feb 2023.

Funding agencies/institutions	Project / Grant
Scientific Research Project Unit of Çankırı Karatekin University	OF061218D07

Competing interests: The authors have declared that no competing interests exist.

Correspondence Sinan Bulut: sbulut@karatekin.edu.tr

Introduction Top

In terrestrial ecosystems, forests have an important role in the global carbon cycle with the functions they undertake. The effect of forest ecosystems on this process is evaluated by the amount of carbon they store. The net amount of carbon added by forests is determined by net primary productivity (NPP). NPP is the difference between the total amount of carbon produced by plants and the amount they spend through autotrophic respiration. NPP researched on global and regional scales is valuable data for the interpretation of carbon emissions and climate change. NPP is also a very important parameter used in ecological-based models (Yan et al., 2021; Quinto-Mosquera et al., 2021; Wang & Yue, 2022). Interactions among variables in terrestrial ecosystems are affected by factors such as temperature and precipitation. NPP is a very successful indicator in evaluating the effects of climate change on vegetation in terrestrial ecosystems which is a complex system (Gower et al., 1999; Pan et al., 2015; Zhang et al., 2019).

The Carnegie-Ames-Stanford Approach (CASA) model, which is a process-based model, is frequently used for the estimation of NPP. This model has four main computational variables. These are: absorbed photosynthetically active radiation, light utilization efficiency, water, and temperature (Potter et al., 1993). The CASA model is remote sensing-based and flexible model. By enabling the use of remote sensing and geographic information systems technologies, databases can be created especially for large-scale areas, and calibration processes can be performed using data related to local areas. In addition, the explanatory variables used can be updated for different regions and different time periods, and calculations can be made in different variations. It allows calculation of NPP for global, national, and local scale areas in monthly, seasonal, and annual periods (Gülbeyaz, 2018).

The physiological processes controlling NPP are directly related to the leaf area index (LAI) and are frequently used in NPP models (Davidson, 2002). LAI is the amount of leaf area per unit area that is actively photosynthesizing. Two methods, direct and indirect, are used in the determination of LAI. The direct method is applied in the form of debris collection and determining the amount of leaves per unit area captured by leaf traps. This method is labor intensive and time consuming. The indirect method is the estimation of LAI using optical sensors (hemispherical photography). The indirect method is both easy to apply for large areas and advantageous in terms of time and cost (Chen et al., 1997; Jonckheere et al., 2004; Chen et al., 2022).

LAI, which is an important indicator of growth and productivity of forest ecosystems, is also an effective indicator of stand crown closure. Any change in the LAI as a result of natural disasters such as storms, snow and wind-induced overturns, drought and management interventions directly affect growth and productivity. The fact that a forest stand has a high LAI value indicates that the site productivity is good. Therefore, it is widely used in modeling stand parameters and shows high correlations with stand characteristics such as increment, growth, and volume. It is also a variable that affects photosynthesis, interception, evapotranspiration, carbon, nutrient, and energy balance. With long-term monitoring of LAI, changes in forest ecosystem productivity and climate dynamics can be better interpreted (Bonan, 1993; Sprintsin et al., 2011; Yan et al., 2012; Günlü et al., 2017; Bahrami et al., 2022).

NPP and LAI are widely used to assess ecosystem structures and dynamics. Luo et al. (2004) investigated climatic and soil controls on NPP and LAI in the Tibetan Plateau. They reported that the climatic conditions highly limited the NPP, and increased precipitation generally increased LAI. Chen et al. (2020) explored the spatiotemporal NPP and LAI of natural and plantation forests in China. They quantified relative contributions and assessed the interannual variabilities-trends of NPP and LAI. Pan et al. (2021) investigated the relationships between NPP and LAI on a global scale. They calculated the NPP/LAI ratio of global vegetation to determine the net primary production per unit leaf area in interpretation of interactions. In addition to these studies, there are also studies in which NPP is modeled with LAI or stand parameters and the relationships between them are investigated. Gholz (1982), Scurlock et al. (2001) and Bulut et al. (2019) modeled NPP using LAI and stand parameters for different study areas and tree species.

In this study, we focused on pure Anatolian black pine stands, which shows an important spreading area in the region. According to the Turkish National Forest Inventory, the Anatolian black pine stands cover an area of 4,199,623 hectares and 18.31% of the total forest areas in Türkiye (GDF, 2021). As a result of the ecological resilience of the Anatolian black pine, it is one of the most widely used tree species for reforestation worldwide (Cenni et al., 1998). For instance, it is widely used in afforestation especially in semi-arid climatic regions in Türkiye. Also, it is the most common and dominant tree species in the study area. Thus, we selected pure Anatolian black pine stands for modelling the relationships between NPP-LAI, NPP-stand parameters, and LAI-stand parameters. The main objectives of the study were to: i) predict the NPP with the CASA model, ii) calculate the stand parameters and LAI values from ground measurements in pure Anatolian black pine stands, and iii) explore and model the relationships between NPP-LAI, NPP-stand parameters, and LAI-stand parameters, using multiple regression analysis (MRA), support vector machines (SVM) and deep learning (DL) methods. We hypothesized that there would be a non-linear relationship between NPP and LAI, and the SVM and DL modeling techniques will increase the model success levels compared to MRA modeling technique.

Material and methods Top

Study area

This study was carried out in pure Anatolian black pine stands, which are distributed within the borders of Ankara Regional Directorate of Forestry (RDF) (Fig. 1). Ankara RDF consists of 11 Forest Enterprises (FE) and 44 Forest Planning Units (FPU). The study area is located between 313096-598185 east longitudes and 4414251-4544558 north latitudes (WGS 1984 UTM Zone 36N). The total area of Ankara RDF is 4,563,247.07 ha.

Figure 1. Location of the study area.

The mean elevation is 1071 m, and the mean slope is 17%. The annual mean temperature is 11.6 °C, and the mean annual total precipitation amount, 399.75 mm. The month with the highest mean monthly precipitation is May (70.85 mm), and the lowest month is July (20.07 mm). Anatolian black pine stands constitute the majority in the study area among the other tree species. Anatolian black pine is also frequently used in plantation areas in the region, since it is well adapted to the semi-arid climatic conditions. The stands under the dominance of Anatolian black pine in the study area are 213,868.56 ha, and pure Anatolian black pine stands cover an area of 136,985.96 ha. Other tree species in the study area are Pinus sylvestris, Pinus brutia, Abies nordmanniana, Fagus orientalis, Carpinus betulus, Quercus cerris, Quercus infectoria, Quercus frainetto, Quercus robur, Quercus petraea, Cedrus libani, Juniperus sp. and Populus tremula.

Field measurement

The field work was carried out in FEs where pure Anatolian black pine stands are distributed in Ankara RDF. The spatial distribution of sample plot numbers was arranged by proportioning the areas of pure Anatolian black pine covered in the relevant FPUs. Accordingly, a total of 180 sample plots were taken. Taking into account the stand crown closure in the field studies, sample plots were taken in the form of a circle, 400 m² in 3 (71-100%), 600 m² in 2 (41-70%) and 800 m² in 1 (11-40%) closed stands. Diameter values from the breast height level were measured bilaterally in trees with a diameter of less than 8 cm in young stands, and 8 cm and thicker in other stands, and the average values were taken. To determine the stand mean age, age was measured on 4-5 trees close to the stand mean diameter. Tree heights were measured on the dominant trees according to the method of 100 trees/ha to determine the site productivity. As a result of the inventory study, stand median diameter (dg), dominant height (h_top), number of trees (N), stand basal area (BA) and stand volume (V) values were calculated for each sample plot (Table 1). In the calculation of stand volume, the diameter and height-dependent double-entry volume equation (Eq. 1) developed for pure and natural Anatolian black pine stands distributed in Ankara RDF was used (Ercanlı, 2020):

v = 0.000097 \times d^{1.740878} \times h^{1.005521}

(1)

Table 1. Stand parameters, equations, and abbreviations.

Stand parameter	Equation	Abbreviation
d_g	${d̄}_{g} = \sqrt{\frac{{d̄}_{1}^{2} + {d̄}_{2}^{2} + ... + {d̄}_{2}^{n}}{n}}$	d_g= stand median diameter (cm)
d_g		d₁...d_n= measured tree diameter
h_top	$h_{top} = \frac{h_{1} + h_{2} + ... + h_{n}}{n}$	h_top= stand dominant height (m)
h_top	$h_{top} = \frac{h_{1} + h_{2} + ... + h_{n}}{n}$	h₁...h_n= measured tree height
N	$N = \frac{10000}{SSA} \times n$	N = number of trees/ha
N	$N = \frac{10000}{SSA} \times n$	SSA= size of sample plot
BA	$BA = \frac{10000}{SSA} \times (\frac{π}{4} x d^{2})$	BA= stand basal area (m²/ha)
BA	$BA = \frac{10000}{SSA} \times (\frac{π}{4} x d^{2})$	d= measured tree diameter
V	$V = \frac{10000}{SSA} \times \sum_{n = 1}^{n} v_{i}$	V= stand volume (m³/ha)
V	$V = \frac{10000}{SSA} \times \sum_{n = 1}^{n} v_{i}$	v_i= volume of tree

Calculation of net primary productivity

The NPP can be determined directly by making local measurements or can be predicted with statistical models or process-based models in forest ecosystems. In this study, CASA (Potter et al., 1993), which is a process-based model (Eqs. 2-9), was used for the estimation of NPP:

NPP = ε_{n} \times APAR

(2)

ε_{n} = ε_{\max} \times f (T) \times f (W)

(3)

f (T) = \frac{(T - T_{\min}) \times (T - T_{\max})}{(T - T_{\min}) \times (T - T_{\max}) - {(T - T_{opt})}^{2}}

(4)

f (W) = \frac{1 + LSWI}{1 + {LSWI}_{\max}}

(5)

LSWI = \frac{P_{nir} - P_{swir}}{P_{nir} + P_{swir}}

(6)

APAR = FPAR \times PAR

(7)

FPAR = \frac{(NDVI - {NDVI}_{\min}) \times 0.95}{({NDVI}_{\max} - {NDVI}_{\min})} + 0.05

(8)

PAR = R_{s} \times 0.50

(9)

where ε_n is light use efficiency; PAR is the photosynthetically active radiation; APAR is absorbed PAR; ε_max represents the maximum light use efficiency, and the ε_max values of young, middle-aged, and old stands used were 0.72, 0.57 and 0.52, respectively (Li & Zhou, 2015); f(T) is the temperature function; T is the atmospheric temperature; T_min, T_max and T_opt are minimum, maximum and optimal temperatures for photosynthetic activities, respectively; f(W) is the precipitation function; LSWI is the land surface water index; LSWImax is the maximum of LSWI values; the p_nir and p_swir variables represent the surface reflectance of the NIR (near-infrared) and SWIR (short-wavelength infrared) bands in Sentinel-2 satellite image, respectively; FPAR is the fraction of PAR absorbed by the vegetation cover; NDVI is the normalized difference vegetation index; NDVI_max and NDVI_min are the maximum and minimum of NDVI, respectively; and Rs denotes total solar radiation (MJ m^-2 day^-1). At the end of this process, maps with a spatial resolution of 10 m were produced for the NPP, which was calculated on a monthly (gC m^-2 month^-1) and annual (gC m^-2 year^-1) scale for the study area.

LAI values

A camera (Canon EOS 50 D) with a fisheye (Sigma 8 mm Fisheye lens) was used to determine the LAI value in each sample plot. Hemispherical photographs covering an area of 180o were taken. These photographs were taken from 5 different points to represent the sample plot, from a height of approximately 1.5 m, and towards the top of the stand. However, care was taken not to take photos when direct sunlight comes in. Then, the digitized photographs were analyzed using HemiView, H. version: 2.1 SR4 software program. By taking the average of the values obtained as a result of the analysis, the LAI values for each sample plot were calculated.

Modeling NPP and LAI with stand parameters

NPP was calculated using the CASA model, stand parameters such as LAI, dg, htop, N, BA, and V were obtained from the field inventory studies, and datasets were made ready for the modeling process. To test the created models, the data sets were divided into two groups: 70% were reserved for the model (126 items) and 30% for the test (54 items). MRA, SVM and DL modeling techniques were applied to the separated data sets. The modeling process using these techniques was carried out in three stages. In the first stage, normality tests of the NPP, LAI, dg, htop, N, BA, and V data sets were performed using Kolmogorov-Smirnov normality test (p > 0.05). In the second stage, the methods were applied using the model data set. In the third stage, the test data were run in the models created with the model data set.

Modeling techniques

The MRA method was performed to determine the linear relationships among NPP, LAI and stand parameters. IBM SPSS 23 software was used for the implementation of the multiple linear regression (MLR) method with stepwise procedure.

The SVM method is based on structural risk minimization approach and statistical learning theory. The purpose of this method is to obtain the optimal hyperplane that can separate the existing classes from each other. There are four different kernel functions such as linear, radial, polynomial, and sigmoid that can be used to create these limits (Zhao et al., 2019). The SVM method was applied in “e1071” package within the “R” software (R Development Core Team, 2014; Meyer et al., 2015).

The DL method is one of the machine learning algorithms and has a more complex structure than artificial neural networks (ANN). While ANN are implemented with one hidden layer, the DL method is implemented using more than one hidden layer. The flowchart of the DL method consists of the input layer, hidden layers, and output layer. Inputs are first transferred to the input layer. Each node provides an output value through its activation function. The outputs of the input layer are used as inputs in the next hidden layer. The neuron numbers in the input layer are the number of inputs in the processed data, and the neuron numbers in the output layer are equal to the number of outputs associated with inputs. The DL method was applied in “H2O” package within the “R” software (R Development Core Team, 2014; Aiello et al., 2015). Early stopping, network reduction and regularization arguments were used in “H2O” package. Early stopping argument was used so that the learning speed does not slow down. Network reduction argument was used to avoid noise learning and to reduce noise by eliminating irrelevant data. The regulation argument ensures that only the useful variable is added the model. On the other hand, it reduces the weights of the variables with less effect.

Model evaluation criteria

Different criteria were used to compare the observation data through the inventory study and the estimation results using modeling techniques. The performance measures used were correlation coefficient (r), coefficient of determination (R²), and root-mean squared error (RMSE).

r = \frac{\sum XY - (\sum X) (\sum Y) / n}{\sqrt{\sum (X^{2} - (\sum X^{2}) / n) (Y^{2} - (\sum Y^{2}) / n)}}

(10)

R^{2} = 1 - \frac{\sum_{i=1}^{n} {(Y_{g} - Y_{t})}^{2}}{\sum_{i=1}^{n} {(Y_{g} - Y_{og})}^{2}}

(11)

RMSE = \sqrt{\frac{1}{n} \sum_{i=1}^{n} {(h_{n})}^{2}}

(12)

where h_i, error; X, primary variable; Y, secondary variable; y_g, observation; y_t, estimation; y_og, mean observation; and n, number of samples. For high model success, r and R² values are expected to be high; RMSE value is expected to be small.

ResultsTop

Net primary productivity

NPP was calculated monthly and annually for 2019 using the CASA model (Fig. 2). The minimum annual NPP amount of pure Anatolian black pine stands in the sample plots was 284 gC m^-2 and the maximum amount was 960 gC m^-2. The annual average NPP amount of the stands was calculated as 602 gC m^-2.

Figure 2. Annual total net primary productivity (NPP) map of the pure Anatolian black pine stands.

Spatial distribution of annual total NPP amounts of Anatolian black pine stands was mapped (Fig. 2). Average NPP amounts of FEs were 664 gC m^-2 year^-1 in Nallıhan, 651 gC m^-2 year^-1 in Beypazarı, 623 gC m^-2 year^-1 in Çamlıdere, 646 gC m^-2 year^-1 in Kızılcahamam, 635 gC m^-2 year^-1 in Eskipazar, 562 gC m^-2 year^-1 in Çerkeş, 594 gC m^-2 year^-1 in Ilgaz, 621 gC m^-2 year^-1 in Çankırı and 630 gC m^-2 year^-1 in Ankara. The highest average NPP in the Regional Directorate was found in Nallıhan FE (664 gC m^-2 year^-1). It should be noted that these data only belong to pure Anatolian black pine stands. It would be inappropriate to make a more efficient/inefficient interpretation among FEs based on this information alone. There are many factors that can cause differences in NPP, and stand characteristics, topographic and edaphic factors in FEs are the most important ones. In addition, the timely and correct implementation of the silvicultural interventions on the stands can greatly contribute the productivity in the stands.

Leaf area index

Correlation levels of dg, h_top, N, BA, V and LAI for sample plots were -0.30, 0.11, 0.58, 0.35 and 0.36, respectively. LAI was grouped according to developmental stages and crown closure classes, and the relationships between NPP and these classes were examined. When classified according to their developmental stages, the highest average LAI value (2.729 m² m^-2) was found in the young stands (a developmental stages). The crown closure classes are determined according to the coverage ratio of the top roof to the ground surface, the increase in the amount of overlap is due to the increase in the amount of needles per unit area. Therefore, the highest average LAI value was determined in the 3rd (1.494 m² m^-2) and the lowest average value was determined in the 1st (0.654 m² m^-2) crown closure class.

Relationships among NPP, LAI, developmental stages, and crown closure classes

A parabolic relationship was determined between NPP and LAI (Fig. 3). The maximum value of NPP determined for pure Anatolian black pine stands was 960.038 gC m^-2 year^-1. NPP reached its maximum value at the breaking point where the LAI was 2.456 m² m^-2. A decreasing trend in NPP was determined when the LAI value exceeded 2.456 m² m^-2. While the relationship between NPP and LAI was positive until the breaking point (r = 0.42), it was found to be negative after the breaking point (r = 0.64).

Figure 3. Relationship between net primary productivity (NPP) and leaf area index (LAI) including graphic illustrating the maximum NPP value on LAI threshold value.

While examining the relationships between NPP and LAI, they were evaluated separately according to developmental stages and crown closure classes (Table 2). The LAI used as an independent variable was included in the analysis by applying transformations. For this purpose, the transformations of LAI², 1/LAI, 1/LAI², lnLAI, logLAI,

\sqrt{LAI}

and 1/

\sqrt{LAI}

were applied to the data sets of LAI. Since there were two stands belonging to the “a” developmental stage, they were combined with the “b” developmental stage in the analysis process. The highest mean LAI and NPP were found in “a” (young stands) developmental stage and 3-closed stands. The lowest NPP and LAI were at “d” developmental stage and 1-closed stands.

Table 2. Relationships between net primary productivity (NPP) and leaf area index (LAI) according to developmental stages and crown closure classes.

	Class	N	Mean LAI (m² m^-2)	Mean NPP (gC m^-2)	NPP model^[a]	R²	r
		180	1.254	601.610	Eq. 13	0.22	0.47
Developmental stage
a	0-7.9 cm	2	2.729	884.995	Eq. 14	0.70	0.84
b	8-19.9 cm	21	1.276	685.000
c	20-35.9 cm	88	1.202	617.822	Eq. 15	0.52	0.72
d	>36 cm	69	1.012	547.341	Eq. 16	0.20	0.45
Crown closure
1	11-40%	26	0.654	569.744	Eq. 17	0.38	0.62
2	41-70%	68	0.974	572.195	Eq. 18	0.21	0.45
3	71-100%	86	1.494	660.714	Eq. 19	0.37	0.61

NPP = 543.949 + (36.119 \times {LAI}^{2})

(13)

NPP = 1093.488 - (649.391 \times \frac{1}{LAI} + (183.788 \times \frac{1}{{LAI}^{2}})

(14)

NPP = 894.244 - (322.221 \times \frac{1}{LAI}) - (350.944 \times lnLAI) + (55.021 \times {LAI}^{2})

(15)

NPP = 456.589 + (91.165 \times LAI)

(16)

NPP = 521.380 + (128.738 \times {LAI}^{2})

(17)

NPP = 388.628 + (171.250 \times \frac{1}{{LAI}^{2}}) + (569.962 \times ln LAI)

(18)

NPP = 639.162 - (64.486 \times \frac{1}{{LAI}^{2}}) + (24.652 \times {LAI}^{2})

(19)

The model with the highest relationship between NPP and LAI was produced from stands belonging to “a” and “b” developmental stages (Fig. 4). In the scatterplots, the slope of the red line was expected to be close to 45 degrees, and the black points were expected to spread close to red line. When the scatterplots of NPP estimations vs CASA model NPP values were examined, it was seen that the best spread was obtained for “a-b” developmental stages (r = 0.84). Comparing the whole data set and the relationships that were divided into different groups, it was observed that certain classes gave higher correlations than the total data set. It was also observed that the relations between LAI and NPP were stronger in the groups devoted for “a-b” and “c” developmental stages and “1” and “3” crown closure stands.

Figure 4. Scatterplots of NPP (g C m^-2 year^-1) based on i) a-b developmental stages, ii) c developmental stage, iii) d developmental stage, iv) 1 crown closure, v) 2 crown closure, and vi) 3 crown closure.

Modeling of NPP and LAI with stand parameters

The terrestrial data sets obtained from the inventory study consist of the calculated stand parameters (dg, h_top, N, BA, and V). NPP and LAI were modeled using stand parameters. In addition, the NPP-LAI model was created to determine the level of modeling of NPP using only LAI. For modelling, 70% of the data sets were devoted for the model and 30% for testing. In the modeling phase, MRA, SVM and DL techniques were used.

As the first step, the MRA method was applied. Thus, the independent variables that were significant in the model to be created for the dependent variable were determined. Stand parameters included in the model created for NPP are h_top, N, BA, and V (Fig. 5, Eq. 20). The coefficient of determination of the resulting model was found to be 0.50, and when the test data were run in the same model, this value decreased to 0.30 (Fig. 5). A slight decrease in test success was expected, because the model was tested with a separate data set that was never encountered during the modeling phase. The test error of the NPP model was found to be 93.757 gC m^-2 year^-1. In the model created for LAI, dg, N and BA were entered as significant parameters (Fig. 5, Eq. 21). Although the model success was not high, there was no difference in success between the model and the test (model R² = 0.48, test R² = 0.47). This shows that the model can be applied to different data according to the level of success and that it has a stable structure. Although there was no difference between the model and test in the NPP model created with LAI (Fig. 5, Eq. 22), the model success was observed at a low level (R² = 0.20). The error value of the NPP-LAI model created was 99.807 gC m^-2 year^-1.

Figure 5. Observed and predicted values obtained from NPP-stand parameters, LAI-stand parameters and NPP-LAI models using the MRA method. Fitted model equations (20, 21, 22) are shown.

Significant independent variables were determined with the MRA models created. Relationships between dependent and independent variables were modeled using SVM and DL methods. In modeling NPP with stand parameters, the performance of the SVM method was found to be higher than the MRA method (Fig. 6). When the similarities between the predictions and observations of the models were compared, the similarity of the MRA method was 0.71, while this value was 0.79 for the SVM method. As a result of the test, SVM was more successful (r = 0.66). A similar situation was observed in terms of success in modeling LAI. While the MRA method modeled the LAI with an error of 0.383 m² m^-2 for the model data set, the error of the SVM method was 0.322 m² m^-2. In addition, no difference was observed between the model and test success of the SVM method (Model R² = 0.64, Test R² = 0.64). In the SVM model between NPP and LAI, a notable increase in success was not obtained compared to the MRA method.

Figure 6. Observed and predicted values obtained from NPP-stand parameters, LAI-stand parameters and NPP-LAI models using the SVM method.

With the independent variables taken from the MRA models, the DL method was applied for different neuron and layer numbers (Fig. 7). The layers were tried as 3, 7 and 10, and the number of neurons starting from 10 up to 100 in ten intervals. Thus, 30 model network trainings were conducted for each model. Among the models developed, the network with the closest model-test successes and the highest test success was selected and marked with black circles on the success graphs in Fig. 7. Findings of the selected best fit model and other methods were compared (Fig. 8). In the estimations made, it was expected that there will not be a high difference between the model and the test. However, due to the overfitting problem frequently encountered in the DL method, the model showed high success, but when the same model was tested with an independent data set, the success went down to very low levels. Model estimations trained with stand parameters of NPP achieved high success (r = 0.95). However, test success was low and not equivalent to model success (r = 0.54). As seen in the success graphs prepared, the difference between model and test success levels for NPP-stand parameters was quite large (model R² = 0.90, test R² = 0.29). For LAI-stand parameters, this difference was lower (model R² = 0.72, test R² = 0.46), while for NPP-LAI success results were similar (model R² = 0.24, test R² = 0.24).

Figure 7. Correlation changes of NPP-stand parameters, LAI-stand parameters and NPP-LAI models and tests according to the changing number of neurons and layers using the DL method.

Figure 8. Observed and predicted values obtained from NPP-stand parameters, LAI-stand parameters and NPP-LAI models using the DL method. N: number of neurons, L: number of layers.

In the applied DL method, overfitting problem was observed especially in the models of NPP and LAI created with stand parameters, and sequential predictions were observed in the modeling of NPP with LAI. To eliminate the problems that occurred, the method was rerun by adding early stopping, network reduction and regularization arguments. By using these arguments, we aimed to learn the model data set of the network, to achieve the desired level of modeling success in independent data groups and to create a model network with high representation ability. The findings obtained by adding the relevant arguments are presented in Fig. 9. The overfitting problems encountered in the DL method were improved with the arguments arranged for the model network. In modeling NPP with stand parameters, while the coefficients of determination for model-test were 0.90-0.29, these values increased to 0.64-0.51 with the regulation arguments (Fig. 9). As a result of these findings, a better performance was obtained than with the SVM method. The coefficients of determination for the model-test, which were 0.72-0.46 in the model network applied for the stand parameters with LAI, also reached the level of 0.68-0.53. For the NPP-stand parameters model, DL performed better than the SVM method. However, the SVM method was found to be more successful for the LAI-stand parameters model because its success in the test data set was higher (SVM test R² = 0.64, DL test R² = 0.53). But there was no difference in the performance criteria obtained for NPP and LAI.

Figure 9. Observed and predicted values obtained from NPP-stand parameters, LAI-stand parameters and NPP-LAI models run by adding arguments to the DL method.

DiscussionTop

Relationships of NPP and LAI

The relationships between NPP and LAI are highly affected by environmental and climatic conditions (Pan et al., 2015; Yan et al., 2021). Subtropical forests have a high NPP but a low LAI (Luo et al., 2004). NPP/LAI ratio is higher in arid than humid regions due to high net production and small leaf surfaces (Pan et al., 2021). These conditions also affect the function of the relationship between NPP and LAI. There are studies in which the NPP-LAI relationship is positive in different study areas (Gholz, 1982; Scurlock et al., 2001; Luo et al., 2004; Kushida et al., 2007). Pan et al. (2021) investigated the relationship between NPP and LAI on a global scale and found that there is a parabolic relationship between NPP and LAI. They reported that parabolic relationships between NPP and LAI are rarely in regions with poor site index conditions. Especially, it has been reported that relation of NPP and LAI is parabolic in humid areas and positively related in dryland areas. In this study, sample plot data up to the breaking point made up most of the data set, and the relationship between NPP-LAI was positive in this section. NPP tended to decrease especially in sample plots where LAI was high after the breaking point. This caused the overall relationship between NPP and LAI to be parabolic. This study was carried out in a semi-arid region backward the Black Sea region. Therefore, it can be said that the parabolic relationship between NPP and LAI is due to the semi-arid region.

Evaluation of modeling findings on net primary productivity

NPP was modeled using stand parameters with MRA method. Significant independent variables in this model were htop, N, BA, and V. The coefficients of determination for the model-test were found to be 0.50-0.30. The coefficients of determination for the model-test for the SVM and DL methods were 0.62-0.44 and 0.64-0.51, respectively. In addition to stand parameters, NPP was also modeled with LAI. The coefficient of determination for the model-test was obtained as 0.24 using DL method. The change in the level of success was controlled by testing the regulation arguments in the DL technique. However, no change was observed in model-test successes as a result of the added arguments (model-test R² = 0.24). Bulut et al. (2019) modeled the NPP with stand volume, volume increment and stand carbon separately using the MRA method in their study. The coefficients of determination obtained by stand volume and volume increment were 0.78 and 0.75, respectively. The highest success was obtained with stand carbon (R² = 0.85). Since the independent variables used represent the accumulation in the stand, the correlation levels were obtained at a high rate. Since NPP gives the net produced carbon amount in the stands, it is expected that the highest level of relationship between them will be determined with stand carbon. In this study, stand parameters such as htop, N, BA, and V were used to model NPP. Although a parameter characterizing the accumulation in the stand, such as V, was used in this study, stand parameters were able to predict only 50% of the explainable variance in NPP modeling. Differences between model performances may be due to the different characteristics of the datasets.

Gholz (1982) examined the relations of NPP with BA and LAI in Tsuga heterophylla-Picea sitchensis, Tsuga mertensiana, Pinus ponderosa and Juniperus occidentalis stands. The coefficient of determination of the MRA model between NPP and BA was found to be 0.95. This value was 0.96 between NPP and LAI. The findings show that there is a high level of linear relations between the parameters. There are quite high differences between our study and Gholz´s (1982). In our study, MRA method and stand parameters (h_top, N, BA, and V) were used in the modeling of NPP and a determination coefficient of 0.50 was obtained. When LAI was used to model the NPP, a determination coefficient of 0.20 was found. The different results between Gholz (1982) and our study may be due to the number of samples used. As the number of data increases, the error in the model is expected to decrease (Sterenczak et al., 2018). Gholz (1982) examined the relationships with 8 sample data taken from 12 different main vegetation types, and the model success was high. However, the low number of data also reduces the representativeness of the model.

The performances of parametric MRA method and non-parametric SVM and DL methods in modeling NPP with LAI were found to be close to each other, because the correlation between NPP and LAI was not high (r = 0.47). Artificial intelligence techniques may be insufficient in terms of modeling success among data that do not have a significant relationship and trend (Aertsen et al., 2010). This may be due to the inclusion of all data from sample plots that different stand structures in the modeling (Çil, 2014). To see this effect, the data were modeled separately according to stand developmental stages and crown closure classes. When classification was made according to stand developmental stages and crown closure, which have a great effect on NPP and LAI, the strength of the relations between the two increased and affected the success of the model. The highest correlation in the classes devoted to stand developmental stages was obtained for the stands at the “a-b” developmental stages (R² = 0.70, r = 0.84). In the crown closure classes, although the highest relationships were close to each other, 1 (R² = 0.38, r = 0.62) and 3 (R² = 0.37, r = 0.61) closed stands gave better results. As a result, it was observed that the success rate increases when the parameters that do not have a direct strong relationship with each other are grouped in terms of the factors that are decisive for them. It is important for such groupings that the number of sample plots in the inventory studies is similar in terms of the study-specific criteria. Thus, each class can have enough data for modelling. In particular, the fact that the number of data belonging to the classes is separable for model-testing also allows the application of machine learning algorithms.

Evaluation of modeling findings on LAI

In the modeling in which stand parameters were used as the independent variable for LAI, the significant parameters selected according to the MRA method were dg, N and BA. Among the applied techniques, the SVM method showed the highest test success. The error of the estimation made by the DL method is 0.407 m² m^-2 after adding arguments, while the error of the SVM method was 0.371 m² m^-2. Considering the model performances, the DL method (Model R² = 0.68) is more successful than the SVM method (Model R² = 0.64). Khosravi et al. (2012) made an estimation of LAI for oak forests using stand parameters. The coefficients of determination obtained were 0.36 with dg, 0.36 with BA, and 0.45 with mean stand height (h). Ercanlı et al. (2018) determined the LAI for 108 sample plots taken from pure Anatolian black pine stands by hemispherical photographic acquisition. Stand parameters were also measured in the sample plots and the LAI was modeled using MRA and ANN techniques. The coefficients of determination obtained for the MRA, and ANN methods are 0.5431 and 0.6392, respectively. With the ANN method, the estimation was improved, and the error amount was reduced from 0.3935 m² m^-2 to 0.3497 m² m^-2.

In our study, coefficient of determination of 0.48 was obtained for LAI model with the MRA method. Coefficients of determination of LAI model by using machine learning algorithms such as SVM and DL methods were 0.64 and 0.68, respectively. It has been observed that more accurate predictions can be made with nonlinear techniques. Here, if the aim is to analyze linear relationships and interpret simple relationships, MRA method may be appropriate. Linear models produced with easy-to-measure parameters that can be obtained from stands are considered useful in estimating LAI. In cases where model success is important in estimating LAI, non-parametric techniques such as SVM and DL can be used. Machine learning methods are thought to be useful in making predictions and mappings with higher success rates and usable levels. In addition to the applied techniques, the independent variables used also affect the model success and representation.

In our study, the relationships between LAI and stand characteristics were investigated. Considering that the stand age will increase as the diameter increases, there is a decrease in the LAI values inversely proportional to the stand age. In young stands, an increase in LAI may occur in the early years during the canopy formation process (Jagodzinski & Kalucka, 2008). However, a gradual decrease in LAI occurs with the increase of stand age (Clough et al., 2000; Kushida et al., 2007). However, it is insufficient to use only stand features in linear models to be produced, because besides stand characteristics, all edaphic conditions affect allometric relations with LAI. Therefore, the use of soil properties, which are edaphic factors, as well as stand parameters in linear models to be produced for LAI estimation, contributes to the success and representativeness of the model (Gower et al., 1999; Asner et al., 2003).

Evaluation of overall findings

Model-test successes were compared for the overfitting problem, which is frequently encountered in SVM and DL methods. According to the findings we obtained, the model-test success levels for the SVM method were found to be compatible with each other. NPP-stand parameters (Model R² = 0.90, Test R² = 0.29) and LAI-stand parameters (Model R² = 0.72, Test R² = 0.46) models were found to have an overfitting problem due to the difference between model and test success. Some of the reasons for this problem are the insufficient training data set, poor representation ability and high outlier data rate (Ying, 2019). To eliminate the problem, early stopping (Raskutti et al., 2014), network reduction (Bramer, 2002) and regularization (Warde-Farley et al., 2013) arguments were added to the model (Candel et al., 2016; MathWorks, 2019). As a result of this process, the effects of the overfitting problem were largely eliminated. Findings obtained were model R² = 0.64 and test R² = 0.51 for NPP-stand parameters, model R² = 0.68 and test R² = 0.53 for LAI-stand parameters.

In the DL method, the trend of the relationship between dependent and independent variables is determined during the training process of the network. When this trend which is established with the existing data reflects the reality, the model accuracy rate and the representation ability can be sufficient. Thus, when the trained model network is applied to independent data groups, satisfactory results can be obtained. In this case, the characteristics of the data sets to be used for modeling are of great importance. There is no obvious trend in the modeling of parameters that do not have any organic link between them (Aertsen et al., 2010). It is very difficult to generalize the model network created in this way. Thus, a consistent estimation cannot be made with the independent data run in the model network created. Although the success in the training process is high, the problem of overfitting, which is frequently encountered in machine learning algorithms, arises. As a result, it causes inconsistent estimates that are not representative and do not allow the application of independent data sets (Ying, 2019).

Overfitting problem was encountered in the models produced for NPP and LAI. NPP is a parameter that is estimated using many data groups (photosynthetic data, meteorological data, remote sensing data), while LAI is a parameter that can vary greatly according to stand structures. Since there is no direct relationship and connection between the independent variables, a clear trend could be formed in the trained network. Thus, the test data also fails in the generated model network. The application of these methods in parameters such as age and site index (Ercanlı et al., 2018), which have a clear relationship, can reduce the incidence of this problem and increase model-test successes to higher levels.

With the application and development of artificial intelligence techniques in different fields, modeling studies with high success rates are carried out. Although these techniques have high success potential, they are not simple and user-friendly (Aertsen et al., 2010). In artificial intelligence studies, it is very difficult and time-consuming to find the necessary arguments and the most appropriate number of neurons and layers for the training process of the model network. In addition, the fact that it has uninterpretable weights originating from the hidden layer and nonlinear activation functions limits the interpretation of ecological processes (Peng & Wen, 1999). As a solution to this problem, there are techniques in the literature to determine the contribution of independent variables (Maier & Dandy, 2000). The use of simpler, transparent, and interpretable techniques can be useful and explanatory to reveal the relationships in the parameters related to the productivity and development of stands such as NPP and LAI.

ConclusionsTop

In this study, the relationships between NPP-LAI, NPP-stand parameters and LAI-stand parameters using MRA, SVM and DL methods in pure Anatolian black pine stands were investigated. A non-linear parabolic relationship was found between NPP and LAI. Compared to the MRA method, the SVM and DL were found to be more successful in modeling NPP-LAI, NPP-stand parameters, and LAI-stand parameters. These models, especially developed with SVM and DL modeling techniques, may considerably contribute to the understanding of stand structure and dynamics. As a result, they will contribute to sustainable forest management planning, developing forest growth and yield models, and determine the most suitable silvicultural prescriptions to the forest stands.

RESEARCH ARTICLE

Assessment of the interactions among net primary productivity, leaf area index and stand parameters in pure Anatolian black pine stands: A case study from Türkiye

Introduction Top

Material and methods Top

ResultsTop

DiscussionTop

ConclusionsTop

AcknowledgmentsTop

Credit author statementTop

References Top