Assessment of inter-annual forest production variations in Italy by the use of remote-sensing and ancillary data

ABSTRACT The application of Monteith’s models to estimate the gross primary production (GPP) of forest ecosystems is a common and efficient practice. A modified version of one of these parametric models, C-Fix, has been widely tested in Italy at daily to annual temporal scales. The current paper is a step forward demonstrating the capability of this model to reproduce also inter-annual forest GPP variations. To this aim, dendrochronological measurements collected in two study areas were utilized to validate the model. The first study area is a coastal Mediterranean pine forest and the second is a more temperate ecosystem dominated by the presence of European beech. The experiment was organized into the following main steps. Modified C-Fix was applied to simulate annual GPP over a 15-year period for both sites. Both the main driving factors used by the model and the estimated GPP were then analysed and assessed against the detrended tree-ring width series of the collected cores. The study yields indications on the factors most influential on GPP for the two forest ecosystems and supports the capability of the proposed modelling approach to reproduce relevant inter-annual GPP variations, which is useful for predicting the impacts of ongoing climatic changes.


Introduction
The estimation of forest production over large areas is important for both scientific and practical purposes. Such an assessment is, in fact, necessary to study the global carbon cycle and produce information useful for planning and managing forest resources (Waring & Running, 2007). Additionally, to understand the effects of the ongoing climatic changes, it is important to detect and follow changes of forest gross primary production (GPP) occurring over time.
Conventionally, forest GPP is estimated by using the eddy-covariance technique as proposed by Aubinet, Vesala and Papale (2012). Regional-scale applications of this technique, however, are expensive and not representative for large areas due to the high spatial variability of the factors which affect forest production, such as, for example, climate, topography, soil fertility and management practices. Moreover, flux tower measurements started only in the 1990s and are difficult to obtain over long time periods due to technical and budget limitations (Baldocchi, 2003). Thus, a procedure mostly independent of field measurements and driven by remote-sensing data would be extremely useful to yield spatial estimates of forest production (Running et al., 2004).
Remote-sensing techniques, in fact, provide information on forest conditions at different spatial (from local to global) and temporal scales (from daily to annual) (Waring & Running, 2007). In particular, the remote-sensing estimation of forest GPP is usually based on Monteith's approach (Monteith, 1972), which exploits the direct linkage between remotely sensed vegetation indices and the fraction of photosynthetically active radiation (fAPAR). Such approach calculates GPP as the product of absorbed photosynthetically active radiation (PAR) absorbed by vegetation and radiation use efficiency (ε), typical for each forest type.
This approach is the basis for many widely utilized parametric models, such as CASA (Field, Randerson, & Malmstrom, 1995), 3-PGS (Coops, Waring, & Landsberg, 1998), MODIS GPP/NPP (Heinsch et al., 2006;Running et al., 2004) and C-Fix (Veroustraete, Sabbe, & Eerens, 2002). These models mostly differ in how to account for thermal and water limitations to potential production, namely in the modulation of the radiation use efficiency (Sanchez-Ruiz et al., 2016). In particular, water stress can be estimated by the consideration of vapour pressure deficit (VPD), i.e. in the MODIS algorithm, or soil moisture (Potter et al., 2003). The use of VPD, however, presents some theoretical and practical drawbacks in arid and semi-arid environments, where it underestimates the effect of water stress (Mu, Zhao, & Running, 2011). On the other hand, soil moisture is not easily estimable over large areas, at proper temporal scales (Zhang & Zhou, 2016).
C-Fix was originally applied in some European countries (Veroustraete et al., 2002) and then adapted to the Mediterranean environment by Maselli, Papale, Puletti, Chirici, & Corona (2009). These authors introduced an additional factor, the coefficient of water stress (Cws), which accounts for short-term water stress and can be easily obtained from a simplified water balance using conventional meteorological data. The modified model has been widely evaluated in terms of daily values and total annual amounts against eddy covariance tower data (e.g. Chirici et al., 2015;Maselli, Barbati, Chiesi, Chirici, & Corona, 2006). The results were promising both in temperate and Mediterranean areas, but these investigations could not assess the model performances in terms of inter-annual GPP variations, due to the mentioned lack of long-term reference observations (see the Fluxnet website, http://www. europe-fluxdata.eu/). This situation calls for a careful testing of the model performances against multi-year production observations. The current work was therefore aimed at assessing the capability of modified C-Fix to follow the inter-annual GPP variations of two forest areas in Italy, representative of different environmental conditions and functional forest types. The first is a Mediterranean semi-arid area mostly covered by an evergreen coniferous species, the umbrella pine (Pinus pinea L.). The second is a mountain temperate-cold area, dominated by a deciduous broadleaved species, the European beech (Fagus sylvatica L.). Different ecological factors constrain the growth of the two forest ecosystems: in the Mediterranean coastal pinewood, summer water availability at root level is the main limiting factor (Teobaldelli, Mencuccini, & Piussi, 2004), whereas temperature and radiation mainly drive the development of the mountain beech forest (Lombardi et al., 2016). The analysis was therefore extended to consider the main drivers utilized by the model, i.e. meteorological factors and fAPAR. The assessment of these variables and of the final GPP estimates was performed against annual tree-ring widths collected in the same areas. This assessment relies on assuming an approximately constant ratio between GPP and net primary production (NPP), which is the variable more directly related to tree-ring widths (Bascietto, Cherubini, & Scarascia-Mugnozza, 2004). While this assumption is unrealistic for single trees in specific growing phases, it can hold for relatively large forest areas and long time scales (Waring, Landsberg, & Williams, 1998).
The current paper is organized as follows. The main features of the two Italian forest sites are described together with those of the ground and remotely sensed data utilized. Next, the basis of the C-Fix methodology and its application in the two areas are presented, together with the analysis of inter-annual GPP variations. This is followed by a description of the results obtained and by a final section where the possibilities and limitations of the proposed approach are discussed and commented.

Maremma and Molise
Two study areas characterized by different environmental and vegetation features ( Figure 1; Table 1), here referred as Maremma and Molise, were selected in Central Italy. The first area, Maremma, is situated in a coastal area in the southern part of Tuscany. The climate is typically Mediterranean, with mean annual temperature of about 15.5°C and limited rainfall (around 650 mm), distributed from mid-autumn to early spring (Mazza & Manetti, 2013). The area, covered by a Mediterranean pine (umbrella pine, P. pinea L.) forest, is characterized by sandy soils and a high salty water table level (ranging between 1 and 2 m) (Pranzini, 1996). During the period of short fresh-water supply (from mid-summer onwards, depending on annual meteorology), the pines might have to cope with the underlying salty water, which can determine a clear reduction of tree growth (Teobaldelli et al., 2004). The intensity of this phenomenon is highest near the shoreline and decreases going towards the inner hills. The second study area, Molise, is located in the central-southern Appenines. The sampled sites can be grouped into three clusters following the main mountain groups that are called Mainarde, Matese and Alto Molise (Lombardi et al., 2016). The climate of the area is temperate, with mean annual temperature of about 14°C and rainfall of about 1600 mm. All forests are mostly situated on high elevations, along the Appenine chain, and are characterized by the presence of European beech.

Meteorological data
For the Maremma study area, daily meteorological data (i.e. minimum and maximum temperature and rainfall) were derived from the nearby meteorological station of Alberese. These data were collected over the years 1986-2000.
For the larger area of Molise, data were derived from the Pan-European E-OBS dataset downscaled at 1-km spatial resolution (Fibbi et al., 2016;Maselli et al., 2012). This meteorological dataset consisted of complete daily series of temperatures (both minimum and maximum) and precipitation collected over the period 1999-2013.
For both study sites, daily solar radiation was calculated applying the MT-CLIM algorithm (Thornton, Hasenauer, & White, 2000) and converted into PAR by a multiplying factor (0.464) (Iqbal, 1983).

Satellite data
For the Maremma study site, National Atmospheric and Oceanic Administration Advanced Very High Resolution Radiometer (NOAA-AVHRR) NDVI data for the years from 1986 to 2000 were derived from the archive of the University of Berlin. As fully explained in Koslowsky, Billing and Eckard (2001), the original band images, having a 1-km pixel size, were acquired by several NOAA missions and were not corrected for the atmospheric effect. The final product consisted of thirty-six 10-day NDVI maximum value composite images for each of the 15 available years (Holben, 1986).
Concerning the Molise study area, Satellite Pour l'Observation de la Terre-VEGETATION (Spot-VGT) NDVI images were considered for the period 1999-2013. These data, which have spatial and temporal resolutions similar to those taken by the NOAA-AVHRR, were freely downloaded from the VITO archive (http://free.vgt.vito.be) in a preprocessed format, which comprises both atmospheric and geometric corrections (Maisongrande, Duchemin, & Dedieu, 2004).

Dendrochronological data
Tree cores were taken during field campaigns performed in the two study areas. More specifically, the ground samples were collected in 2003 and 2014 in Maremma and Molise, respectively (Lombardi et al., 2016;Teobaldelli, Cherubini, & Piussi, 2008). Concerning the first study area, only the samples most distant from the shoreline were currently considered to minimize the effect of the salty water table. All tree-ring width data series were collected and preprocessed following standard methods, useful to build a tree averaged series and the mean chronologies of the sampled sites (Fritts, 1976). Thus, linear annual increments were finally obtained over the two 15-year periods for both study areas.

Data processing
Computation of forest GPP C-Fix is a Monteith type model driven by temperature, radiation and the fAPAR (Veroustraete et al., where ε is the maximum radiation use efficiency (1.2 g C MJ −1 ) (Chirici et al., 2015), Tcor is the MODIS temperature correction factor and Cws is the water stress factor, all referred to day i. While fAPAR can be derived from NDVI, the two correction factors are computed from daily meteorological data. Tcor, which takes into account the effect of temperature on photosynthesis, is calculated as a linear ramp function of daily minimum air temperature following Running and Zhao (2015). This scalar is adopted also in the MODIS GPP algorithm and differentiates between biome types (i.e. evergreen conifers and deciduous broadleaves in the case of Maremma and Molise, respectively).
Cws acts reducing photosynthesis in case of shortterm water shortage, while long-lasting dry periods are captured by the NDVI profile (Maselli et al., 2009). This scalar is computed from a simplified water balance as where AET and PET are actual and potential evapotranspiration, respectively, both cumulated over 2 months for forest ecosystems. AET is assumed to equal precipitation when this is lower than PET and to equal PET in the other cases; Cws can, therefore, vary between 0.5 and 1, meaning that GPP can be reduced up to a half of its potential value. Equation (1) was applied to both Maremma and Molise forest sites using the meteorological and NDVI data previously described. Daily PET was computed following Jensen and Heise (1963). The NDVI values of the two study forests were extracted from the available NOAA-AVHRR and the Spot-VEGETATION archives for the two study periods. Additionally, NOAA-AVHRR NDVI values, which are affected by several disturbances (Koslowsky et al., 2001), were converted into radiometrically and atmospherically corrected NDVI values by applying the linear regression equation reported in Maselli (2004) (Figure 3b). For both sites, the corrected NDVI values were linearly interpolated on a daily basis and converted into fAPAR using the generalized linear equation proposed by Myneni and Williams (1994). Daily GPP was finally obtained for the two study areas over the 15-year periods considered..

Trend and correlation analyses
The collected meteorological and NDVI data were first analysed to characterize the evolution of the driving factors during the respective study periods. To this aim, linear trend analyses were applied to the annual series of the most influential model drivers. In the semi-arid area of Maremma, attention was focused on the drivers which account for water limitation, i.e. NDVI and the water stress factor. The first was aggregated on an annual basis, while Cws was averaged over the period from March to July, which is the most influential on inter-annual GPP variations of Mediterranean forests (Maselli et al., 2014a). In the mountainous areas of Molise, attention was focused on the model drivers accounting for thermal limitation, i.e. NDVI and the thermal scalar. Due to the deciduous nature of these forests, NDVI and PAR were averaged over the growing season, from May to October, while Tcor was averaged over the period January-April, which is most influential in temperate-humid zones Kim, Kimball, Didan, & Henebry, 2014). PAR was also analysed for both sites, averaged on an annual basis for the evergreen forest of Maremma and over the growing season (May-October) for the deciduous beech ecosystems of Molise.
All averaged driving factors (i.e. Tcor, Cws, NDVI and PAR), as well as annual GPP, were finally assessed against the available tree-ring widths. Treering width chronologies, in fact, although possibly biased (Babst et al., 2014), may be used as indicators of the woody biomass accumulated by trees during different growing seasons (Bascietto et al., 2004). Under the previously mentioned assumption of constant ratio between NPP and GPP, these chronologies can, therefore, be subjected to correlation analyses aimed at evaluating both the effect of the main model drivers and the reliability of the final estimates.
Tree-ring width series, however, are affected by age trends and stand disturbances, such as management operations, whose effects can be reduced by the application of a detrending operation (Cook & Peters, 1997). This consists of transforming the annual treering widths series as where Y detr is the detrended value of tree-ring width, Y is the original value, while Y and σ Y are the moving average and the standard deviation of order three, respectively (Rodolfi, Chiesi, Tagliaferri, Cherubini, & Maselli, 2007). This operation does not alter the major statistical properties of the original data and, in particular, it does not introduce artificial correlations between independent data series (see e.g. Cook & Peters, 1997). The use of detrended tree-ring width series for correlation analysis, however, requires a similar detrending of all other variables to be correlated (i.e. meteorological factors, NDVI and GPP estimates), because the presence of trends in these variables can confuse the possible inter-relationships (Rodolfi et al., 2007). For this reason, all data series to be used in the correlation analyses were previously detrended by the application of Equation (3).

Case study of Maremma
The mean monthly values of precipitation, PET and NDVI obtained by averaging the data over the years 1986-2000 are shown in Figure 2. During this period, the annual precipitation is about 580 (±150) mm year-−1 , distributed with two maxima in spring and autumn and a minimum in the summer season. On the opposite, PET reaches the highest values during summer, following temperature and solar evolution. Its annual amount is about double of precipitation (about 1000 mm year −1 ). Water shortage periods are, therefore, noticeable and prolonged from March-April until September; this determines dry conditions for vegetation, which is evidenced in the corresponding NDVI evolution ( Figure 2). In particular, NDVI shows a reduction during summer but, due to the evergreen nature of the considered species (i.e. umbrella pine trees), its annual range is very limited. Figure 3a shows the evolution of the main drivers which contribute to determine the annual GPP based on the above described Monteith's approach. The water stress factor, Cws, computed over the months from March to July, decreases following precipitation during the 15 years; it also shows high variability, with minima in 1987, 1990 and 1997, which partly contribute to determine the NDVI evolution. The same figure highlights the variation of mean annual NDVI data, which partly follows the same pattern of Cws, decreasing from about 0.6 in 1986 to about 0.45 at the end of the period. Figure 3b reports the evolution of annual PAR. This driver also decreases during the observed period contributing, with the two previously described factors (i.e. NDVI and Cws), to determine a reduction of annual GPP. More specifically, the latter decreases from about 1100 g C m −2 year −1 in 1986 to about 800 g C m −2 year −1 in 2000.
The influence of all driving factors in determining the inter-annual variations can be quantified by Table 2, where the correlation coefficients between detrended drivers and tree-ring widths are given. As expected, the highest influence on production is due to the water stress factor, followed by the thermal factor: in both cases, the correlations are highly significant. Figure 4 shows the final assessment of the modelled inter-annual GPP variations against the measured tree-ring widths. The accordance between the two detrended data series is moderately high and statistically significant.

Case study of Molise
The annual evolution of precipitation and potential evapotranspiration computed over the 15-year period in Molise is shown in Figure 5; when compared to the Maremma study area, the climatic condition is different, being annual rainfall about 1600 (±425) mm and PET about 900 mm. As expected, slightly dry conditions are evident only during the summer period (from June to August), but the high soil water recharge from autumn to spring prevents the occurrence of heavy stress in the investigated forests. This is evident analysing the mean annual NDVI evolution, which is typical for a deciduous species with a minimum in winter and a nearly stable plateau during the growing season ( Figure 5). This last feature indicated an only marginal occurrence of summer water stress and provided the basis for the deactivation of the water stress factor (Maselli et al., 2014b).
The main drivers affecting the inter-annual GPP variations are shown in Figure 6a,b. The annual spring Tcor shows an increase from 1999 towards 2013. The same evolution, even more marked, is evident in the average NDVI of each year, meaning that the higher spring temperature enables green biomass to increase during the same period (Figure 6a).
The growing season PAR has a general tendency to decrease (Figure 6b). This obviously affects the GPP of the same season that, during the 15 years, is nearly stable. Its average is around 1500 g C m −2 year −1 and shows inter-annual variations which are in good accordance with those of the above-mentioned PAR.    Table 2 provides indications on the influence of all model drivers on forest growth during the investigated period. Both PAR and NDVI are highly correlated with tree-ring widths, while the opposite is for Tcor and Cws. The increasing temperature and the water stress factors act positively on NDVI but have only a limited effect on forest production due to the concomitant reduction of solar radiation.
The final evaluation of the detrended inter-annual GPP variations against tree-ring widths is shown in Figure 7. The accordance between the two data series is very good and the correlation is highly significant.

Discussion and conclusions
The estimation of forest production is a task achieved by numerous investigations carried out at different spatial and temporal scales. Most of the studies that use remotely sensed data are based on Monteith (1972), starting from the most traditional models, i.e. CASA (Field et al., 1995), to the MODIS method,  which routinely provides a 1-km operational product (Running et al., 2004). These models adopt different correction factors to account for the main environmental constrains, i.e. thermal and water stress, which down-regulate the maximum radiation use efficiency. The main model drivers can act differently, according to the eco-climatic characteristics of the investigated areas (Nemani et al., 2003). This is evident from global studies on the estimation of forest productivity (Running et al., 2004) and particularly from studies highlighting the impacts of recent changes on different model drivers (e.g. Boisvenue & Running, 2006). It has in fact been recognized that vegetation is generally adapted to the local environment in which it grows and its production is mostly constrained by water availability and temperature. For example, Nemani et al. (2003) found that the temperature increase has a positive impact on the forest production, especially at higher latitudes, while for Mediterranean areas, Cherubini et al. (2003), analysing tree rings, highlighted the importance of rainfall, especially if fallen in periods during which soil can recharge its water content. The importance of the water stress factor in Mediterranean areas was further highlighted by Maselli et al. (2014a) and Gilabert et al. (2015), who also indicated a reduced influence of the thermal factor. This issue has currently been investigated by analysing two case studies representative of different environmental conditions and climatic evolutions. In Maremma, which is a water-limited area, the precipitation fallen in late winter and spring is particularly important in determining annual GPP. These results are in accordance with those of Maselli et al. (2014a), who found that the influence of spring Cws increases along a gradient of aridity. The role of spring rainfall as a major ecological driver of annual tree growth is confirmed by Gandolfo (1999), who showed that earlywood is affected primarily by the precipitation occurring during the months preceding the start of vegetative season. More generally, the importance of including water stress factors for the simulation of GPP in semi-arid areas is supported by other studies which investigated the possibility of deriving them from remote-sensing imagery (Coops, Hilker, Hall, Nichol, & Drolet, 2010;Moreno et al., 2014;Sanchez-Ruiz et al., 2016). This, however, does not exclude the influence of other environmental factors which can be locally very relevant. An example is given by the salty water table and the pathogens' attacks affecting the Maremma pine forest (Teobaldelli et al., 2004), which can explain the only moderate correlation found between detrended treering widths and modelled GPP. Globally, the annual GPP of this forest was decreasing during the investigated 15-year period, as a consequence of the decreasing spring rainfall and summer NDVI.
In Molise, which corresponds to a temperature/ radiation limited area, summer water stress has only marginal importance in determining annual tree growth. This confirms results found by Tognetti, Lombardi, Lasserre, Cherubini and Marchetti (2014) in European beech forests across a latitudinal transect in Italy over the twentieth century. These authors observed that the enhancement of basal area increment (derived from tree-ring analysis) and intrinsic water-use efficiency (derived from stable isotope analysis) were uncoupled with the estimated drought index (based on precipitation sums and potential evapotranspiration). A similarly low influence of Cws on the annual GPP of European beech forests was found by Chiesi et al. (2016).
The applied GPP estimation approach requires the deactivation of Cws when water is supplied by sources additional to precipitation, such as irrigation or a water table. In the current case, Cws was deactivated after analysing the mean NDVI evolution of the beech forest, which, being not affected by summer Figure 7. Scatter plot between detrended tree-ring widths (X axis) and estimated annual GPP (Y axis) for Molise (**highly significant correlation, P < 0.01).
dryness, was indicative of the persistence of soil water. A similar deactivation was applied by Maselli et al. (2014b) to simulate the actual evapotranspiration of a presumably irrigated cropped field in Central Italy where the fractional vegetation cover was high during the period from May to September. Other criteria could be used as alternative indicators of unstressed conditions, still based on the uncoupling between rainfall occurrence and seasonal NDVI development.
The consequence of deactivating Cws is currently verified in other Italian forest areas by testing modified C-Fix against the eddy covariance observations described in . Better results are achieved in most cases where water remains available for vegetation also during summer (unpublished data). For example, in the mountain beech forest of Collelongo (central Italy), which is in an eco-climatic condition similar to Molise, the r 2 between observed and estimated daily GPP improves from 0.802 to 0.921. In the coastal plain of Castelporziano, where a water table underlies a holm oak forest, the deactivation of Cws improves the same coefficient from 0.815 to 0.850.
In these conditions, the temperature scalar becomes more influential, particularly during spring, when it has maximum effect on the start of the growing season. The same conclusion was obtained by Chiesi et al. (2016), who demonstrated that the inter-annual GPP variations of European beech forests are mainly dependent on minimum temperature around the beginning of the growing season. More generally, the influence of spring temperature on the annual growth of temperate-humid environments is supported by numerous studies, e.g. Delpierre et al. (2009) and Muraoka et al. (2010). A direct link between spring temperature and phenology is recognized, since warmer temperatures determine an earlier leaf unfolding and new wood formation, globally leading to an increase in the photosynthetic activity. In the current case, the increasing trend of mean spring air temperature determined a lengthening of the growing season, which was indicated by an increase of annual NDVI. These phenomena, however, were concomitant with a negative radiation trend presumably related to increasing cloud cover, which led to a general stability of estimated GPP.
Globally, the significant correlations found in the two case studies between GPP and tree-ring widths demonstrate the capability of the applied modelling approach to reproduce the inter-annual growth variations of forests limited by different environmental constraints. The ecological significance of the identified GPP evolutions is obviously limited to the periods examined, which are consequent on the length of the available data series. Important implications can, however, be drawn within the current debate about the impacts of climatic changes on Mediterranean forest ecosystems. The study, in fact, indicates that a warmer climate can have different effects depending on the eco-climatic characteristics of the observed areas and the interrelations between meteorological factors. More specifically, the expected temperature increase can act differently in environments limited by water (e.g. Mediterranean areas) and radiation (e.g. mountain zones), and these patterns can be further complicated by possible changes of rainfall and cloudiness, which are difficult to predict at appropriate spatial and temporal scales.