Modeling the Radial Stem Growth of the Pine (Pinus sylvestris L.) Forests Using the Satellite-Derived NDVI and LST (MODIS/AQUA) Data

The paper considers a new approach to modeling the relationship between the increase in woody phytomass in the pine forest and satellite-derived Normalized Difference Vegetation Index (NDVI) and Land Surface Temperature (LST) (MODIS/AQUA) data. The developed model combines the phenological and forest growth processes. For the analysis, NDVI and LST (MODIS) satellite data were used together with the measurements of tree-ring widths (TRW). NDVI data contain features of each growing season. The models include parameters of parabolic approximation of NDVI and LST time series transformed using principal component analysis. The study shows that the current rate of TRW is determined by the total values of principal components of the satellite indices over the season and the rate of tree increment in the preceding year.


Introduction
Research addressing the growth and changes in productivity of boreal forests is useful for understanding and predicting sustainability of these ecosystems affected by diverse environmental stresses associated with extreme weather events and climate change [1][2][3][4][5][6][7]. Most studies dealing with this problem predict that, because of climate change, tree growth will be affected by more frequent and violent extreme events (ice storms, wetland formation, droughts, windblows, etc.). However, there is considerable uncertainty about the processes underlying tree response to the increase in the frequency of extreme events.
It is well known that tree growth results from complex physiological processes in tree organs, and availability of resources for tree growth is largely determined by climate [8]. Tree growth processes and weather events can be understood using methods of dendrochronology, which studies time series of tree ring widths (TRW) over many years [9][10][11][12].
Data on long-term TRW variations are used to determine how much wood is formed in the tree stem, which can be directly associated with the biomass increase and carbon consumption [13][14][15][16]. In recent decades, dendrochronology methods have been used to study interactions between species-specific growth and temperature and humidity variations at the local, regional, and continental scales [17][18][19][20][21].
Although the data on tree rings provide the basis for accurate and retrospective estimates of tree growth, to collect wood samples in the field and to measure them in the laboratory are painstaking tasks. Therefore, the use of this approach in real-time monitoring of forest growth at large spatial scales is limited [22]. Satellite remote sensing data, however, are successfully used to face these challenges. Over the past 30 years, extensive research efforts have been devoted to determining biophysical variables of the forest based on satellite measurements [23][24][25][26][27][28][29][30].
Normalized Difference Vegetation Index (NDVI) is the most commonly used vegetation index, which measures fractional absorbed photosynthetically active radiation [31]. The data derived from MODIS (Moderate-Resolution Imaging Spectroradiometer aboard the Aqua satellite) observations are used to produce NDVI time series at a spatial resolution of 0.25 km. However, the time period covered by these data is too short (since 2000) compared to the far longer series of the data obtained from tree rings.
Various vegetation indices determined using satellite data are effective indicators of tree photosynthetic activity and productivity, and a number of studies reported attempts to relate directly the data on TRW with the satellite data [32][33][34][35][36][37][38]. Those studies found a moderate positive correlation between the NDVI interannual variations and annual growth of the tree. There were other studies, though they did not find any significant positive correlation [35,36,[39][40][41]. The relationship between the satellite vegetation indices and TRW is not obvious; it is determined by the species composition of the forest, landscape, and the choice of the proper initial satellite and meteorological data [42].
In most of the studies cited above, the authors analyzed correlations between the time series of TRW data and NDVI composite data for certain periods (e.g., 2 weeks or 1 month). The use of statistical methods, though, is not the best approach to revealing the relationship between NDVI and tree growth. For instance, annual wood formation is more strongly related to the net primary production (NPP) calculated over long periods, as the increase in tree biomass is not directly determined by the rate of photosynthesis but is rather associated with stems, branches, and roots growing and forming wood [43,44]. Thus, the storage of the long-lived phytomass in the woody parts of the trees is physically and temporally separated from the growth of leaves, identified by NDVI. Therefore, the different climatic conditions necessary for wood formation and for photosynthetic activity may disconnect the NDVI-growth relationship, as higher temperatures are needed for growth than for photosynthesis [45,46]. However, in boreal forests, NDVI indicates the physiological status of the trees, and, thus, the state of the tree crowns in summer (June and July) is a more important factor for tree growth than the length of the growing season [47].
As regards tree growth research, it can be conducted at a tree level, using, e.g., dendrochronology methods, or at a level of the tree stand or plot-to determine total biomass of trees or stands using various remote sensing techniques [48]. However, in the attempt to combine both approaches, researchers are faced with the challenge of differences in the time scales of crown formation processes, leaf photosynthesis (observed by satellites), and wood formation processes (integrated in tree rings). The response of photosynthesis to weather is relatively instantaneous, but there is a well-documented time lag between carbon absorption and increase in biomass of photosynthetic carbon [49]. Therefore, the TRW data and vegetation indices cannot be expected to be completely consistent with each other, as the space and time scales of TRW time series and satellite indices are different [46,50].
The reason for the spatial mismatch of the scales is that satellite data for NDVI (MODIS) are a pixel of land surface of an area of 6250 m 2 , whereas the TRW data are obtained for one tree with a crown projection area of up to 250 m 2 , and these data may vary considerably within the pixel. Another problem arises from the temporal mismatch between photosynthesis processes observed by satellite and wood formation processes integrated in TRW. For these reasons, no direct correspondence between TRW data and satellite vegetation indices can be expected. These processes will vary between ecosystems, tree species, and climate regions [2][3][4].
The search for interactions between the data on radial stem growth and the satellitederived NDVI and LST (Land Surface Temperature) gives rise to the problem of selecting significant variables for analysis. For each year, radial stem growth is characterized by one value, while annual satellite data on NDVI and LST are time series. For instance, 8-day seasonal composites of NDVI and LST data consist of 40-50 points. It is rather difficult to Atmosphere 2021, 12, 12 3 of 15 find correlations between one number and function. Attempts to choose a certain number of variables in NDVI and LST time series and construct multiple regression equations of relationships between tree ring widths and the selected remote sensing data are hampered by the uncertainty about the choice of remote sensing variables for analysis. The problem is that for the time series of remote sensing data, one can choose from a great number of key variables or introduce additional variables such as the sums of temperatures for a certain period (but which?) of the season. The present study deals with the methods of constructing models of the relationship between the increase in woody phytomass and remote sensing indicators based on using integrated parameters of NDVI and LST time series transformed by principal component analysis.

Materials and Methods
Dynamics of radial stem growth was studied and phenological data were collected on sample plots in pine (Pinus sylvestris L.) forests located around the city of Krasnoyarsk. Trees aged about 100 years were selected to analyze radial stem growth on all sample plots. The trees had the stems of diameter at breast height ranging from 12 to 32 cm, and the average stand density was about 600-800 trees per ha. Time series of radial stem growth between 2003 and 2015 were analyzed. This time interval was chosen because the remote sensing data were only available from 2003. Coordinates of sample plots and characterization of tree stands are provided in Table 1.
An increment borer was used to extract core samples from all trees on the sample plots. The tree ring width d(t) in year t was measured on the polished surfaces of the core samples in the field of view of the microscope in the "Lintab 5 Tree-RingStation" (RINNTECH ® ) (version No. 5) with an accuracy of 0.1 mm. The dating of tree rings was carried out visually, using graphs in the TSAP-WinTM (RINNTECH ® ) software, and control of the dating was performed using the COFECHA program (version 6.0P) from the dendrochronology program library, DPL [51]. Figure 1 shows the positions of sample plots. The search for interactions between the data on radial stem growth and the sa lite-derived NDVI and LST (Land Surface Temperature) gives rise to the problem of lecting significant variables for analysis. For each year, radial stem growth is charac ized by one value, while annual satellite data on NDVI and LST are time series. For stance, 8-day seasonal composites of NDVI and LST data consist of 40-50 points. I rather difficult to find correlations between one number and function. Attempts choose a certain number of variables in NDVI and LST time series and construct mult regression equations of relationships between tree ring widths and the selected rem sensing data are hampered by the uncertainty about the choice of remote sensing va bles for analysis. The problem is that for the time series of remote sensing data, one choose from a great number of key variables or introduce additional variables such as sums of temperatures for a certain period (but which?) of the season. The present stu deals with the methods of constructing models of the relationship between the increas woody phytomass and remote sensing indicators based on using integrated parame of NDVI and LST time series transformed by principal component analysis.

Materials and Methods
Dynamics of radial stem growth was studied and phenological data were collec on sample plots in pine (Pinus sylvestris L.) forests located around the city of Krasnoya Trees aged about 100 years were selected to analyze radial stem growth on all sam plots. The trees had the stems of diameter at breast height ranging from 12 to 32 cm, a the average stand density was about 600-800 trees per ha. Time series of radial s growth between 2003 and 2015 were analyzed. This time interval was chosen because remote sensing data were only available from 2003. Coordinates of sample plots characterization of tree stands are provided in Table 1. The search for interactions between the data on radial stem growth and the sa lite-derived NDVI and LST (Land Surface Temperature) gives rise to the problem of lecting significant variables for analysis. For each year, radial stem growth is charac ized by one value, while annual satellite data on NDVI and LST are time series. For stance, 8-day seasonal composites of NDVI and LST data consist of 40-50 points. I rather difficult to find correlations between one number and function. Attempts choose a certain number of variables in NDVI and LST time series and construct multi regression equations of relationships between tree ring widths and the selected rem sensing data are hampered by the uncertainty about the choice of remote sensing va bles for analysis. The problem is that for the time series of remote sensing data, one choose from a great number of key variables or introduce additional variables such as sums of temperatures for a certain period (but which?) of the season. The present stu deals with the methods of constructing models of the relationship between the increase woody phytomass and remote sensing indicators based on using integrated paramet of NDVI and LST time series transformed by principal component analysis.

Materials and Methods
Dynamics of radial stem growth was studied and phenological data were collec on sample plots in pine (Pinus sylvestris L.) forests located around the city of Krasnoyar Trees aged about 100 years were selected to analyze radial stem growth on all sam plots. The trees had the stems of diameter at breast height ranging from 12 to 32 cm, a the average stand density was about 600-800 trees per ha. Time series of radial st growth between 2003 and 2015 were analyzed. This time interval was chosen because remote sensing data were only available from 2003. Coordinates of sample plots a characterization of tree stands are provided in Table 1. An increment borer was used to extract core samples from all trees on the sam plots. The tree ring width d(t) in year t was measured on the polished surfaces of the c samples in the field of view of the microscope in the "Lintab 5 Tree-RingStatio (RINNTECH ® ) (version No. 5) with an accuracy of 0.1 mm. The dating of tree rings w carried out visually, using graphs in the TSAP-WinTM (RINNTECH ® ) software, and c trol of the dating was performed using the COFECHA program (version 6.0P) from dendrochronology program library, DPL [51]. Figure 1 shows the positions of sample pl For the variables analyzed over the time interval selected in this study, no stro trend of the NDVI and LST data was observed, but there was a temporal trend of TRW typical TRW time series with trend is shown in Figure 2 (curve 1). An increment borer was used to extract core samples from all trees on the sam plots. The tree ring width d(t) in year t was measured on the polished surfaces of the c samples in the field of view of the microscope in the "Lintab 5 Tree-RingStatio (RINNTECH ® ) (version No. 5) with an accuracy of 0.1 mm. The dating of tree rings w carried out visually, using graphs in the TSAP-WinTM (RINNTECH ® ) software, and c trol of the dating was performed using the COFECHA program (version 6.0P) from dendrochronology program library, DPL [51]. Figure 1 shows the positions of sample pl For the variables analyzed over the time interval selected in this study, no stro trend of the NDVI and LST data was observed, but there was a temporal trend of TRW typical TRW time series with trend is shown in Figure 2 (curve 1). An increment borer was used to extract core samples from all trees on the sample plots. The tree ring width d(t) in year t was measured on the polished surfaces of the core samples in the field of view of the microscope in the "Lintab 5 Tree-RingStation" (RINNTECH ® ) (version No. 5) with an accuracy of 0.1 mm. The dating of tree rings was carried out visually, using graphs in the TSAP-WinTM (RINNTECH ® ) software, and control of the dating was performed using the COFECHA program (version 6.0P) from the dendrochronology program library, DPL [51]. Figure 1 shows the positions of sample plots For the variables analyzed over the time interval selected in this study, no strong trend of the NDVI and LST data was observed, but there was a temporal trend of TRW. A typical TRW time series with trend is shown in Figure 2 (curve 1). For the variables analyzed over the time interval selected in this study, no strong trend of the NDVI and LST data was observed, but there was a temporal trend of TRW. A typical TRW time series with trend is shown in Figure 2  The procedure of estimating the relationships between the tree growth and NDVI and LST data consisted of several stages. In Stage 1, the temporal trend of TRW was excluded from the data on radial stem growth. Therefore, for each tree on each sample plot, time series of the TRW first differences (TRW FD) between years were calculated: TRW FD(t) = TRW(t) − TRW(t−1). In subsequent calculations, the TRW FD time series, which were stationary time series, were used instead of the time series of tree ring widths (Figure 1, curve 2). The TRW FD time series was actually the rate of change of TRW. The time series of the first differences were calculated for each tree, and then average first differences were calculated for each sample plot for each year between 2003 and 2015. Table S1 provides data on the average radial stem growth for all sample plots during 1987-2015.
To estimate the relationships between the current value of radial stem growth and the corresponding values for the previous years, these relationships will be considered as the autoregressive equation: where k is the order of the autocorrelation.
To estimate the order of the autocorrelation for the time series of the radial stem growth on the sample plots between 1987 and 2015, partial autocorrelation function (PACF) was calculated. Figure 3 shows the PACF for the TRW FD data on sample plot 1. The procedure of estimating the relationships between the tree growth and NDVI and LST data consisted of several stages. In Stage 1, the temporal trend of TRW was excluded from the data on radial stem growth. Therefore, for each tree on each sample plot, time series of the TRW first differences (TRW FD) between years were calculated: TRW FD(t) = TRW(t) − TRW(t − 1). In subsequent calculations, the TRW FD time series, which were stationary time series, were used instead of the time series of tree ring widths ( Figure 1, curve 2). The TRW FD time series was actually the rate of change of TRW. The time series of the first differences were calculated for each tree, and then average first differences were calculated for each sample plot for each year between 2003 and 2015. Table S1 provides data on the average radial stem growth for all sample plots during 1987-2015.
To estimate the relationships between the current value of radial stem growth and the corresponding values for the previous years, these relationships will be considered as the autoregressive equation: where k is the order of the autocorrelation.
To estimate the order of the autocorrelation for the time series of the radial stem growth on the sample plots between 1987 and 2015, partial autocorrelation function (PACF) was calculated. Figure 3 shows the PACF for the TRW FD data on sample plot 1. The procedure of estimating the relationships between the tree growth and NDVI and LST data consisted of several stages. In Stage 1, the temporal trend of TRW was excluded from the data on radial stem growth. Therefore, for each tree on each sample plot, time series of the TRW first differences (TRW FD) between years were calculated: TRW FD(t) = TRW(t) − TRW(t−1). In subsequent calculations, the TRW FD time series, which were stationary time series, were used instead of the time series of tree ring widths (Figure 1, curve 2). The TRW FD time series was actually the rate of change of TRW. The time series of the first differences were calculated for each tree, and then average first differences were calculated for each sample plot for each year between 2003 and 2015. Table S1 provides data on the average radial stem growth for all sample plots during 1987-2015.
To estimate the relationships between the current value of radial stem growth and the corresponding values for the previous years, these relationships will be considered as the autoregressive equation: where k is the order of the autocorrelation.
To estimate the order of the autocorrelation for the time series of the radial stem growth on the sample plots between 1987 and 2015, partial autocorrelation function (PACF) was calculated. Figure 3 shows the PACF for the TRW FD data on sample plot 1. PACF is significant only at k = 1, i.e., the current value of TRW FD is only determined by its previous value ( Figure 3). Calculations of PACF for other sample plots result in the same value of k = 1, i.e., the current value of TRW FD(t) correlates with the value of TRW FD(t − 1) in the previous year.
In Stage 2, integrated parameters characterizing NDVI and LST were considered. In analysis of remote sensing data, NDVI was determined using the standard formula as: where NIR (0.620-0.670 µm) and Red (0.841-0.876 µm) are normalized values of the intensity of near-infrared and red light reflected by the land surface at a given point. NDVI values were determined using the MYD09Q1 product of the MODIS (or Moderate Resolution Imaging Spectroradiometer), a key instrument aboard the Terra (originally known as EOS AM-1) and Aqua international research satellite. The MYD09Q1 product included values of the red and infrared components of the reflected radiation for a 250 m × 250 m pixel. For each pixel, an 8-day composite was used, taking into account atmospheric conditions. The data for 2003-2017 were obtained from the Earth Observing System Data and Information System EOSDIS resource [52]. The data on the land surface temperature (LST) were obtained for the same period (the MYD11A2 product has 8-day composite at a spatial resolution of 1 km).
Complete seasonal data of the MYD11A2, MYD09Q1, MOD11A2, and MOD09Q1 products on the Terra and Aqua satellites are available for the study area from 2003.
We assumed that NDVI data contain features of phenological dynamics of each growing season [53] and LST data characterize weather changes. NDVI and LST seasonal dynamics will be described using nonlinear regression equations, y(n) = a 1 n 2 + a 2 n + a 3 , where n is the first 8-day composite from the beginning of the year. This approximation is possible for boreal forests if the winter period is excluded. The annual curve of NDVI is used to find the n 1 and n 2 dates, which correspond to the beginning and end of the growing season determined for the study tree stands (when NDVI was stably higher than 0.35 of the NDVI maximum). Approximation of the NDVI and LST time series is done on the (n 1 − n 2) section. The NDVI time series was preliminarily smoothed using a rectangular filtering window of width of 3 composites to remove measurement artifacts.
A typical seasonal NDVI time series is shown in Figure 4 and an LST seasonal time series in Figure 5. Quadratic regression equations of the NDVI and LST time series for the growing season were calculated for each sample plot using the data for the 2003-2015 period. The accuracy of the parabolic approximation is described by the coefficient of determination, R 2 , which characterizes the proportion of the variance in the time series of remote sensing data accounted for by regression. Parameters of approximation of the NDVI and LST seasonal curves NDVI(k,j,n) = a 1 (k,j)n 2 + a 2 (k,j)n − a 3 (k,j) and LST seasonal curves LST(k,n,j) = b 1 (k,j)n 2 + b 2 (k,j)n − b 3 (k,j) for sample plot k (k = 1, . . . 4), in season j (j = 2003, . . . , 2015) and dates n can be used to estimate the relationships between the remote sensing data and the growth of the tree ring width. Typical time series of NDVI and LST are shown in Figures 4 and 5.
Calculations of parameters of the NDVI and LST seasonal time series and validation of parabolic models using Fisher's test for all study years and all sample plots are provided in Tables S2 and S3.          Table 2 provides characterization of principal components for sample plot No. 1. Based on Table 2, the value of the first principal component, G1(t), can be expressed in the following equation: The values of other principal components can be determined in the same way from Table 2. Table 3 provides contributions of the first eigenvalues to the total NDVI and LST variance for all sample plots. In all cases, the first principal components for NDVI and LST account for at least 72% of variance, and the first two principal components taken together account for at least 96% of variance (Table 3); thus, two variables, G1 and H1, or four variables, G1, G2, H1, and H2, can be used instead of eight. Figure 6 shows the curves indicating temporal dynamics of the first principal components G1(t) for the NDVI matrices (curve 1), H1(t) for LST matrices (curve 2), and the time series of the first differences of radial stem growth, TRW FD(t) (curve 3) for sample plot 1.
After such transformations, one can construct linear regression models of relationships between the first differences of radial stem growth and values of the principal components of matrices A and B. The quality of the proposed models can be assessed by t-test for coefficients of the model, coefficient of determination, R 2 , of the model series, and cross-correlation function of the TRW FD time series and the model series, characterizing synchrony of change of these series. Results of the calculations were used as the basis for selecting variables for constructing models of the relationship between TRW FD, NDVI, and LST. After such transformations, one can construct linear regression models of relationships between the first differences of radial stem growth and values of the principal components of matrices A and B. The quality of the proposed models can be assessed by t-test for coefficients of the model, coefficient of determination, R 2 , of the model series, and cross-correlation function of the TRW FD time series and the model series, characterizing synchrony of change of these series. Results of the calculations were used as the basis for selecting variables for constructing models of the relationship between TRW FD, NDVI, and LST.

Results
To construct the model of radial growth, we take into account the relationship of the current year growth, TRW(t), to the growth of the previous year, TRW(t−1), determined by calculating PACF, and the effect of the weather will be characterized by the values of the first principal components, G1(t) and H1(t). This model was tested using the ADL (autoregressive distributed lag) model: (2) As model (1) had passed the test, it was used to make calculations for all sample plots over the period considered in this study. Results are given in Table 4.

Results
To construct the model of radial growth, we take into account the relationship of the current year growth, TRW(t), to the growth of the previous year, TRW(t − 1), determined by calculating PACF, and the effect of the weather will be characterized by the values of the first principal components, G1(t) and H1(t). This model was tested using the ADL (autoregressive distributed lag) model: TRW FD(t) = c 0 + c 1 G1(t) + c 2 H1(t) + c 3 TRW FD(t − 1). (2) As model (1) had passed the test, it was used to make calculations for all sample plots over the period considered in this study. Results are given in Table 4. Cross-correlation function, CCF, was constructed to test time synchrony of the field and model data (Figure 7). The CCF maximum was observed at the delay k = 0, suggesting the high degree of synchrony between the field and model data series.
Almost all coefficients for sample plots were significant ( Table 5). The significance level of p < 0.1 was used, as some of the factors influencing phytomass increase (precipitation and soil conditions) cannot be estimated using remote sensing measurements, although they do affect the phytomass increase. Although satellite microwave data ASMR-E [54] can be used to estimate precipitation and soil moisture, the low spatial and temporal resolution and accuracies of such data are not appropriate for this research.
Cross-correlation function, CCF, was constructed to test time synchrony of the field and model data (Figure 7). The CCF maximum was observed at the delay k = 0, suggesting the high degree of synchrony between the field and model data series. Almost all coefficients for sample plots were significant ( Table 5). The significance level of p < 0.1 was used, as some of the factors influencing phytomass increase (precipitation and soil conditions) cannot be estimated using remote sensing measurements, although they do affect the phytomass increase. Although satellite microwave data ASMR-E [54] can be used to estimate precipitation and soil moisture, the low spatial and temporal resolution and accuracies of such data are not appropriate for this research.
Susceptibility of the mean radial stem growth to the change of the first principal component for the NDVI data is positive, i.e., as G1 increases, TRW FD(t) increases as well. By contrast, susceptibility of TRW FD(t) to H1 is negative. For almost all sample plots, susceptibility of TRW FD(t) to the change of the mean radial growth in the preceding year is negative.
The LST and NDVI data for the period of tree dormancy were excluded from the time series of remote sensing data and were not used in calculations. Only the growing season data were taken into account. The reduction in NDVI and LST date is needed to make it possible to compare the time series of TRW with time series of remote sensing data.
Thus, when only satellite data were used to investigate growth processes in tree stands, no direct relationship was found between integrated parameters of NDVI and LST and radial stem growth rate.

Discussion
The signs of coefficients of Equation (2), which characterize susceptibility of the first differences of radial growth to changes in the values of independent variables in model  Susceptibility of the mean radial stem growth to the change of the first principal component for the NDVI data is positive, i.e., as G1 increases, TRW FD(t) increases as well. By contrast, susceptibility of TRW FD(t) to H1 is negative. For almost all sample plots, susceptibility of TRW FD(t) to the change of the mean radial growth in the preceding year is negative.
The LST and NDVI data for the period of tree dormancy were excluded from the time series of remote sensing data and were not used in calculations. Only the growing season data were taken into account. The reduction in NDVI and LST date is needed to make it possible to compare the time series of TRW with time series of remote sensing data.
Thus, when only satellite data were used to investigate growth processes in tree stands, no direct relationship was found between integrated parameters of NDVI and LST and radial stem growth rate.

Discussion
The signs of coefficients of Equation (2), which characterize susceptibility of the first differences of radial growth to changes in the values of independent variables in model (2), differ between sample plots. For example, for H1 on all sample plots (except No. 3), the variables have the negative sign, i.e., the larger G1, the lower the first differences of radial growth.
It could be assumed that TRW FD, D(t), might be related to not only the TRW FD value of the preceding year, D(t − 1), but also phenological and weather parameters of the preceding year G1(t − 1) and H1(t − 1). However, for model TRW FD(t) = c 0 + c 1 G1(t − 1) + c 2 H1(t − 1) + c 3 TRW FD(t − 1), coefficients c 1 and c 2 are nonsignificant even at p = 0.50, and the proportion of variance accounted by this model is only 0.30, i.e., the current TRW FD is not dependent on these parameters.
As the sample plots are located rather close to each other, the tree stands have comparable species compositions, and weather conditions are similar as well; therefore, the average values of G1, H1, and TRW FD(t) can be used in the model. Parameters of the averaged model (2) are listed in Table 5.
Coefficients of model (2) for averaged data are significant at p < 0.03, coefficient of determination reaches 0.70, and CCF(k = 0) = 0.84, i.e., the field and model data series are practically time synchronous (Table 5).
Both amplitudes and phases of the averaged field data and model data are in good agreement (Table 5). For some of the sample plots, model (2) accounts for more than 50% of the variance of long-term data (Table 4). For the data averaged over all sample plots, model (2) accounts for 70% of the data variance. Coefficients of Equation (2) are significant for both individual sample plots and the averaged data, characterizing susceptibility of growth processes to changes in independent variables of this equation.
To determine properties of tree assemblages on the sample plot, it is important to assess temporal synchrony between tree growth processes in the assemblage. The synchrony of changes in TRW FD values of individual trees could be assessed using the CCF(x,y) matrix of the series of radial growth first differences between trees x and y. Radial growth can be regarded as synchronous if the highest values of CCF in this matrix are close to 1 and are observed at delay k = 0. For a sample of 25-30 trees, it would be necessary to calculate about 400 paired CCF, and, thus, it may be more convenient to use the correlation matrix of the first differences of radial stem growth. The correlation coefficient r(x,y) of the TRW FD series for trees x and y corresponds to the value of CCF(x, y) at delay k = 0. If r(x,y) is close to 1, these time series can be considered as synchronous. The average value of correlation coefficients r(x,y) or the function of distribution f (r) of correlation coefficients for the tree assemblage can be used as an indicator of the synchrony of radial stem growth in the tree assemblage. For the tree assemblage with the high level of synchrony, the function of distribution of correlation coefficients will peak close to r = 1. Figure 8 shows the function of distribution f (r) of correlation coefficients for different sample plots.  Choosing relevant parameters for models connecting ecological characteristics (such as growth processes) with the environmental factors (such as weather) is always a challenge. The choice is difficult because of the diversity of the parameters, e.g., temperature and precipitation for a certain period of the season or the totals for the entire season. The approach to choosing independent model parameters proposed in this study makes it Choosing relevant parameters for models connecting ecological characteristics (such as growth processes) with the environmental factors (such as weather) is always a challenge. The choice is difficult because of the diversity of the parameters, e.g., temperature and precipitation for a certain period of the season or the totals for the entire season. The approach to choosing independent model parameters proposed in this study makes it possible to limit the number of parameters to the values of two first principal components of approximated NDVI and LST series, considerably simplifying the analysis. The test of the accuracy of a model can be based on, first, coefficient of determination, R 2 , characterizing the similarity between the field and model data by amplitudes, and, second, cross-correlation function, characterizing the synchrony between the dynamics of the field and model time series.
The relationships between remote sensing data (NDVI) and radial stem growth were studied by a number of researchers [39,47]. For example, NOAA AVHRR PAL time series were used to analyze evolution of vegetation in the Komi Republic (the Northwest of Russia) between 1982 and 2001. The authors revealed statistically significant correlation (R 2 = 0.44-0.59) between the Normalized Difference Vegetation Index (NDVI) and tree ring width [37]. At the same time, comparison of the chronology of tree ring width (TRW) and remote sensing data did not reveal any significant correlation between NDVI and TRW [35].
There may be various reasons for the disagreement between the conclusions reached by different researchers. For instance, most of the studies used the NDVI data obtained for a 64-km 2 pixel, and a synchronous NDVI signal can hardly be detected on such an area. The forest vegetation on the study plot was not examined for uniformity, and the synchronicity of the time series of the radial stem growth was not estimated for individual trees. However, the present study demonstrates that as the level of synchronicity of the radial stem growth of individual trees decreases, correlation between the NDVI and TRW data becomes weaker. Moreover, analysis of the correlation between the time series of NDVI and TRW seasonal dynamics should take into consideration that if the dynamics are similar, correlation may be strong while there may be no functional relationship between the parameters.
The present study considered uniform tree stands on a 250 × 250 m 2 pixel and estimated the synchronicity of the time series of radial stem growth in individual trees on the sample plot. Results of the study suggest that even on the small pixel area, the relationship between NDVI and TRW may depend on the synchronicity of the processes of radial stem growth, and both strong and rather weak correlations can be obtained between NDVI and TRW. Thus, research on relationships between woody plant growth and remote sensing data should take into account specific spatiotemporal dynamics of tree growth.

Conclusions
The model proposed in this study relates the average value of TRW FD for the tree stand to the current seasonal integrated remote sensing parameters NDVI and LST. Parameters NDVI and LST were determined from the curves of the temporal dynamics of time series using transformation of principal components. However, no direct relationship was revealed between TRW FD and remote sensing parameters for the tree stands on the study sample plots. The current modelled value of TRW FD was only determined by the total principal components of remote sensing parameters for the season and TRW FD of the preceding year. The first-order ADL-model, i.e., the model with the discrete-time Markov process, was used for analysis. This model is able to find the relationships between tree growth and NDVI and LST.
The present study shows that approximately 50-70% of the current radial stem growth is determined by the growth of the previous year (i.e., has Markov properties) and integrated properties of NDVI and LST. Both the current year growth and the growth of the previous year are certainly affected by weather factors, but the effect of the weather of previous years in the proposed ADL equation is actually represented by the parameter of radial stem growth of the previous year. Moreover, rather than determining the period in the season when the effect on the radial stem growth is the greatest (which may be a problem), one can use a very simple procedure of calculating principal components of the entire seasonal curve. The value of (1 − R 2 ) can be used to estimate the contribution of other, unknown, components to the growth. In such calculations, it is important to take into account the synchrony of radial stem growth processes.