A New Wetness Index to Evaluate the Soil Water Availability Influence on Gross Primary Production of European Forests

Abstract

Rising temperature, drought and more-frequent extreme climatic events have been predicted for the next decades in many regions around the globe. In this framework, soil water availability plays a pivotal role in affecting vegetation productivity, especially in arid or semi-arid environments. However, direct measurements of soil moisture are scarce, and modeling estimations are still subject to biases. Further investigation on the effect of soil moisture on plant productivity is required. This study aims at analyzing spatio-temporal variations of a modified temperature vegetation wetness index (mTVWI), a proxy of soil moisture, and evaluating its effect on gross primary production (GPP) in forests. The study was carried out in Europe on 19 representative tree species during the 2000–2010 time period. Results outline a north–south gradient of mTVWI with minimum values (low soil water availability) in Southern Europe and maximum values (high soil water availability) in Northeastern Europe. A low soil water availability negatively affected GPP from 20 to 80%, as a function of site location, tree species, and weather conditions. Such a wetness index improves our understanding of water stress impacts, which is crucial for predicting the response of forest carbon cycling to drought and aridity.

1. Introduction

Climate change results in rising air temperature and changes in the frequency and distribution of precipitation [1]. More-frequent extreme events are projected for the next decades in many regions around the globe [2]. Earlier studies have suggested a significant increase in frequency and severity of droughts [1,3], especially in the Mediterranean region [4,5,6,7,8,9]. Knowledge of mechanisms through which climate variations affect forest ecosystem processes, such as productivity, phenology, regeneration, and mortality, is still partial [10,11,12]. Consequently, there is an urgent need to realize a permanent assessment of forest sensitivity to climate change [13,14].

Mediterranean forests are considered particularly susceptible to future warming and drought conditions due to their location within the transition zone between arid and temperate climates [2,15,16,17]. Several studies have reported increased plant mortality rates and die-off events, reduced seedling recruitment, long-term shifts in vegetation composition, reduced radial growth, and increased crown defoliation responses [18,19,20,21,22]. Similarly, studies carried out in boreal forests, located in regions less constrained by drought, have reported reductions in productivity and widespread increases in tree mortality because of increasing drought stress [23]. On the one hand, warming conditions might promote drought stress, reducing gross primary production (GPP, representing the sum of gross carbon uptake by plant photosynthesis) and enhancing root and tree mortality [24,25,26,27]. On the other hand, warming conditions could increase: (i) forest growth from CO₂ fertilization, (ii) water use efficiency, and (iii) growing seasons duration [28,29,30,31,32].

Soil moisture plays a pivotal role in vegetation dynamics [8,33], considering that soil water availability is a crucial limiting factor for plant photosynthesis [34,35]. However, soil water content (SWC) is a critical variable depending on many factors such as precipitation amount, the duration and intensity of drought, level of soil erosion, surface runoff impacts, and rate of evapotranspiration [36,37]. Consequently, SWC is difficult to estimate both with field measurements and model estimations [38]. There are different methods to estimate SWC, such as volumetric and gravimetric methodologies [39], or remote sensing techniques [40,41]. The Soil Moisture and Ocean Salinity (SMOS) mission conducted by the European Space Agency (ESA) has the capacity to estimate soil moisture globally and routinely based on a microwave radiometer system [42]. To date, SMOS soil moisture data have provided particularly suitable long-term records of global soil moisture, even if limited to shallow soil [43].

Remote sensing can be used to generate a vast amount of information about the Earth’s surface and can be useful for suitable land-management planning. A strong relationship was found between SWC and values derived from surface temperature–vegetation index (Ts–VI) methods using optical and thermal remote sensing data as inputs [44]. Based on this method, a temperature–vegetation wetness index (TVWI) was developed to estimate SWC in a forest-dominated and topographically variable landscape in eastern Canada [45]. The integration of the normalized difference vegetation index (NDVI) and potential surface temperatures (θs) was used for an indirect assessment of soil moisture [45]. Following the concepts described above, a modified version of the temperature vegetation wetness index (mTVWI) was proposed using the complementary value of crown defoliation instead of NDVI [46]. In Europe, tree crown defoliation is monitored by the International Cooperative Programme on Assessment and Monitoring of Air Pollution Effects on Forests (ICP Forests) and is regarded as morphological damage affecting the green biomass of plants [47].

The driving hypothesis in this study is that a large variation of GPP, in response to soil water availability, depends on species, location, and weather characteristics, helping to better understand the relationship between GPP and mTVWI. Using the mTVWI, and considering GPP as a proxy for photosynthesis at the canopy level, we analyzed spatio-temporal variations of mTVWI across Europe, quantifying the impact of mTVWI and thus soil water availability on the GPP of European forests.

The analysis was carried out over the time period 2000–2010 for 19 representative tree species distributed across four bio-geographic areas using a MODIS GPP at 1 km of spatial resolution (Collection5 MOD17).

2. Materials and Methods

2.1. Data Collection

Crown defoliation data were provided by the ICP Forests network (Level I monitoring plots) for the main European tree species over the time period 2000–2010. In Europe, crown condition is the most widely applied indicator for forest health and vitality assessments [43,44,45,46,48,49,50,51]. Crown defoliation of the dominant tree species is visually estimated as the relative amount of foliage (expressed as percentage) in the tree crown compared to a reference tree of the same species with full foliage, growing in the same stand and site conditions [52,53]. Crown defoliation is estimated annually for 1623 monitoring plots of Level I according to the Manual on methods and criteria for harmonised sampling, assessment, monitoring, and analysis of the effects of air pollution on forests [53]. Mean annual values of crown defoliation were calculated for each Level I monitoring plot and for all considered species as reported in Figure S1 in the supplementary material. In particular, we considered nineteen tree species: six conifers (Larix decidua Mill., Picea abies L., Pinus halepensis Mill., Pinus nigra J.F. Arnold, Pinus pinaster Aiton, and Pinus sylvestris L.), twelve broad-leaf deciduous trees (Acer pseudoplatanus L., Alnus glutinosa Gaertn, Betula pendula Roth, Betula pubescens Ehrhart, Carpinus betulus L., Castanea sativa Miller, Fagus sylvatica L., Fraxinus excelsior L., Populus tremula L., Quercus petraea Liebl., Quercus pubescens Willd., and Quercus robur L.) and one broad-leaf evergreen species (Quercus ilex L.) distributed across four European bio-geographic areas (Boreal, Atlantic–Temperate, Continental–Temperate, Mediterranean). Furthermore, the crown defoliation database included site characteristics such as latitude, longitude, and altitude.

Mean annual temperature data (T mean, °C) were obtained from the European Climate Assessment & Dataset (ECA&D) for the time period 2000–2010. The ECA&D consists of a daily series of observations collected at several meteorological stations (>2000) throughout Europe. The nearest meteorological station to each ICP Forests site was used to provide temperature data.

GPP data (gC m⁻² day⁻¹) were provided by MOD17 product of the MODIS (MODerate Resolution Imaging Spectroradiometer) instrument aboard the Terra satellite, as modified by [54]. MOD17 is the first continuous satellite-driven dataset monitoring global vegetation productivity at 1 km of spatial resolution, providing 8-day, monthly, and annual GPP averages. The GPP component of the algorithm is estimated by a light use efficiency model and requires daily inputs of incoming photosynthetically active radiation (PAR) and climate variables [55,56,57]. The Collection5 MOD17 products, temperature, and crown defoliation data were used with a temporal resolution of 1 year, over the time period 2000–2010. All analyses were separately performed for all considered years in the 1623 plots.

2.2. Modified Temperature Vegetation Wetness Index (mTVWI)

The mTVWI [45,46] was defined as an indirect assessment of soil moisture related to NDVI and potential surface temperatures (θs). NDVI is one of the most extensively applied vegetation indices, closely associated to leaf area index (LAI), primary production [58], and the biochemical processes underlying primary production (photosynthetic activity). Following Vitale et al. (2014) [46], we used the complementary values of crown defoliation Γ, defined as 100 minus defoliation values as a proxy of NDVI.

Γ = 100 − Defoliation_plot

(1)

where Defoliation_plot is the mean crown defoliation value in one specific plot (in %). Γ was calculated for all considered plots and years. The potential surface temperature (θs) was calculated in two steps, as in [45]. The first step concerned the calculation of atmospheric pressure (p, in kPa), given by:

$$p=101.3 {\left[\frac{293-0.0065z}{293}\right]}^{5.26}$$

(2)

where z (in m) is the plot elevation relative to sea level. This equation is based on a simplified form of the ideal gas law for a neutrally-stratified atmosphere and a temperature of 20 °C at a standard atmosphere (i.e., 101.3 kPa).

Then, surface potential temperature θs (in K) was calculated as:

$$\mathsf{\theta }\mathrm{s}=T_s {\left[\frac{p_0}{p}\right]}^{R/c_p}$$

(3)

where T_S is the mean annual temperature (in K), p₀ is the average pressure at mean sea level (set at 100.0 kPa), R is the gas constant (i.e., 287 J kg⁻¹ K⁻¹), and c_p is the specific heat capacity of air (~1004 J kg⁻¹ K⁻¹). A scatter plot between θs and Γ was generated to derive two parameters, θ_dry and θ_wet (Figure 1), for each considered year. θ_dry (in K) indicates low levels of water availability limiting evapotranspiration; it was calculated as the straight line passing through points A and B, where A is the point with coordinates (max Γ; max θs in max Γ) and B is the point with coordinates (max Γ in the 99th percentile of θs; 99th percentile of θs) as depicted in Figure 1 and Figure S2.

θ_wet (in K) indicates an unlimited amount of water available to sustain evapotranspiration, and was calculated as the 1st percentile of θs. To derive the percentiles we considered all analyzed species together (1623 monitoring plots) for each year.

Finally, we calculated the mTVWI as follows:

$$\mathrm{mTVWI}= \frac{{\mathsf{\theta }}_{dry}-{\mathsf{\theta }}_s}{{\mathsf{\theta }}_{dry}-{\mathsf{\theta }}_{wet}}$$

(4)

where mTVWI is a dimensionless quantity varying from 0 (the case where water is not available for evapotranspiration, θ_dry) and 1 (the case where water is freely available for evapotranspiration, θ_wet). Consequently, the mTVWI is an index strictly linked to forest crown defoliation (θ_dry), surface potential temperature, and vapor pressure deficit (VPD) conditions [46].

All the above parameters were calculated on an annual basis during the time period 2000–2010.

2.3. Statistical Analyses

Under the assumption that partial correlation is a measure of association between two variables controlling or adjusting the effect of one or more additional variables [59], partial correlation coefficients define the degree of existing correlation between two variables when the effects of one or more variables are removed [59]. Since temperature strongly affects GPP and mTVWI, possibly causing confounding effects [60,61,62], we evaluated the partial correlation between annual data of GPP and the mTVWI after that temperature effect was removed [63] over the time period 2000–2010. We also performed a pair-wise Spearman’ rank correlation analysis using r_s coefficients [64], a non-parametric test to investigate both linear and non-linear correlations [59,60,65,66]. Significance was tested at p < 0.05 after Bonferroni’s correction for multiple comparisons [67]. One-way analysis of variance (ANOVA) was used to determine if there was at least one significant difference among the means of three or more independent groups. Specifically, we hypothesized that the means of mTVWI and GPP reduction in broad-leaf deciduous species significantly differed from mTVWI and GPP reduction in conifer and broad-leaf evergreen species. We used the Tukey post hoc test to find significant patterns and/or relationships (p < 0.05) between subgroups of broad-leaf deciduous, conifer, and broad-leaf evergreen species.

2.4. Impact of mTVWI on GPP

The impact of mTVWI on GPP was estimated similarly to the method proposed by [62,68] to quantify the impact of tropospheric ozone on GPP. The quantification of mTVWI influence on GPP (GPP related to water availability, GPP_Water) was expressed as follows:

$$GP{P}_{Water}=GPP\times mTVWI$$

(5)

where GPP is the original GPP derived from the MODIS product.

Finally, we calculated the percentage variation of the original GPP in respect to GPP_Water to evaluate the impacts of aridity on European tree species considered in this study for each year from 2000 to 2010.

2.5. Anomalies Analysis

The term anomaly is defined here as the difference between the value of a variable (e.g., GPP and mTVWI) for a single temporal period (year) and the grand mean of the same variable over a given time period (2000−2010). The analysis of anomalies was performed for each analyzed tree species, displaying GPP and mTVWI anomalies together in order to evaluate a recurrent pattern and to identify all deviations from the grand mean.

3. Results

3.1. Space–Time Distribution Patterns of mTVWI in Europe

The spatial distribution of the mTVWI values for the 19 plant species over the time period 2000–2010, is shown in Figure 2A–C. A north–south gradient in Europe is highlighted. Minimum mTVWI values (0–0.20) were mainly found in Southern Europe (Italy and Spain) while the maximum (0.80–1) were observed in Northeastern Europe (Belarus and Norway). In Central Europe, we found mTVWI values ranging from 0.40 to 0.80. Mediterranean evergreen oak (Q. ilex L.) showed a low mTVWI value (0.17) whilst B. pubescens, a typical Northern European species growing in wet conditions, showed a higher mTVWI value (0.87). In general, we found that broad-leaf deciduous species showed higher mTVWI values than conifers and the broad-leaf evergreen (Table S1). The analysis of variance (ANOVA) corroborated this evidence. In fact, the ANOVA and Tukey post hoc test indicated significant differences (p < 0.05) among mean mTVWI values referring to conifers, broad-leaf deciduous, and evergreen species for all years (Figure S3, Table S2). Taking into account temporal variations, Figure 3 shows the mean annual mTVWI values (±standard error, SE) over the time period 2000–2010 when considering all tree species together. We found the lowest mTVWI values in 2007 and 2008 (0.49 ± 0.01 and 0.33 ± 0.01, respectively) whilst in 2000, 2002, and 2009 we found mTVWI values ranging from 0.60 to 0.63.

3.2. Correlation between GPP and mTVWI

Partial correlation coefficients between GPP and mTVWI were significant (p < 0.05) in most of the European tree species (

Table 1. Partial correlation coefficients (controlling for: temperature (T)) for the correlation between gross primary production (GPP) and the modified temperature vegetation wetness index (mTVWI). Bold indicates significance at p < 0.05.

SPECIES	r (GPP, mTVWI)
	2000	2001	2002	2003	2004	2005	2006	2007	2008	2009	2010
Conifers	0.15	0.15	0.13	0.09	0.15	0.13	0.08	0.12	0.12	0.10	0.10
Larix decidua (N = 72)	0.27	0.34	0.35	0.19	0.40	0.38	0.40	0.38	0.33	0.39	0.35
Picea abies (N = 97)	−0.03	0.03	−0.04	−0.20	0.04	−0.10	−0.09	−0.02	−0.01	0.001	−0.11
Pinus halepensis (N = 82)	0.11	0.09	0.11	0.08	0.07	0.07	0.09	−0.02	0.05	0.11	0.12
Pinus nigra (N = 89)	0.26	0.25	0.23	0.17	0.18	0.26	0.15	0.16	0.17	0.21	0.15
Pinus pinaster (N = 87)	0.46	0.47	0.46	0.43	0.45	0.50	0.44	0.46	0.48	0.38	0.36
Pinus sylvestris (N = 77)	−0.39	−0.40	−0.38	−0.38	−0.36	−0.32	−0.39	−0.39	−0.34	−0.32	−0.29
Broad-leaf deciduous	−0.20	−0.20	−0.27	−0.25	−0.22	−0.21	−0.25	−0.22	−0.18	−0.10	−0.15
Acer pseudoplatanus (N = 100)	−0.30	−0.36	−0.33	−0.37	−0.28	−0.27	−0.36	−0.32	−0.31	−0.25	−0.25
Alnus glutinosa (N = 78)	−0.44	−0.44	−0.46	−0.43	−0.43	−0.34	−0.34	−0.43	−0.09	−0.09	−0.18
Betula pendula (N = 97)	0.45	0.56	0.46	0.46	0.61	0.42	0.24	0.22	0.09	0.12	−0.05
Betula pubescens(N = 71)	0.03	0.07	0.02	−0.004	−0.10	−0.21	−0.19	−0.01	−0.16	−0.003	−0.13
Carpinus betulus (N = 92)	−0.25	−0.36	−0.30	−0.41	−0.34	−0.32	−0.26	−0.29	−0.09	0.13	−0.15
Castanea sativa (N = 73)	0.01	−0.11	−0.04	−0.09	−0.08	−0.05	−0.07	−0.03	0.06	0.16	0.02
Fagus sylvatica (N = 82)	−0.14	−0.10	−0.12	−0.21	0.03	−0.01	−0.04	−0.15	−0.02	0.03	0.05
Fraxinus excelsior (N = 88)	−0.27	−0.30	−0.35	−0.28	−0.30	−0.28	−0.30	−0.34	−0.29	−0.24	−0.30
Populus tremula (N = 94)	−0.55	−0.54	−0.56	−0.42	−0.53	−0.48	−0.38	−0.69	−0.51	−0.25	−0.16
Quercus petraea (N = 85)	0.03	−0.04	−0.08	−0.13	−0.06	−0.10	−0.16	−0.06	−0.12	−0.03	−0.05
Quercus pubescens (N = 75)	−0.19	−0.29	−0.32	−0.36	−0.29	−0.18	−0.12	−0.20	−0.18	−0.03	−0.11
Quercus robur (N = 94)	−0.27	−0.34	−0.33	−0.35	−0.34	−0.29	−0.35	−0.32	−0.20	−0.09	−0.23
Broad-leaf evergreen	0.19	0.27	0.24	0.25	0.28	0.27	0.33	0.26	0.28	0.20	0.19
Quercus ilex (N = 90)	0.19	0.27	0.24	0.25	0.28	0.27	0.33	0.26	0.28	0.20	0.19

). In general, GPP and mTVWI correlated positively in conifers and broad-leaf evergreen species and negatively in broad-leaf deciduous species. We found positive correlation coefficients in L. decidua, P. pinaster, P. nigra, B. pendula, and Q. ilex ranging from 0.22 to 0.56 (p < 0.05). By contrast, P. sylvestris, A. pseudoplatanus, A. glutinosa, C. betulus, F. excelsior, P. tremula, Q. pubescens, and Q. robur showed negative coefficients varying from −0.24 to −0.56 (p < 0.05). Other species such as P. abies, P. halepensis, B. pubescens, C. sativa, F. sylvatica, and Q. petraea did not show significant relationships between GPP and mTVWI. Mediterranean species such as P. pinaster and Q. ilex showed the highest positive correlation coefficients between GPP and mTVWI (0.36–0.50 and 0.24–0.33, respectively). Conversely, species growing in wet conditions like P. sylvestris and P. tremula showed negative correlation coefficients between GPP and mTVWI.

Spearman rank-order correlation coefficients between GPP and mTVWI (

Table 2. Spearman rank-order correlation (rs) for the correlation between gross primary production (GPP) and the modified temperature vegetation wetness index (mTVWI). Bold indicates significance at p < 0.05.

SPECIES	rs (GPP, mTVWI)
	2000	2001	2002	2003	2004	2005	2006	2007	2008	2009	2010
Conifers	0.26	0.28	0.23	0.16	0.24	0.29	0.21	0.26	0.23	0.22	0.16
Larix decidua (N = 72)	0.12	0.17	0.16	−0.05	0.26	0.24	0.22	0.20	0.13	0.23	0.16
Picea abies (N = 97)	−0.19	−0.13	−0.17	−0.34	−0.06	−0.15	−0.16	−0.08	−0.08	−0.04	−0.22
Pinus halepensis (N = 82)	0.19	0.22	0.18	0.20	0.19	0.14	0.17	0.05	0.12	0.19	0.16
Pinus nigra (N = 89)	0.25	0.22	0.20	0.13	0.15	0.34	0.23	0.18	0.27	0.31	0.21
Pinus pinaster (N = 87)	0.47	0.45	0.45	0.46	0.47	0.47	0.48	0.46	0.51	0.40	0.42
Pinus sylvestris (N = 77)	−0.34	−0.38	−0.40	−0.40	−0.37	−0.36	−0.49	−0.28	−0.41	−0.48	−0.43
Broad-leaf deciduous	−0.58	−0.61	−0.66	−0.63	−0.65	−0.57	−0.63	−0.59	−0.53	−0.43	−0.44
Acer pseudoplatanus (N = 100)	−0.43	−0.54	−0.49	−0.52	−0.39	−0.44	−0.47	−0.38	−0.44	−0.39	−0.35
Alnus glutinosa (N = 78)	−0.60	−0.48	−0.67	−0.54	−0.63	−0.52	−0.59	−0.61	−0.40	−0.22	−0.21
Betula pendula (N = 97)	−0.39	−0.27	−0.54	−0.42	−0.45	−0.42	−0.48	−0.35	−0.26	−0.26	−0.08
Betula pubescens(N = 71)	−0.56	−0.55	−0.50	−0.50	−0.52	−0.48	−0.53	−0.60	−0.60	−0.44	−0.36
Carpinus betulus (N = 92)	−0.35	−0.49	−0.48	−0.40	−0.43	−0.26	−0.39	−0.30	−0.13	0.05	−0.20
Castanea sativa (N = 73)	−0.15	−0.32	−0.28	−0.29	−0.34	−0.20	−0.20	−0.19	−0.15	0.04	−0.22
Fagus sylvatica (N = 82)	−0.56	−0.59	−0.59	−0.73	−0.49	−0.47	−0.42	−0.56	−0.50	−0.30	−0.42
Fraxinus excelsior (N = 88)	−0.49	−0.52	−0.56	−0.53	−0.59	−0.47	−0.52	−0.55	−0.48	−0.46	−0.44
Populus tremula (N = 94)	−0.45	−0.34	−0.44	−0.29	−0.32	−0.43	−0.45	−0.31	−0.19	−0.18	−0.11
Quercus petraea (N = 85)	0.01	−0.22	−0.19	−0.21	−0.24	−0.19	−0.26	−0.09	−0.17	−0.08	−0.07
Quercus pubescens (N = 75)	−0.25	−0.31	−0.37	−0.43	−0.39	−0.25	−0.20	−0.20	−0.18	−0.07	−0.09
Quercus robur(N = 94)	−0.35	−0.51	−0.54	−0.53	−0.57	−0.47	−0.57	−0.38	−0.44	−0.33	−0.35
Broad-leaf evergreen	0.23	0.22	0.24	0.25	0.22	0.21	0.24	0.21	0.20	0.09	0.13
Quercus ilex (N = 90)	0.23	0.22	0.24	0.25	0.22	0.21	0.24	0.21	0.20	0.09	0.13

) were in agreement with partial correlation analysis results. The GPP and mTVWI relationships were positive and statistically significant (p < 0.05) in conifers and broad-leaf evergreen species and negative in broad-leaf deciduous species. We found positive correlation coefficients in L. decidua, P. nigra, P. pinaster, and Q. ilex (the same species highlighted in the partial correlation analysis except for B. pendula) and negative coefficients in most of the analyzed species (e.g., P. sylvestris, A. pseudoplatanus, B. pubescens, F. excelsior).

3.3. Impact of mTVWI on GPP

The percentage change of GPP due to mTVWI is shown in Figure 4A–C. The largest decrease occurred in Mediterranean countries, whereas the smallest decrease occurred in Belarus and Norway for all considered years. Italy and Spain showed a GPP reduction ranging from 40 to 100%. Central Europe presented a GPP reduction ranging from 20 to 40% for all years except for 2007 and 2008, where the GPP reduction reached 80%. Northern and Northeastern Europe were characterized by a small GPP reduction (0–40%).

Table 3. Species-specific GPP reduction (%) due to mTVWI over the time period 2000–2010 (±standard deviation).

SPECIES	GPP Reduction (%)
	2000	2001	2002	2003	2004	2005	2006	2007	2008	2009	2010
Conifers	−49.3 ± 21.88	−56.9 ± 25.21	−52.70 ± 23.57	−55.53 ± 24.95	−55.61 ± 25.27	−59.55 ± 25.31	−60.16 ± 26.51	−63.76 ± 20.52	−75.96 ± 15.19	−48.34 ± 19.79	−58.28 ± 20.57
Larix decidua (N = 72)	−47.03 ± 26.80	−52.55 ± 30.60	−49.39 ± 29.58	−52.17 ± 30.87	−51.48 ± 30.86	−54.92 ± 30.29	−53.94 ± 29.90	−61.32 ± 25.85	−73.36 ± 18.01	−44.48 ± 22.74	−54.16 ± 22.50
Picea abies (N = 97)	−46.81 ± 21.14	−53.01 ± 23.95	−50.89 ± 22.88	−52.73 ± 24.68	−52.40 ± 24.11	−57.16 ± 24.67	−57.40 ± 26.80	−62.71 ± 19.87	−74.80 ± 15.04	−48.39 ± 23.29	−55.68 ± 20.79
Pinus halepensis (N = 82)	−60.97 ± 13.87	−72.27 ± 15.21	−64.11 ± 15.04	−69.11 ± 13.14	−70.40 ± 14.33	−71.42 ± 15.97	−72.19 ± 17.94	−74.34 ± 11.54	−83.21 ± 9.55	−51.95 ± 13.39	−67.27 ± 16.19
Pinus nigra (N = 89)	−57.18 ± 19.34	−66.86 ± 21.72	−60.91 ± 20.76	−65.13 ± 21.87	−65.46 ± 22.80	−69.38 ± 23.43	−69.50 ± 25.00	−71.85 ± 18.04	−80.98 ± 14.90	−54.54 ± 19.29	−64.51 ± 20.11
Pinus pinaster (N = 87)	−49.86 ± 21.08	−58.87 ± 22.29	−54.76 ± 21.36	−57.23 ± 21.51	−56.95 ± 22.31	−62.92 ± 22.41	−65.48 ± 23.40	−62.71 ± 18.68	−77.25 ± 13.88	−51.19 ± 17.28	−62.38 ± 18.85
Pinus sylvestris (N = 77)	−32.52 ± 16.28	−36.38 ± 20.16	−34.10 ± 18.45	−34.71 ± 20.37	−34.96 ± 19.87	−39.02 ± 19.79	−39.88 ± 21.36	−47.88 ± 16.70	−64.88 ± 11.76	−37.46 ± 15.50	−44.0915.62
Broad-leaf deciduous	−30.13 ± 17.64	−33.47 ± 20.71	−31.57 ± 19.15	−32.25 ± 21.08	−31.87 ± 20.26	−35.97 ± 20.44	−36.09 ± 21.42	−43.21 ± 19.37	−60.62 ± 17.51	−34.71 ± 17.31	−41.54 ± 18.45
Acer pseudoplatanus (N = 100)	−39.80 ± 16.70	−43.04 ± 19.07	−41.37 ± 18.60	−42.05 ± 19.07	−40.33 ± 19.43	−43.62 ± 22.10	−45.78 ± 23.62	−50.00 ± 18.88	−67.49 ± 15.08	−41.71 ± 20.24	−46.72 ± 18.06
Alnus glutinosa (N = 78)	−23.37 ± 16.78	−25.45 ± 19.36	−23.83 ± 17.43	−22.99 ± 19.36	−23.79 ± 18.10	−28.17 ± 17.80	−27.56 ± 18.45	−37.90 ± 16.66	−58.46 ± 12.43	−30.27 ± 11.82	−38.93 ± 12.90
Betula pendula (N = 97)	−19.16 ± 14.85	−20.29 ± 16.48	−19.71 ± 14.01	−17.57 ± 17.24	−18.65 ± 16.06	−23.91 ± 15.97	−22.69 ± 15.05	−25.51 ± 21.38	−40.87 ± 28.07	−20.91 ± 15.69	−26.73 ± 19.35
Betula pubescens (N = 71)	−14.70 ± 8.83	−13.15 ± 8.81	−13.21 ± 8.49	−13.61 ± 8.86	−15.21 ± 9.29	−19.95 ± 10.49	−19.36 ± 10.92	−26.56 ± 9.98	−51.50 ± 8.50	−23.63 ± 9.12	−27.30 ± 23.74
Carpinus betulus (N = 92)	−34.32 ± 10.33	−35.58 ± 10.93	−32.82 ± 11.04	−34.92 ± 11.65	−33.04 ± 10.87	−37.06 ± 10.69	−40.06 ± 14.15	−46.35 ± 9.34	−64.33 ± 9.12	−41.71 ± 15.11	−48.02 ± 12.63
Castanea sativa (N = 73)	−37.85 ± 19.07	−44.29 ± 21.25	−41.53 ± 20.94	−43.46 ± 20.73	−42.68 ± 21.29	−47.42 ± 23.40	−48.30 ± 25.22	−52.09 ± 18.51	−67.89 ± 15.11	−44.45 ± 21.08	−50.65 ± 20.11
Fagus sylvatica (N = 82)	−50.00 ± 17.73	−55.39 ± 22.99	−53.18 ± 20.40	−55.21 ± 22.71	−53.60 ± 23.53	−58.35 ± 22.39	−58.94 ± 22.73	−64.32 ± 18.64	−75.85 ± 15.85	−49.01 ± 19.95	−58.94 ± 20.43
Fraxinus excelsior (N = 88)	−28.83 ± 15.02	−31.53 ± 16.76	−30.12 ± 16.06	−30.41 ± 16.85	−30.13 ± 16.09	−35.10 ± 17.82	−34.70 ± 17.48	−42.52 ± 13.44	−61.23 ± 11.33	−36.98 ± 15.12	−41.48 ± 13.84
Populus tremula (N = 94)	−16.82 ± 13.54	−18.97 ± 16.69	−18.38 ± 15.56	−16.13 ± 17.36	−17.38 ± 16.26	−21.96 ± 15.41	−20.94 ± 16.51	−30.15 ± 16.07	−52.00 ± 18.05	−26.42 ± 13.37	−34.31 ± 14.00
Quercus petraea (N = 85)	−30.58 ± 10.98	−34.87 ± 12.75	−32.67 ± 11.30	−34.48 ± 12.66	−33.18 ± 12.63	−36.39 ± 12.74	−36.91 ± 13.19	−45.67 ± 10.94	−62.12 ± 8.50	−33.61 ± 9.81	−41.08 ± 9.14
Quercus pubescens (N = 75)	−42.77 ± 13.46	−53.33 ± 16.52	−47.00 ± 15.40	−52.05 ± 16.14	−50.81 ± 16.49	−52.35 ± 16.88	−50.06 ± 18.09	−60.44 ± 14.63	−70.28 ± 10.76	−37.86 ± 11.42	−48.77 ± 13.50
Quercus robur (N = 94)	−25.50 ± 10.64	−29.02 ± 11.41	−27.54 ± 10.87	−28.07 ± 12.46	−27.62 ± 11.86	−31.21 ± 12.19	−31.02 ± 13.02	−41.24 ± 10.72	−59.25 ± 10.41	−31.63 ± 13.41	−37.82 ± 10.97
Broad-leaf evergreen	−57.63 ± 14.65	−71.26 ± 15.08	−64.20 ± 15.15	−68.88 ± 14.72	−69.85 ± 15.62	−73.12 ± 16.53	−72.16 ± 18.08	−74.46 ± 12.82	−83.41 ± 10.01	−54.08 ± 12.52	−67.98 ± 15.86
Quercus ilex (N = 90)	−57.63 ± 14.65	−71.26 ± 15.08	−64.20 ± 15.15	−68.88 ± 14.72	−69.85 ± 15.62	−73.12 ± 16.53	−72.16 ± 18.08	−74.46 ± 12.82	−83.41 ± 10.01	−54.08 ± 12.52	−67.98 ± 15.86

shows a plant species-based GPP reduction (%) for all years. Conifers showed a mean GPP reduction of nearly 50% ranging from 40% in P. sylvestris to 80% in P. halepensis, P. nigra, and P. pinaster. The broad-leaf deciduous species showed an average GPP reduction close to 40%, with some species showing higher reductions (54%, F. sylvatica and Q. pubescens) and other species with limited GPP reduction (20% B. pubescens, B. pendula, P. tremula). By contrast, Q. ilex presented values of GPP reduction close to 70% for all considered years. Figure 5, Figure 6, Figure 7, Figure 8 and Figure 9 show the average GPP reduction values over the time period 2000–2010 (standard deviations are shown in Figure S4). The largest decrease occurred in F. sylvatica, P. halepensis, P. nigra, P. pinaster, and Q. ilex, whereas the smallest decrease occurred in B. pendula and B. pubescens.

In general, we found that broad-leaf evergreen (Q. ilex) and conifer species showed higher GPP reduction values than broad-leaf deciduous (Table 3), which was also confirmed by analysis of variance (ANOVA) and Tukey post hoc tests (Figure S5, Table S3). The anomalies analysis is shown in Figure 10. We found a species-specific pattern of GPP and mTVWI anomalies in the time period 2000–2010. In general, highly negative anomalies of mTVWI (a condition of severe aridity) often correspond to a reduction of GPP. Indeed, in 2008 a highly negative deviation from the grand mean corresponds to a GPP reduction in most of the considered species (e.g., A. pseudoplatanus, F. sylvatica, P. abies, L. decidua, P. tremula). Taking into account the 2003 heat wave that affected much of Europe from June to September [69], we found a highly negative deviation from the grand mean of GPP not matching with high negative mTVWI anomalies.

4. Discussion

4.1. Spatial and Temporal Distribution of mTVWI

Our study assessed the impact of water availability in 19 tree species distributed across four European bio-geographic areas by using a modified temperature vegetation wetness index (mTVWI). The study focused on the impact of mTVWI, a proxy for soil water availability, on GPP. Our analyses revealed a north–south gradient of the mTVWI with minimum values in Southern Europe and maximum values in Northeastern Europe over the time period 2000–2010. Precipitation over the Northern Europe increased between 10 and 40% in the 20th century [70], whereas some parts of Southern Europe dried by up to 20% [70], reporting a geographic pattern consistent with our results. However, a north–south gradient of the mTVWI is not surprising as soil moisture and soil water availability affect vegetation growth and productivity especially in arid or semi-arid ecosystems [34]. In fact, the impact of water deficit on growth is much higher in arid and semi-arid sites in which water availability principally affects the main physiological processes of vegetation (growth, photosynthesis, carbon, and nitrogen use) than in sites with a moderate or well-balanced water supply [71,72,73,74]. A significant decrease of the mTVWI (−40% respect to others years) was observed in 2007 and 2008 over Europe. This reduction is likely related to high values of crown defoliation [75], especially in Mediterranean areas, impacting mTVWI values. Unfortunately, due to the lack of measurements of soil water content into the field conditions in our sites, it is unfeasible to compare the estimation of the index with real world measurements.

4.2. Correlations, Impacts, and Anomalies of GPP and mTVWI

Significant correlations between GPP and mTVWI were found in most of the European tree species. The highest positive coefficients were found for two Mediterranean species, Q. ilex and P. pinaster. Greater water availability allows higher stomatal conductance and thus a positive relationship between mTVWI and GPP [76]. Conversely, we found negative GPP–mTVWI relationships in other species (e.g., P. sylvestris, P. tremula) mainly growing in Central–Northern Europe. These results were well corroborated by Pasho et al. [77] in 2011, where they found a strong association between tree growth and a drought index (the standardized precipitation index) in P. sylvestris. On the other hand, other studies in Northern European countries did not find any association between SPI and crown defoliation (e.g., [78]). The different responses to drought from tree species are related to different strategies in coping with drought and local conditions of the forest stands [77]. As a function of site location and tree species, mTVWI negatively affected GPP from 20 to 80%. We found that broad-leaf evergreen (Q. ilex) and conifer species showed higher GPP reduction than broad-leaf deciduous species. A generalized increase in defoliation and mortality occurring was observed during 1987–2007 in Southern European forests (principally in Q. ilex stands), largely related to a severe drought regime [79], thus confirming our results. In agreement with our results, a high sensitivity to drought was found for broad-leaf Mediterranean forests [24]. P. pinaster, a Mediterranean species, showed high GPP reductions (50–80%) over the years, which are in line with earlier studies suggesting a high dependence of water availability in mixed pine–oak Mediterranean forests [22,80,81,82]. Probably, forests located in sites experiencing strong drought events (e.g., Southern European forests) could reduce their carbon sink efficiency and contribute to the reduction in carbon sink in the Northern Hemisphere [82,83]. However, the response of GPP to water availability significantly varied among species and locations. For example, B. pendula and B. pubescens, mainly distributed in northern temperate regions [84], showed low GPP reductions (10–30%). Conversely, F. excelsior was characterized by moderate-to-high GPP reductions (30–60%) being a water-demanding species suffering severe damages under increasing temperature [85]. Similarly, F. sylvatica showed high GPP reductions (up to 70%) likely due to its high sensitivity to drought [24]. Furthermore, an increase of crown defoliation was reported for the beech, especially in the Mediterranean area [50]. In fact, sensitivity to extreme temperature in F. sylvatica increased with the decreasing of soil water availability in the Mediterranean area [86].

Anomalies analysis showed a similar pattern of mTVWI and GPP anomalies over the time period 2000–2010, highlighting that negative anomalies of mTVWI (more aridity) induced a reduction of GPP. This result is in line with several studies (e.g., [87,88]) that reported a reduction of GPP NPP related to drought. Taking into account the 2003 heat wave, we found a highly negative deviation from the grand mean of GPP, whereas mTVWI did not show a corresponding variation. However, it should be noted that mTVWI was calculated considering temperature of the overall year (2003), whilst the peak temperatures in 2003 were found only during the growing season [68], thus impacting mainly GPP.

5. Conclusions

Models such as CMIP5 running under a higher temperature increased the potential evapotranspiration and changes in seasonal precipitation patterns and highlighted an intensification of drought events both in frequency and intensity [1,3], especially in the Mediterranean region [5,6]. Increasing climate and soil aridity are expected to cause growth decline and enhancement of mortality particularly in drought-sensitive forest species [15,16,89,90,91]. Thus, understanding and predicting the impacts of climate change on ecosystems is one of the pivotal challenges for the global change science [92,93]. As climate change is strictly dependent on the increase of the emissions of air pollutants, and consequently on air pollutant concentrations [1,3], we cannot study the impacts of climate change on soil aridity without considering air pollution impacts. Our study contributes to an assessment of GPP responses to soil water availability based on a dimensionless index—the modified temperature vegetation wetness index (mTVWI)—which is an indirect measure of soil moisture tightly linked to VPD conditions. Empirical results of this study outline a negative impact of low soil water availability on GPP, although it varied according to plant functional types, years, and bio-geographic areas. In fact, tree species can vary their ability to adjust physiological functions to soil water deficit [94,95]. Moreover, mTVWI could be used to estimate soil water availability when direct measurements and/or model estimations are lacking. Thus, its application can improve quantification of the carbon gain of forests undergoing drought promptly, to inform forest management planners of short-term climate change and its influence on the status of forest health.