The intensive plots of the European Forest Monitoring Programme (ICP-Forests) provide high quality data on tree vitality and adaptability, nutrient cycling, critical loads, and water relations (de Vries et al., 2003). For this study we used data from four ICP-Forests core sites in Austria (Neumann and Kindermann, 2016), namely four Picea abies stands located between 715m and 1,540 m above sea level with a area of 2,500m² per site (Figure 1, Table 1). According to Köppen climate classification the sites are located in temperate and continental climatic regions. To give an overview of the site conditions climate diagrams for each sites are included in the (Supplementary Figure 1). The data used in this study were collected directly at the measurement sites in the period 2010-2020. Information on measurement sensors can be found in Supplementary Table 1. At an interval of 15 minutes open land parameters, namely temperature, relative humidity, wind speed, wind direction, solar radiation and precipitation were measured. Within the stand, ST and SM were measured at hourly intervals in 15cm depth. To ensure representativeness of the measurement, SM was measured in parallel at three different pits. The average value of these three pits is used for this analysis. At each site 10 trees are equipped with automatic band dendrometers manufactured by UMS München.

Dendrometers have a high measurement accuracy of ≤ 1 μm, which allows them to detect even small changes in tree girth, but this also means that the device reacts very sensitive to disturbances (e.g. litter, disturbances by animals) that can distort the measurement. Dendrometer data were cleaned using the R package treenetproc (Haeni et al., 2020; Knüsel et al., 2021; Wickham et al., 2019). Additional manual corrections were made where necessary, and in some cases years had to be removed completely because the amount of missing or corrupted data was too high. In total, we used 254 individual years from 40 tree individuals. The dendrometer value of day 91 (1. April) was subtracted from each individual year and individual tree in order to be able to better compare stem radius increments from individual years.

The target variable for modeling in this study was the stem radius (SR). To model SR we employed generalized additive models (GAMs) developed using the mgcv package framework in R (Wood, 2011; Wood et al., 2016). This framework allows to include smooth terms, which capture non-linear relationships between predictors and the response variable; parametric terms, which model linear effects; and random effects, which account for variability within hierarchical or grouped data structures. The smooth terms are composed by k simpler basis functions, which are summed to represent the non-linear function. Different types of basis functions are available such as thin plate regression splines, cubic regression splines and cyclic splines. The thin plate regression splines are the default splines in the mgcv package for one-dimensional smooth terms and the default number of basis functions is k=10. Both isotropic and anisotropic smooths can be used to include interactions in generalized additive models. The later, so called tensor interactions, can be used to model interactions between climate effects. Interaction terms use cubic regression splines as default. Moreover, the following random effects were included in all models:

Here TreeID is the individual and unique tree identification number, year is the current year and DOY is a variable containing the day of the year and the nearest hour (calculated by: DOY = Day of the year + Hour/24). The first random effect, s(TreeID, year, bs = “re”), accounts for variation between trees and years (nested structure where tree-specific responses are measured within each year). The second term, s(TreeID, year, DOY, bs = “re”), further captures any nested variation in tree responses within specific days of the year. In our modeling approach different climatic covariates were tested. As a significant amount of SM data was missing (3%), we imputed the data using also a GAM that contains the site ID and year as a random effect and the day of the year and precipitation (hourly cumulative sum, 12 hour cumulative sum, 24 hour cumulative sum, 48 hour cumulative sum). In addition to the parameters obtained directly at the ICP Forests sites, we included the VPD calculated from the measured parameters. We obtained ST and SM values at a depth of 15 cm as deeper layers might not respond as quickly to environmental changes.

As different meteorological parameters are measured at the sites, we tested different model configurations, each containing a parameter representing energy availability combined with a parameter describing moisture availability. Additionally we implemented models that contained both atmospheric and soil climate information to depict wether the combination yields better results. Furthermore, we adapted six models based on findings of recent studies from Etzold et al. (2022); Oogathoo et al. (2022); Rossi et al. (2006); Vospernik et al. (2020) and Zweifel et al. (2021) to our GAM framework and explored the results (Supplementary Table 2). Since soil water potential (SWP) data were not available at the measurement sites we used SM instead. The best fitting model was selected on the basis of AIC. The analysis of the model residuals indicated hourly autocorrelation and the final model was fit including an AR1 process.

To investigate the interaction between SM and ST, we created an artificial dataset representing various combinations of SM and ST values. These values were based on the observed 5%, 25%, 50%, 75%, and 95% percentiles for a specific day within the measurement period. To illustrate particularly extreme scenarios, we augmented the 95% percentile of ST by 1°C and 2°C to represent extremely warm conditions, and reduced the 5% percentile of SM by 2% and 4% to represent extremely dry conditions. These extreme values however remain within the range of historical climate. This artificial dataset was then used as input for the selected GAM to visualize how different combinations of SM and ST influence the prediction of SR. To showacse how these interactions might vary depending on the season we included visualizations for three different times during the growing season: the beginning of the growing season (May 25^th), the middle of the growing season (June, 21^st) and toward the end of the growing season (August, 10^th). The days were selected based on the average curves of stem radius. For the early growing season a day was selected where growth at all four stands has certainly started. The middle of the growing season was selected to be summer solstice. The selected day for the end of tree growing season was chosen when when stem radius reaches a stable value.

Given that the four measurement sites are spread throughout Austria with varying stand ages and growing conditions, in order to characterize the variations in annual radial increment between sites as well as between years, annual growth increment was extracted from the dendrometer timeseries (Table 2). On average, the annual growth increment across the four measurement sites ranged from 0.60 mm to 1.66 mm. Compared to the average site specific growth increment, deviations ranged from -50% (MZ in 2015) to +50% (JB in 2017). On average the first site to exceed 5% of annual growth was MR on April, 17th, followed by MS and JB on May, 2^nd and May, 3rd respectively, and latest for MZ on May, 11th. The first site to exceed 95% of annual growth was JB on July, 29th, followed by MR on August, 5th and MS on August, 6th, with MZ again being the last site on August, 17th (Figure 2).

Overall all models consisting of combinations of different energy and moisture limitating parameters yielded a good R² ranging from 0.97 to 0.98. R² of models based on literature varied between 0.85 and 0.98 (Table 3). The best performing model, P1, was selected based on AIC and consisted of the random effects TreeID, year and DOY (Section 2.2) and the following fixed effects:

and the tensor interaction:

The smooth terms in Equation 2 model non-linear relationships. s(DOY) captures the overall effect of the day of the year (DOY), while s(DOY, SM) and s(DOY, ST) allow for interactions between DOY and SM as well as between DOY and ST, respectively. They account for seasonal dynamics and how they interact with environmental conditions. The tensor product interaction ti(DOY, SM, ST) models the combined, non-linear effect of DOY, SM, and ST. The tensor interaction allows for flexible modeling of three-way interactions without assuming that the effects are strictly additive, enabling us to capture complex relationships between these variables. The variable DOY appears as both random effect (see Equation 1) and a fixed effect in all tested models to capture different aspects of its influence on tree growth. As a random effect, it accounts for the specific deviations in growth patterns for individual trees and years, allowing the model to capture tree-specific and year-specific variations in growth dynamics. As a fixed effect, DOY models the overall seasonal trend in growth across all trees and years. This formulation further helps in accurately capturing how the effects of SM and ST levels may vary depending on timing. A comparison of observed and predicted SR for the year and individual with the highest and lowest growth can be found in the (Supplementary Figures 2, 3).

In order to investigate the effects of different SM and ST levels on growth and how these effects may change seasonally we created conditional dependence plots for May 25th, June 21st, and August 10th (Figure 3). It is important to note that interpretation of the warm/dry conditions (shown in the graphs as either red shaded areas or dashed lines) should be done with caution, as there is little observational data available in this range. The overall distribution of SM and ST records can be found in Supplementary Figure 4.

At the beginning of the growing season SR increases with increasing soil temperatures up to 7-13°C, a further increase in ST leads to a decline in SR prediction (left graph in Figure 3A). An increase in SM yielded higher results in SR prediction, this increase was more pronounced when soil temperature was lower (right graph in Figure 3A).

During the middle of the growing season SR increased with increasing ST and SM respectively. Within lower temperature ranges an increase in SM has a bigger positive impact on SR in comparison to higher temperature ranges above 12°C. Strong positive effects of ST on SR are found under drier conditions whereas predictions under conditions with high moisture vary only minor in comparison with changes in ST (Figure 3B).

Toward the end of the growing season increasing ST has a constant or positive effect on growth however a strong variation with SM can be found. In general SM has a strong positive effect on growth which is more pronounced in the lower temperature ranges (Figure 3C).

A more detailed look into how variation of the two parameters (SM and ST) can be given when keeping one of the two constant (at the median value for the given day) and assessing the change in growth prediction when the other varies between the 5% and 95% percentile. Starting with constant SM, early in the growing season growth decreased by 0.03mm when ST increased, in the middle of the growing season it increases by 0.08mm and later in the growing season it increases by 0.08mm. If ST is fixed and SM increased, early in the growing season growth increased by 0.08mm, in the middle of the growing season it increases by 0.17mm and toward the end of the growing season it increases by 0.17mm. Hence both an increase in ST and an increase in SM has the biggest overall positive effect in the middle and toward the end of the growing season. When comparing the effect of variation of either ST or SM on the SR (-0.03mm vs. 0.08mm early in the growing season, 0.08mm vs. 0.17mm in the middle of the growing season, 0.08mm vs. 0.17mm toward the end of the growing season) it can be seen that SM variation has a stronger impact on the growth prediction than ST variation.

The derived model effectively captures the overall annual increments, including annual patterns; however, its ability to represent swelling and shrinking remains limited (Supplementary Figures 2, 3). We explored several approaches to better capture the processes of swelling and shrinking, including the incorporation of temperature moving averages and lagged precipitation sums into the model. However, these modifications did not lead to significant improvements in the model's ability to represent these dynamics. One possible explanation for this is that, while SM and ST are good predictors of SR increment, other environmental factors cloud be more closely related to the processes driving shrinking and swelling. Another possible explanation is the substantial variation in the magnitude of shrinking and swelling observed across the four research sites. This variability may cause the model to adopt a more conservative representation of these processes, particularly at sites where they are more pronounced. When evaluating the average RMSE throughout the course of the year, we observed seasonal changes. RMSE fluctuates around 0.09 mm at the start of the year, rises around the beginning of July (fluctuating around 0.11 mm), declines in the middle of August (down to 0.04 mm), and then increases once more at the start of September (up to 0.15 mm) (Supplementary Figure 5).

Previous studies reported that atmospheric conditions had a greater influence on tree growth than soil climate conditions (Oberhuber et al., 2014; Lupi et al., 2011). In contrast, our model, which incorporates both SM and ST, outperformed those relying solely on atmospheric parameters, as indicated by AIC. ST can be highly correlated with air temperature (which is used to calculate VPD), though the variability in soil temperature is typically dampened. This raises the possibility that the effects we observe might also be influenced indirectly by VPD. However, the correlation between ST and air temperature can vary considerably depending on site-specific factors and measurement depth of the soil temperature sensor (Zhang et al., 2022; Islam et al., 2015). At our study sites, we found a moderate Pearson correlation of 0.6 between ST and air temperature. Additionally, the correlation between VPD and ST was only weak (0.2), suggesting that the observed effects of ST may not primarily operate through VPD.

Furthermore, we observed that the effect of increasing ST on SR shifted early in the growing season from being beneficial to being detrimental (left graph in Figure 3A). This is in contrast to previous findings on effects of elevated air temperature (Deslauriers et al., 2008; Lupi et al., 2010). In principle, higher temperatures can enhance metabolic activity and radial stem growth (Mäkinen et al., 2001, 2003). Especially at the beginning of the growing season when SM levels are high due to snow melt. However, our findings suggest that other factors may modulate this response when ST gets to high. One possible explanation is that trees allocate more resources to root development at the expense of radial stem growth during this period. This hypothesis aligns with Lahti et al. (2005) who found that diameter growth in five-year-old Picea abies peaked at approximately 9°C and declined with higher ST.

Later in the season, increasing ST appears to have mostly positive effects on growth, as indicated by the model's higher predictions under moist conditions (left graph in Figures 3B, C). This is consistent with studies that report temperature-induced enhancement of photosynthesis and growth when soil moisture is adequate. Overall, our findings highlight the complexity of the relationship between ST and SR and suggest that local environmental factors, such as soil properties, seasonal dynamics, or tree physiological strategies, may play a significant role. While further research is required to fully elucidate the mechanisms behind these observations, our results underline the need to understand SR responses to temperature variations. Overall, our findings highlight the complexity of the relationship between ST and SR and suggest that local environmental factors, such as soil properties, seasonal dynamics, or tree physiological strategies, may play a significant role. While further research is required to fully explain the mechanisms behind these observations.

Growth declines observed during the study period correspond to years identified as dry in Europe, even though the yearly climate diagrams for the study sites do not indicate evidence of drought during these periods. This suggests that drought effects at the study sites are visible in the form of growth declines, but a finer temporal resolution is needed to fully capture these dynamics. This supports the hypothesis of Zweifel et al. (2021) that drought should be assessed at the seasonal and daily level and is in line with Oberleitner et al. (2022) who report drought limitation also for sub-alpine sites. The observed growth decreases are however smaller than on more xeric sites (Table 2).

When studying sites located at lower altitudes (between 319m and 425m) D'Andrea et al. (2023) found that Picea abies is affected more by precipitation (and therefore water availability) than temperature. This is in accordance to our findings as it was also evident that SM had a bigger impact on SR than ST. However Alavi (2002) reported that even in areas with high annual precipitation, water availability significantly affects tree development, with drier sites experiencing greater growth variation. Our study is consistent with these observations, as none of our sites were water limited, yet we still observed significant SR variation associated with changing SM conditions (Figure 3). Furthermore, Ježík et al. (2015) showed that non-irrigated trees had significantly lower growth compared to irrigated trees, highlighting the positive influence of adequate SM on growth. Our results echo these findings, as increased SM generally led to increased predictions in our model.

It is important to acknowledge constraints given by the selected modeling approach. Among these is the model's exclusion of legacy effects. Previous research found that especially precipitation from the previous year had a strong impact on the growth of the radial stem (Miller et al., 2023; Zweifel et al., 2020). Moreover age effects on growth are only taken into account indirectly by the random effects implemented within the model. The higher average growth for MS, seen in Figure 2, is likely caused by such effects as the stand age of 26 is lower than for the other three stands (ranging from 76 and 133 years).

Despite these limitations, we observed a strong relationship between SR of Picea abies and both SM and ST. This relationship may be stronger compared to other tree species, as previous research has shown that Picea abies tends to react more sensitively to environmental conditions (Vacek et al., 2019). Our findings underscore the importance of incorporating SM and ST measurements into future study designs.