Phenological Response of an Evergreen Broadleaf Tree, Quercus acuta, to Meteorological Variability: Evaluation of the Performance of Time Series Models

Abstract

Phenological events are key indicators for the assessment of climate change impacts on ecosystems. Most previous studies have focused on identifying the timing of phenological events, such as flowering, leaf-out, leaf-fall, etc. In this study, we explored the characteristics of the green chromatic coordinate (GCC) values of the evergreen broadleaf tree (Quercus acuta Thunb.), which is a widely used index that serves as a proxy for the seasonal and physiological responses of trees. Additionally, we estimated their relationship with meteorological variables using time series models, including time series decomposition and a seasonal autoregressive integrated moving average with exogenous regressors (SARIMAX). Our results showed that the GCC values and the meteorological variables, which were collected at daily intervals, exhibited a strong autocorrelation and seasonality. This suggests that time series analysis methods are more suitable than ordinary least squares (OLS) regression methods for the fulfillment of statistical assumptions. The time series analysis results highlighted a strong association between precipitation and GCC variation in evergreen broadleaf trees, particularly during the dry season. These results improve our understanding of the response of plant phenology to climate change.

1. Introduction

The importance of plant phenological events in the context of climate change is undeniable [1], particularly in regions with distinct seasonality, as phenological mismatches can have severe consequences for agriculture and ecosystems [2]. It is well known that the timing of biological events is strongly controlled by meteorological variations. A previous study reported that the flowering timing of 90% of 243 studied plants was significantly related to temperature [3], and other researchers also argued that temperature is one of the most important factors for plant phenology [4,5]. However, few studies have assessed the importance of multiple meteorological factors such as temperature, precipitation, and solar radiation on plant phenological responses.

Traditional plant phenologic observations, including flowering, leaf-out, leaf-falling, etc., have been recorded by repeated visual observations. Capturing the exact date of key events is relatively labor-intensive and costly. Recently, repeated photography taken with ground-based cameras has been used to document changes in phenological variation. Specifically, the calculated RGB (red, green, and blue) bands can be manipulated in numerous ways to detect patterns and classify objects. RGB images can be captured inexpensively and effectively estimate percent plant cover and biomass [6], and the phenocams can serve as a bridge between satellite remote sensing and field observation methods in plant phenology research [7,8]. The interpretation of time series RGB bands can be improved by incorporating meteorological data collected alongside them [9], and time series statistical models can help identify phenological triggers hidden within these datasets.

Time series data consists of sequential observations taken at equally spaced intervals. Time series analysis is based on the fundamental assumption that past patterns will continue to repeat in the future [10]. The most commonly used method for time series analysis is an autoregressive integrated moving average model (ARIMA). This time series model is frequently applied in various fields such as meteorology, epidemiology, limnology, etc. [11,12,13,14], but it is seldom applied in the terrestrial ecological field [15]. In particular, ground-based phenology datasets have rarely been analyzed with time series statistical methods.

The primary objective of this study is to evaluate the applicability of ground-based digital photography and time series analysis as tools for estimating relationships between tree phenological responses and meteorological factors. The first specific objective is to demonstrate the utility of time series models for ground-based phenocam datasets. The second specific objective is to identify which meteorological variables, namely temperature, precipitation, and solar radiation, have an influence on the GCC variability of evergreen broad-leaved forests.

2. Data and Methods

2.1. Study Site and Data Collection

Mt. Sangwang is part of Wando-gun, Jeolla province (34.3° N, 126.6° E) (Figure 1). Mt. Sangwang is located southeast of the Korean Peninsula and belongs to the humid subtropical climate zone (Cfa) according to the Köppen climate classification. Our study site supports one of the largest evergreen broadleaf forests in the Korean Peninsula. Dominant tree species are Quercus acuta Thunb., Castanopsis sieboldii (Makino) Hatus., and Quercus glauca Thunb. Its annual average temperature is 13.6°C and the annual average precipitation was nearly 1310 mm over the study period. It has more than 60% of its precipitation during the summer (monsoon season), and winter and spring suffer from drought (Figure 2). Our site is located 240 m above sea level.

The phenocam used in this study was positioned to capture the southwest-facing canopy of a Quercus acuta forest from the ecological research tower of the National Institute of Forest Science (34°21′32.6″ N, 126°40′22.4″ E). The camera, a Nikon D5600 DSLR, was connected to a modem for remote data transmission, and a plastic housing was used to protect the camera and other electronic equipment from harsh environmental conditions. To obtain optimal daily green chromatic coordinate (GCC) values for phenological analysis, we selected one image per day throughout the study period from September 2017 to August 2023. Typically, to avoid shadows and minimize light scatter, this was taken between 12:00 and 13:00. Regions of interest (ROI) were defined for the evergreen broadleaf tree, Q. acuta, as described by [16,17]. The distance from the camera to the ROI is 8 m. The vegetation index of the GCC was calculated by extracting the relative brightness of the red, green, and blue color channels (RGB), a process which is widely used to identify phenological phase changes [18,19].

$$GCC={\displaystyle \frac{{DN}_G}{{DN}_R+{DN}_G+{DN}_B}}$$

where DN is the digital number and R, G, and B denote the red, green, and blue channels. The GCC value acts as a proxy for understanding tree physiological processes, including chlorophyll content, photosynthetic activity, and carbon assimilation. These processes were performed using the R package ‘phenopix’ [20].

The daily meteorological data was provided by the weather station of the Korea Meteorological Administration [21]. Meteorological factors included mean temperature (Tmean), maximum temperature (Tmax), minimum temperature (Tmin), range temperature (Trange), mean humidity (Hmean), maximum humidity (Hmax), minimum humidity (Hmin), precipitation (precip), and solar radiation (radiation) (

Table 1. Summary statistics of the daily observed variables during the sampling period from September 2017 to August 2023.

Variables (Unit)	Mean	Std	Min	Max	Ljung–Box Test		ADF Test
					Statistics	p-Value	Statistics	p-Value
GCC	0.38	0.022	0.33	0.46	47,040	0.000	−3.80	0.003
Tmean (°C)	13.69	8.38	−8.08	28.5	48,094	0.000	−2.05	0.265
Tmax (°C)	16.75	8.08	−5.47	32.6	46,227	0.000	−2.51	0.113
Tmin (°C)	11.05	8.82	−11.4	26.4	48,873	0.000	−2.54	0.105
Trange (°C)	5.70	1.96	0.36	14.0	5492	0.000	−4.86	0.000
Hmean (%)	75.4	14.3	28.9	100	10,533	0.000	−3.83	0.003
Hmax (%)	92.0	9.86	38	100	4278	0.000	−4.73	0.000
Hmin (%)	57.1	18.2	13	100	10,975	0.000	−3.60	0.006
Precipitation (mm)	3.59	13.3	0	165.5	124	0.001	−38.91	0.000
Solar radiation (MJ/m2)	14.67	7.28	0.91	30.6	7310	0.000	−4.49	0.000

). The meteorological observation station was located 4.4 km away from the phenocam observation site. (Figure 1).

2.2. Statistical Analysis

To estimate the relative importance of meteorological variables on the greenness of evergreen broad-leaved forests, we applied several time series analysis methods. Before statistical analysis, the gap-filling process was carried out based on the spline interpolation method (cubic spline) due to missing 23% of the GCC values [22]. We chose the cubic spline interpolation method because this method is generally more robust for unequally spaced points compared to many other interpolation methods [23]. Missing values are primarily caused by mechanical failures in data collection or storage. We also standardized all the collected variables (Z-Score normalization) to compare these variables on the same scale. To examine the characteristics of the data, the Ljung–Box test was used to assess autocorrelation, and the ADF test was conducted to evaluate stationarity in our datasets [24].

As time series data have many underlying patterns, it is necessary to split the time series data into trend, season, and residual components using seasonal-trend decomposition using LOESS (STL) [25]. After decomposing the time series data, we extracted the trend component to compare the distance between the GCC index and meteorological variables using the dynamic time warping (DTW) method [26]. By focusing on the underlying trends rather than short-term fluctuations or noise, the trend component distance helps explore how two different time series evolve together over a long-term period [27].

The basic time series models, ARIMA models, can be divided into autoregressive (AR) and moving average (MA) terms, which are defined as ARIMA (p, d, q), where p and q are the orders of AR and MA models, respectively, and d represents the degree of series integration (differencing) [28,29]. The SARIMA model can analyze seasonal lags by incorporating seasonal parameters (P, D, Q)^S. Additionally, the SARIMAX model can incorporate exogenous variables to reflect the structural changes of the analyzed process. Its general equation can be defined as

$$\begin{gathered} Y\left(t\right)=c+\phi 1Y\left(t-1\right)+\dots +\phi pY\left(t-p\right)+\theta 1\mathsf{ϵ}\left(t-1\right)+\dots +\theta q\mathsf{ϵ}\left(t-q\right) \\ +\phi a\left(Y\left(t-s\right)-Y\left(t-s-m\right)\right)+\beta X\left(t\right)+\epsilon \left(t\right) \end{gathered}$$

where $Y\left(t\right)$ is the value of the time series at time $t$, $c$ is the constant term, $\phi 1$ to $\phi p$ represent the autoregressive parameters, $\theta 1$ to $\theta q$ represent the moving average parameters, $\beta X\left(t\right)$ represents the exogenous variables at time t and $\beta$ is a vector of coefficients for those variables, $\epsilon \left(t\right)$ is the error term at time $t$, $s$ is the seasonal period, and $m$ is the number of seasonal periods.

Our SARIMAX model incorporates meteorological factors as exogenous variables to estimate the effects of environmental variables on the GCC variations. Owing to computation expenses and to minimize the effects of outliers, the weekly mean time series values were used for the SARIMAX model. In this model, the periodicity parameter S was set to 52 weeks (1 year). The procedure for defining the structure of SARIMAX was as follows: (1) Multicollinearity was investigated before entering the exogenous variables into the SARIMAX model. (2) An augmented Dickey–Fuller test (ADF) was performed to verify the stationary and the periodical terms were estimated using the autocorrelation function (ACF) and partial ACF (PACF) diagrams of GCC data. (3) The parameters p, P, q, and Q range from 0 to 4, and the parameters d and D range from 0 to 1 and were constructed to fit the time series of the GCC, and then the lowest Akaike information criterion (AIC) value was selected considering the optimal SARIMAX model. (4) Forecast accuracy was evaluated using the root mean square error (RMSE) and mean absolute percent error (MAPE). We considered the lower values of these statistics to be better for forecasting accuracy. (5) The selected optimal model was diagnosed with the Ljung–Box test for the independence of residuals, and the Jarque–Bera test for the normal distribution of residuals. All time series analyses were performed using the “statsmodels” package in Python 3.6.3. [30].

3. Results

3.1. Time Series Data Decomposition

The summary statistics for the daily GCC values and meteorological factors, as well as the results of the Ljung–Box test and the ADF test, are described in Table 1. The mean GCC value was 0.38 ± 0.02 and was highly autocorrelated over time (p < 0.001) according to the results of the Ljung–Box test. The Ljung–Box test also showed significant autocorrelation in all meteorological factors. The ADF test showed that the temperature variables (Tmean, Tmax, Tmin) were statistically non-stationary.

To explore the time series data in more detail, the GCC time series was decomposed into three components: trend, seasonal, and residual (Figure 3). The trend component had the lowest value in early 2022 and gradually increased. The seasonal component showed a distinct seasonal periodicity of the GCC, typically peaking annually in summer.

This study calculated paired dynamic time warping (DTW) distances to compare the similarity of time series datasets (Figure 4). Although the observed GCC time series had higher similarities with temperature variables, the trend component, with seasonal and random variations removed to reveal the underlying long-term pattern, showed greater similarity to humidity and precipitation variables than temperature (Figure A1).

3.2. SARIMAX

To estimate the relationship between weekly GCC values and meteorological factors, the SARIMAX time series model was developed. The time series GCC data was tested for stationary by an augmented Dickey–Fuller (ADF) test. The test statistics confirmed the stationary after the first difference (−10.86, p < 0.001). The best-fitted SARIMAX model developed for the GCC was selected: SARIMA (0, 1, 1) (0, 1, 1)₅₂ with four meteorological exogenous variables (AIC = −65.8, BIC = −42.2) (

Table 2. The best fitted SARIMAX (0, 1, 1) (0, 1, 1)₅₂ model parameter’s estimate and their significance for the GCC on Mt. Sangwang.

Model Parameter	Estimate	SE	Prob > |Z|	Confidence Interval
				0.025	0.975
Tmax	0.095	0.056	0.091	−0.015	0.205
Hmax	0.039	0.024	0.096	−0.007	0.086
Solar radiation	0.041	0.027	0.135	−0.013	0.094
Precipitation	0.062	0.024	0.009	0.015	0.108
Moving average, lag 1	−0.160	0.052	0.002	−0.262	−0.058
Seasonal moving average, seasonal lag 1	−0.580	0.093	0.000	−0.763	−0.397
$\mathrm{Model}\mathrm{variance}\left({\sigma }^2\right)$	0.037	0.004	0.000	0.030	0.044
AIC = −65.8, BIC = −42.2, Ljung–Box (Prob (Q) = 0.95), Jarque–Bera (Prob (JB) = 0.07)

).

Our SARIMAX model includes the moving average term (lag 1), especially the seasonal moving average term (52-week moving average), which significantly accounted for seasonal variations of the GCC time series. To identify the importance of meteorological factors, we chose Tmax, Hmax, solar radiation, and precipitation as the optimal exogenous variables considering multicollinearity. Among the meteorological variables, Tmax had the highest positive feature importance, with a weight coefficient of 0.095, but it also had a larger standard error (SE) compared to other exogenous variables. Precipitation was the next important exogenous factor with a weight coefficient of 0.062 and a significant positive confidence interval. However, solar radiation and Hmax were not relevant to the prediction of the GCC within a 95% confidence interval.

Figure 5 displays observed and predicted values for visualization of the performance of our SARIMAX (0, 1, 1)(0, 1, 1)₅₂ model. This model was evaluated based on its RMSE and MAPE. When comparing the model performance with and without exogenous factors, the model with exogenous factors shows a slightly better fit to the observed data. The SARIMAX model with exogenous factor has the RMSE value of 0.246, suggesting that the model’s predictions are 0.246 units away from the actual values on average. Model performance was validated using an independent test dataset to ensure robustness. The plot demonstrates that the model’s predicted values closely align with the observed values, providing evidence of its accuracy.

To interpret the results of our model for the GCC, the diagnostic residual plots for the model were generated (Figure A2). Our model showed randomly scattered residuals, indicating that this model correctly captured the underlying patterns in the data. In addition, the normality of the residuals of our model was revealed by a histogram and a Q-Q plot. Finally, the correlogram of our model’s residuals shows no significant correlation at any lag, indicating the removal of any autocorrelation from the data. The developed SARIMAX model was tested for independence and normal distribution of the residuals using the Ljung–Box test (p = 0.95) and the Jarque–Bera test (p = 0.07), which showed the adequacy of the fit of the estimated model.

4. Discussion

This study explored the effects of meteorological factors on the GCC value of Q. acuta using time series analysis. Daily or weekly changes in the GCC and meteorological variables are recorded as a sequence of data points over an interval time, meaning that past values can influence future ones due to trends, seasonality, or cyclical patterns in the data [27]. These temporal dependencies should be considered in these time series datasets for the validity of the model [31]. Other researchers have emphasized the importance of time series analysis in the ecological field [15]. However, these methods are rarely applied to terrestrial ecology, particularly to ground-based phenocam datasets.

Most of the GCC value was focused on identifying the seasonal key events, such as the start of season (SOS), end of season (EOS), and peak of season (POS). The timing of these events is a vital mechanism for maintaining species coexistence, since plants are the base of the food chain [32]. Furthermore, it can determine how ecosystems respond to the impact of climate change [33], but it takes at least 10 or 20 years to identify the effect of climate on the tree’s phenological response because these seasonal patterns occur once a year. Additionally, the change of the GCC for the evergreen tree is vague, so detecting seasonal key events is more difficult in the evergreen tree compared to the deciduous tree.

To overcome these weaknesses, we applied a time series statistical method to identify the effect of meteorological factors on GCC variations of the evergreen tree. Time series statistical methods, such as ARIMA, SARIMA, SARIMAX, etc., are commonly applied to analyze or predict meteorological variables due to their data properties [34]. The satellite remote sensing dataset is also frequently explored by time series methods. For example, the ARIMA model was conducted to predict the vegetation health index and to estimate the relation with meteorological variables [35], and to evaluate the danger of forest fires for past and future climatic periods [36]. However, time series methods have not been applied to ground-based phenocam data. The Ljung–Box test and the ACF plots showed strong autocorrelation in the GCC index and meteorological factors, indicating that ordinary least squares (OLS) regression is not suitable, as it could result in unreliable statistical inferences and spurious relationships. However, numerous previous phenocam studies have conducted data analyses using ordinary least squares regression methods [37,38,39]. Our SARIMAX model demonstrated excellent predictions of GCC values and good diagnostic results using phenology data from September 2017 to August 2023. This study demonstrates that time series statistical methods are a valuable approach for analyzing ground-based phenocam data.

When comparing the long-term trends of GCC values with meteorological variables, the time series trends of humidity and precipitation showed greater similarity to GCC variations than those of temperature and solar radiation. Specifically, precipitation at our study site from October 2021 to February 2022 was only 55 mm, just 24% of the average rainfall. During the same period, the GCC values of Q. acuta were also significantly low. Our SARIMAX model also confirmed the relevance of precipitation for the GCC fluctuation. Precipitation was the only statistically significant exogenous variable, suggesting that variability in precipitation may have a greater impact on the phenological variability of Q. acuta than temperature. A previous study also suggested that variation of precipitation is a critical driver in phenological shifts in drylands [40]. Prolonged drought can affect groundwater availability, potentially triggering an unavoidable physiological response in trees or even leading to tree mortality [41]. The results of satellite-derived monitoring also show that the interannual variation in annual rainfall distribution influenced the length of the growing season [42]. Deep soil water can offer a relatively stable water source for trees, but the shallow soil causes water to drain quickly after rainfall at our site [43]. Consequently, prolonged droughts lasting several months can put significant stress on the trees. In particular, evergreen broad-leaved trees may be more sensitive to drought than deciduous broad-leaved trees during winter and early spring, as they continue to lose water through transpiration via their stomata [44,45], even in the drought season on the Korean Peninsula.

While our study provided statistical evidence of the impact of drought on the GCC value of Q. acuta, a limitation is that it did not establish a causal relationship between precipitation and GCC variation. To better understand the effects of meteorological variables on the phenological responses of trees, it is crucial to measure physiological indicators such as sap flow, photosynthetic rate, and osmotic potential. Our study site is situated at the boundary between the temperate and warm–temperate climate zones, where evergreen broad-leaved trees must adapt to low temperatures and drought during the winter season. At our site, Q. acuta may prioritize accumulating soluble carbohydrates, amino acids, and water-soluble proteins in their cells to regulate osmosis over enhancing photosynthesis. Investigating the combined effects of low temperatures and drought could be an important research topic. Additionally, monitoring soil moisture contents, which is directly influenced by precipitation, is crucial for estimating causal relationships.

5. Conclusions

In order to clarify the influence of meteorological variables on the phenological variation in the evergreen broad-leaved forests, we used ground-based digital photographic images and analyzed the collected data by means of time series methods. The results showed strong autocorrelation and seasonality in the GCC value collected at daily intervals. This means that unreliable statistical inferences and spurious relationships can be made using the ordinary least squares (OLS) regression method. It is crucial to carefully examine the characteristics of the data and to use appropriate statistical models to satisfy basic statistical assumptions. We recommend using time series analysis on phenocam data sets that are highly autocorrelated and seasonal. Our time series analysis highlighted the influence of dry season precipitation on the GCC variation of evergreen broad-leaved trees on the Korean Peninsula. Our results help to understand how plant phenology responds to climate change and improve the accuracy of phenological modeling.