Cleland et al. (2007) discussed limitations and opportunities of different methods to monitor plant phenology. The traditional method for detecting phenological phases relies on the utilization of observers (oftentimes volunteers) that register days that are associated with changes in phenological phases such as leaf unfolding during spring, bud formation, flowering, or leaf coloring and defoliation in late summer and fall. The phenological network provides conditions and thresholds to identify the various phases. One of the main problems arising when using data from phenological observations is the low spatial and temporal resolution. Phenological data can also be obtained using eddy-covariance flux measurements. As atmospheric carbon dioxide oscillates on an annual cycle, the amplitude and shift in this cycle can be used as an indicator of plant phenology. The collected data is however site-specific and in general, the setup for eddy-covariance monitoring is very expensive. Other ground and near-surface methods include digital camera sensing, spectral radiometers, and cameras carried by Unmanned Aerial Vehicles (Berra and Gaulton, 2021). The availability of reliable data is especially crucial for the evaluation and improvement of emerging, generally less complex technologies like phenology modeling and remote sensing. Phenology modeling allows not only for future estimates but also for the reconstruction of time series into the past where data is scarce. Limitations arise because the number of physiological processes described in the models is limited. Berra and Gaulton (2021) reviewed monitoring of forest phenology based on remote sensing which is oftentimes referred to as the Land-Surface-Phenology (LSP). Plant phenology shifts occur at a rate of 1 day per decade on average, however, the usual uncertainties from satellite LSP are much larger. This suggests that utilizing satellite data to track phenological changes still has considerable limitations.

The forest canopy extent can be described by numerous different parameters, with the most commonly used being the Leaf-Area-Index (LAI) which is defined as the projected area of leaves per unit of ground (Ross, 2012). LAI was shown to affect local temperature and vapor pressure deficit conditions, creating a characteristic forest microclimate (Von Arx et al., 2013). De Frenne et al. (2013) suggested that microclimatic buffering in temperature might help to decelerate a species turnover to the benefit of species from warmer regions (thermophilization) which is expected to occur due to macroclimate warming caused by climate change whilst mentioning that high closure in forest canopy might decrease the occurrence of light-demanding understory plants in forests. In a more recent approach, De Frenne et al. (2019) described the microclimatic buffering to be a potentially useful functioning to protect biodiversity. When deciding on a location and time for a field campaign, the LAI's dynamical nature must be taken into account: The main driver of spatial LAI variations between stands and years has previously been identified to be forest management, as thinning and harvesting drastically alter stand structure. Le Dantec et al. (2000). Aside from the high spatial heterogeneity, LAI is changing dramatically on various time scales. Besides year-to-year differences in the extent and natural seasonal variations especially in deciduous forests disturbances due to climatic or biotic factors can also have a major influence. Measurements of LAI are classified into direct and indirect methods. Breda (2003), Jonckheere et al. (2004), and Weiss et al. (2004) reviewed advantages and disadvantages of different measurement techniques in detail. Direct methods vary in their sampling destructiveness. Destructive methods can involve the removal of leaves from one or more model trees. Obtained LAI values are then up-scaled under the assumption of homogeneity within the stand. For a less destructive direct measurement littertraps can be placed and evaluated. The main disadvantage of direct methods is that they are in general time-consuming and labor-intensive and therefore only partially applicable for detailed long time monitoring. However, they are still needed for the calibration and validation of indirect methods. Indirect methods can be divided into contact and non-contact measurements. One commonly used indirect contact measurement is the utilization of allometric relationships between LAI and other parameters such as the sapwood area at breast height. These relationships however are not universally applicable since they are dependent on stand characteristics (Mencuccini and Grace, 1995; Le Dantec et al., 2000). An important group of non-contact indirect measurement methods are Beer's law-based optical methods such as hemispherical photographs. The setup for this method is usually a high-resolution digital camera equipped with a wide-angle fisheye lens mounted on a self-leveling system. The manually taken images are analyzed using software packages. Yan et al. (2019) described the three main issues of optical measurements to be the Leaf angle-distribution (LAD), clumping effects, and woody components. Moreover, Zheng and Moskal (2009) noted that a potential error source is the manual setting of exposure values if a camera is being used. Another issue when using photographs is that they only provide two-dimensional information. LiDAR measurements offer point cloud information, that could be used to retrieve three-dimensional data (Yan et al., 2019). It is distincted between three different LiDAR setups: terrestrial laser scanner (TLS), airborne laser scanner (ALS), and space-borne laser scanner (SLS). The latter of which is particularly relevant for generating LAI information on a global scale (Tang et al., 2014). The LAI can also be derived from satellite data utilizing vegetation reflectance properties. Typically, in that case, satellite-derived vegetation indices such as the normalized difference vegetation index (NDVI) (Rouse et al., 1974) or the enhanced vegetation index (EVI) (Huete et al., 2002) are being used to compute LAI estimates. Potithep et al. (2013) stated that seasonal changes are the most likely source of uncertainty in broadleaf forests and that it still requires in situ data of not just LAI but also vegetation indices to improve satellite-based LAI estimates. They further compared the determined LAI values to in situ measurements using NDVI and EVI data from MODIS images and discovered that EVI estimates produced better outcomes, particularly during the leaf growth phase. Non-green leaves cause mistakes in both VIs (from the saturation to the leaf fall period).

The objective of this study is to examine the modification of temperature inside a beech stand in comparison to open land conditions and analyze to which extent microclimatic differences can be used to quantify the status of the forest canopy. The presented method provides an LAI parametrization on a day-to-day basis which allows depicting seasonal changes. Furthermore, year-to-year differences in the occurrence and duration of phenological phases can be assessed. As stand data used in this study is collected as part of the ICP forests monitoring program similar time-series data are available in several European countries. Replicating the analysis might not only allow for detailed comparisons between different regions but moreover, results could be used as a reference when analyzing phenological phases determined using remote sensing technology. The novelty of this method is the generation of detailed phenological data on an operational basis.

This study focuses on a beech stand (Fagus sylvatica L.) in Klausen-Leopoldsdorf (Austria, 48°07′16″N 16°02′52″E, 510 msl) during the years 2011–2019. This site is a core plot of the Level II Monitoring initiative (de Vries et al., 2003; Neumann and Kindermann, 2016). Level II monitoring is an initiative of the ICP Forests program which started in 1995. In Austria 16 intensive observation areas were established. Six of which were selected as core plots, where additional parameters are being measured. Among others, the core plots supply data on meteorological conditions inside the forest stand. We used daily maximum temperature (T_F) inside the forest stand which was measured in agreement with the programme standards using a temperature sensor placed 2 meters above ground at measurement intervals of 15 min.

As a reference for open land conditions we used daily maximum temperature (T_R) provided by the Austrian National Weather Service (ZAMG). The station is part of the national weather monitoring service (TAWES) and located in open land also in Klausen-Leopoldsdorf in proximity to the Level II site (48°03′6.84″N 16°00′2.159″E, 389msl). Measurements are recorded in an interval of 10 min.

An overview of both stations and the research area can be found in Figure 1, maximum temperature characteristics for both stations can be found in Table 1. It is important to note that the lapse rate between the two stations has not been taken into account since it has no significant influence on the results and findings of our study.

For the determination of LAI values we used hemispherical photographs which were also conducted as part of the ICP Forests program. The photographs were taken in the near vicinity of stand temperature measurements. To get more representative LAI estimations, photographs were taken at 16 different photographic points in a 10 × 10 meter grid. We used a Canon EOS 500 D camera fixed on a 1.3 m tripod which was leveled and aligned with the magnetic North. We used the programme WinSCANOPY 2009a with the recommended settings to analyze images. Photographs were obtained once a year between the 3 and 29. August. The obtained LAI (LAI_hem) was on average at 1.96 with an overall standard deviation of 0.33. The maximum LAI_hem of 3.41 was found in 2012 where the standard deviation between the photographic points was with 0.41 also at maximum. The minimum value of 1.26 was found in 2011. The standard deviation between the photographic points was lowest in 2012 with a value of 0.08.

Phenological observations of leaf unfolding and leaf fall for Fagus sylvatica L. were also used for comparison. Since no data from fixed observation points was available we calculated the median day of observation in an area stretching in West-East direction from 14°45′57.9492″O to 17°9′38, 7″O and stretching in North-South direction from 48°30′35.334″N to 46°59′12.0084″N. Observations were provided by the ZAMG and obtained as part of the PhenoWatch programme (Templ et al., 2018). The point in time when at least three parts of the observed tree leaves have completely unfolded and in their final form (but not size) was classified as leaf unfolding. Leaf fall was defined as the point at which roughly half of the crown had defoliated.

For this analysis Python version 3.7 was used (Van Rossum and Drake, 2009) in addition with the packages numpy (Harris et al., 2020), pandas (Wes McKinney, 2010; pandas development team, 2020) and scipy (Virtanen et al., 2020) for data analysis as well as matplotlib (Hunter, 2007) and seaborn (Waskom, 2021) for visualization.

The median course of the year for DMAX can be found in the Supplementary Material (Supplementary Figure S1). The plot was sectioned into two different foliage phases: The "Maximum foliage phase" ranging from 8.June to 27.August (160–day 240) and the "No foliage phase" ranging from the 20.November to the 4.April (325–day 95). This process is relatively subjective and only gives approximate dates since it was based solely on the form of the curve. The median value of DMAX is at −2.55°C with a standard deviation of 0.89°C. The median during the no foliage phase in late autumn, winter and early spring is at −1.92°C with a standard deviation of 0.27°C. At the maximum foliage phase the median is at −3.86°C with a standard deviation of 0.33°C. The median course of the year for LAI_par can be found in Figure 2. Due to the definition of this parameter standard deviation values are the same as for DMAX, median values differ only in sign. The fit determined by applying the gompertz function to LAI_par values including measurement values can be seen in Figure 3, the function parameters can be found in Table 2. Characteristic days namely the day of maximum leaf growth (day 135) and maximum defoliation (day 289) were marked by vertical lines in the plot. LAI_par values and characteristic days for individual years were calculated and plotted. They can be found in the Supplementary Figure S2, Table S1 (for the years 2011–2015) respectively Supplementary Figure S3 (for the years 2016-2019).

During the measurement period the median value for DMAX_N is at −0, 17 with a variance of 0, 42. During the no foliage phase a median of −0, 35 with a variance of 1, 05 was obtained. It is important to note that high variation during this phase can partly be caused by maximum temperature approaching 0°C on some days. During the maximum foliage phase the median of DMAX_N remained relatively constant at –0,15 (with a variance of 0.0002). The constancy during the maximum foliage phase is especially of intereset, since formula 5 can be rearranged to:

Formula 7 emphasizes that DMAX_N describes the effective moderating effect the forest canopy has on forest maximum temperatures. When determining DMAX_N for the maximum foliage phase of individual years, minimal variation was found ranging from –0,17 (2016) to –0,12 (2017). Detailed averages and variances for the maximum foliage phase in each individual year can be found in Table 3. Because the value of DMAX_N during the maximum foliage phase showed minimal variation over the whole measurement period and variation was low between individual years, we tested if we could use formula 7 to approximate T_F using T_R and a constant DMAX_N of –0,15. A comparison of observed and approximated T_F values for the maximum foliage phase can be found in Figure 4.

A comparison between standardized LAI_hem and LAI_par values can be found in Figure 5. Boxplots were used to depict the variation of LAI_hem for the different photographic points. The comparison of annual values can be found in Table 4. To further investigate the influence of the year 2019 three scatter plots were created, showing the relationship between annual average LAI_hem and LAI_par for the timespans 2011–2017, 2011–2018, and 2011–2019, which can be found in Figure 6. Corresponding regression parameters as well as calculated Spearman correlations can be found in Table 5. Spearman correlation between LAI_hem and LAI_par was very strong obtaining a value of 0.93 with a p-value of 0.006 during the time span 2011–2017. A moderate correlation of 0.69 with a p-value of 0.021 was found during 2011–2018. During 2011-2019 correlation became negligible obtaining a value of 0.32 with a p-value of 0.757.

By analyzing the median course of the year for LAI_par (Figure 2) we were able to identify two constant foliage phases: The “Maximum foliage phase” during which the forest canopy is fully developed and no defoliation has yet taken place and the “No foliage phase” for which the constant low temperature difference can be interpreted as the blocking of solar radiation by the branches when defoliation is completed. Comparing the timing of the end of the no foliage phase with observed leaf unfolding showed a delay, with the no foliage phase ending on day 95 and the observed leaf unfolding being marked as day 108. However, this was expected as leaf unfolding was defined within the PhenoWatch program as the day where at least three parts of the crown leaves have already unfolded and have their final shapes and LAI_par will potentially also react to smaller disturbances such as bud formation and unrolling. For leaf fall a similar problem arises. As the observed leaf fall is at day 280 whereas the no foliage phase starts 45 days later. Per definition, this observation occurs when at least 50% of the crown is defoliated. However, the remaining crown is still blocking radiation and therefore affects DMAX, respectively, LAI_par. This could be the reason why leaf fall was observed before the beginning of the no foliage phase. The R² for the Gompertz function fit was at 0.68 during the leaf growth phase and at 0.63 for the defoliation phase. We compared the occurrence of the main phenological events as well as duration of phenological phases determined using the presented method to the results of previous studies on Fagus sylvatica L. foliage throughout the year:

A detailed comparison of characteristic days and duration of different foliage phases obtained via the four methods can be found in Table 6. Since the data sets come from very different climatic regions and from different years, this comparison is only used to check if characteristic days determined using the presented method are plausible. The determined day for the end of the no foliage phase using the LAI_par method is prior to the day determined using hemispherical photographs or radiation profiles and subsequent to the day determined using MODIS. Wang et al. (2005) argued that an early onset date for MODIS is found because it detects the onset of greens rather than a species specific event. It can not be ruled out that greening has an effect on the local temperature conditions and thereby also on LAI_par however the magnitude of this influence is unclear. The start of the maximum foliage phase was late in comparison to the day determined using radiation profiles and early in comparison to the day obtained using hemispheric photographs. The start of the no foliage phase was approximately the same as determined by MODIS. Considering the high variation of the determined phenological events between the different methods, the presented LAI_par method leads to plausible results. Bigger differences occurred when comparing the duration of specific phases. Where duration of the leaf growth phase was remarkably longer than duration obtained using the radiation model and duration of the maximum foliage phase remarkably shorter. Overall our findings suggest that LAI_par can be used to determine when considerable microclimatic changes are arising from foliage change within a deciduous forest.

We found the relationship between the observed and approximated T_F to be highly significant and concluded that DMAX_N can provide insight into stand-specific microclimate conditions and can be used to highlight the recreational function of forest land during summer. As an example, a hypothetical temperature prediction of 35°C for open land corresponds to a lower maximum temperature by 5, 25°C inside the stand. For individual years this reduction would range between 4.2°C (2017) and 5.95°C (2016).

Diminished daytime heating inside forest stands is of particular interest in regions where climate predictions suggest that the frequency of heat days and heatwaves is increasing. DMAX_N is presumed to vary depending on stand characteristics such as forest type, dominant species, altitude, and slope, which are known to have a considerable impact on temperature moderation inside forest stands (Renaud and Rebetez, 2009; Ferrez et al., 2011; Renaud et al., 2011).

Correlations were very strong during the time span 2011–2017. During the timespan 2011–2018 correlation was moderate. We identified a few possible reasons for the decline in correlation when including the year 2019 to our analysis:

Since LAI_hem measurements are only available annually and the period of the analysis is limited to seven years, caution should be exercised when interpreting the results of this section. On the one hand, the correlations during 2011–2017 and 2011–2018 seem to be promising, on the other hand, the drop in correlation when adding data from 2019 cannot be disregarded. To assess how well LAI_par is suitable as an LAI parameterization, the collection of more LAI data is needed, preferably also during different foliage phases. Our results suggest that LAI_par cannot replace conventional LAI methods but can provide interesting insights on stand foliage and microclimate conditions.