Study area

The study area is located in “Sierra de Los Filabres” (37º 22’ N, 2º 50’ W, between 1300 and 2186 m.a.s.l) (hereafter abbreviated as Filabres) (Fig. S1 [supplementary]). Mean annual precipitation is 320 mm, temperatures are moderately mild (13.1ºC, 1000 m a.s.l) for the 1940–2007 period with several drought events during this period (1990–1995, 1999, 2005; see more details in Sánchez-Salguero et al. 2012a). Soils are developed on schists and quartzites and they have a loam and silty loam texture (average composition 30–35% sand, 40–45% silt, 15–20% clay). Soil depth is 45–150 cm and available soil water content between 100 and 150 mm. Soil information was obtained from soil cartography at 1:100.000 scale (Aguilar et al., 1987). This information was used to estimate plot soil data as well as a test soil pit (Rosa et al., 1984) located in the study area. Dominant soils are xerorthents regosols and topography is characterized by steep slopes (>35%) (Aguilar et al., 1987). The forest area is dominated by Pinus sylvestris L. (hereafter as Scots pine; covering 7507 ha) and Pinus nigra Arn. subsp. salzmannii (Dunal) (hereafter as Black pine, covering 5900 ha) plantations. There are a reduced areas of native pine stands (Navarro-Cerrillo et al., 2007). Plantations were established using subsoiling as ground preparation between 1970 and 1976), but 1973 (Scots pine) and 1976 (Black pine)were selected as reference years. The planting density was 2000 trees ha^–1 and the current density ranges between 900 to 1000 trees ha^–1. Basal area range from 20.73 m² ha^–1 to 26.37 m² ha^–1 (Sánchez-Salguero et al. 2012a; Table S1 [supplementary]).

Field data collection

A stratified sampling was carried out on the basis of species and the degree of decline across the whole area of the plantations (Dobbertin & Brang, 2001). Accordingly, 9 plots of Scots pine and 9 plots of Black pine were established in July 2009 (Fig. S1 [supplementary]). Crown defoliation status was estimated according to ICP forests’ crown defoliation guide, aggregated at two different decline classes: class 1 (slightly declining trees and trees without evident defoliation <50%), and class 2 (declining trees 50-70%) (Dobbertin, 2005; Sánchez-Salguero et al., 2012 a, b). This assumption was supported by physiological data since stomatal conductance and water potential were significantly lower in defoliated than in non-defoliated trees in both species (Hernandez-Clemente et al., 2011). In each plot, all the trees with a diameter at breast height (1.3 m above ground level, dbh) greater than or equal to 10 cm were measured in circular plots with a radius of 15 m (706 m²). For every tree in each sample plot, two dbh measurements were taken at right angles with a tree calliper. Total tree height was measured with a Vertex III hypsometer (Table 1).

Table 1. Definitions, symbols, units, values, and sources for the parameters used in 3-PG model for Scots pine (Pinus sylvestris L.)and Black pine (Pinus nigra Arnold) in Sierra de los Filabres (Southern Spain)

Dominant trees of each plot were selected based on diameter at breast height (DBH) ≥20 cm and age ≥30 years, and a total of 45 trees were harvested. One wood disk per tree ca 5 cm of thick was cut at 1.30 m of tree height. After an acclimation period of 4 weeks in a chamber at 25-30 ºC, the wood discs were sanded with progressively finer grades of sandpaper until the wood anatomical elements were visible in transverse section and then scanned at 3200 dpi using an Epson Perfection V750 Pro scanner^© (Seiko Epson Corp., Nagano, Japan). Tree-ring width was measured to the nearest 0.01 mm along the two radii of each section using WinDendro^© (Regents Instruments, Quebec, Canada). Tree-ring series were dated following standard procedures (Stokes & Smiley, 1996). In order to detect dating and measurement errors, width series were checked with COFECHA software (Grissino-Mayer, 2001). We did not detrend chronologies to filter out eventual decrease of increment with age (age trend) since this operation could partly remove growth signal related to more short-term changes in trees` crown status (for further information see Sánchez-Salguero et al., 2012a). For the 3PG model of Basal Area Increment (hereafter abbreviated as BAI), we converted tree-ring width into basal-area increment, overcoming the problem of declining growth in bigger trees, using the formula (BAI = π (R²_t – R²_t–1)) where R is the radius of the tree and t is the year of tree-ring formation. To obtain values of overbark diameter and its evolution over time, a fixed ratio of bark thickness to diameter was considered.

The diurnal time course of stomatal conductance to water vapour (gL, CIRAS–1 instrument, PP Systems, Hitchin Herts, UK; every 90 min between 06.00 and 19.00 h (26^th July 2009) xylem water potential (Ψ) (Predawn Y_pd, 4:00 GTM and midday Y_m, 12:00 GTM; pressure chamber, SKPM 1400, Skye Instruments, UK) were measured in five trees per plot.

In each plot, LAI was non-destructively measured using a widely used optical instrument, the Plant Canopy Analyzer LAI–2000 (LICOR Inc., Lincoln, Nebraska, USA). The LAI–2000 measurements were performed in the sample plots during July 2009 using one intercalibrated instrument for all measurements. The unit was first located in a large open area to obtain the above-canopy measurements, and immediately the instrument was carried along a plot transect (10 m) to make the below-canopy measurements. All measures were done at dawn conditions. To minimize the contribution of the understory, LAI–2000 measurements were taken at 1.5 m above the ground. LAI was calculated as the arithmetic average of four LAI measures taken along the transect. LAI measured using LAI–2000 corresponds to plant area index (PAI) including photosynthetic and non-photosynthetic components. We computed effective LAI according to the methods used in Dufrêne & Breda (1995).

3-PG model and site specific parameterization

For this study, the 3-PG model was used to predict stand growth of Scots pineand Black pine plantations. The data required to run the model are organized in four categories: climate data, site-specific factors, initial conditions in the plantation and 3-PG physiological parameter relating to the species under study (Sands, 2004).

Site-specific climate data required by the model were obtained from thirty local meteorological stations with long, continuous records from the study area (distance less than 30 km in all cases) (Spanish Meteorological Agency; http://www.aemet.es/, see Sánchez-Salguero et al., 2010) including minimum, maximum, and average monthly air temperature (ºC), monthly rainfall (mm), number of rain and frost days per month. Direct measurements of vapour pressure deficit (VPD, mbar) were not available from the weather station used and VPD was estimated using the equations included in the 3-PG model (Landsberg, 1986). Solar radiation data (MJ m^–2 day^–1) were also not available from the permanent weather stations near our study sites, solar radiation was calculated using the Thornton-Running model by the software package developed by MT-CLIM for Excel, (available online: from http://www.ntsg.umt.edu/bioclimatology/mtclim). This produces solar radiation estimates based on site latitude, altitude, and daily precipitation and minimum and maximum temperatures (Thornton & Running, 1999; Amichev et al., 2011). Although the 3-PG model can be run using long-term monthly averages, we used current monthly weather data to account for drought events, in particular, during periods of relatively high mean temperature and high solar radiation that occurred during 1983, 1990–1995, 1999, 2001 and 2005. Climatic and site characteristics datasets were implemented in a Geographic Information System (GIS) environment. Spatial interpolation based on inverse distance weighting and spline, were used to predict monthly and average annual precipitation, as well as average annual, minimum and maximum temperatures with a spatial resolution of 10 m (Sánchez-Salguero et al., 2013).

The soil information needed for 3-PG initialization were obtained from soil cartography at 1:100.000 scale (Aguilar et al., 1987) and a test soil pit (Rosa et al., 1984) located in the study area. This information was used to parameterize the available soil water (ASW), and soil texture class variables required by 3-PG (Table 1). Soil texture was described as clay loam. The range of ASW was from 100 to 150 mm. We consider previous works using empirical expressions based on some measurements of soil depth and texture which have demonstrated to be useful in some recent works with 3-PG in Spain (e.g. Pérez-Cruzado et al., 2011; Vega-Nieva et al., 2013). Fertility rating (FR) was estimated based on the best match between model predictions and the measured average stem diameter, keeping all other parameters unchanged (Sampson et al., 2006; Rodríguez-Suárez et al., 2010). In this study, low values for FR were chosen, FR = 0.55, in healthy and decline stands, because of limitations imposed by soil texture and depth.

Species specific parameterization

Species-specific properties and features are controlled in the 3-PG model by different parameters grouped into five main categories: (i) biomass partitioning and turnover; (ii) growth modifiers; (iii) stem mortality and self-thinning; (iv) canopy structure and processes; and (v) wood and stand properties. We used four methods to parameterize each of the four species-specific sections of the model. First, we used empirical observations from permanent plots established on our study sites (Hernandez-Clemente et al., 2011); second, we used Scots pine and Black pine data from the literature (Patenaude et al., 2008); third, default 3-PG parameters published for Scots pine (Landsberg et al., 2005) and for Black pine (Patenaude et al., 2008); and fourth, we varied the remaining parameters to fit model output to DBH observations (Table 1). Twelve of the 3-PG parameters in Table 1 were parameterized using field stand observations, two with data reported in the literature, and three parameters were parameterized by fitting DBH data. We used the default 3-PG values for the remaining 31 parameters (Table 1).

The initial stem (WS), root (WR), and foliage (WF) biomass were estimated based on seedling weight and plantation density, with values of WF=0.068/0.089 tons ha^–1, WR= 0.058/0.071 tons ha^–1, and WS=0.041/0.050 tons ha^–1 for Scots pineand Black pine respectively (Montero et al., 2005). Tree diameters and BAI were derived from dendrochronological data. We parameterized the 3-PG tree allometric equations using a nation-wide Scots pineand Black pinestem biomass equations (Montero et al., 2005) including stem weight W_s (kg) = 0.0215 * (DBH, cm)^2.7184 for Scots pine(R²=0.981), and W_s (kg) = 0.043808 * (DBH, cm)^2.4975 for Black pine(R²=0.989).

For the simulation model with adaptive silvicultural treatments (e.g. uniform thinning), we considered a density reduction from the initial 2000 trees ha^–1 (in 1973 or 1976) to 1400 trees ha^–1 (in 1980) and to 1000 trees ha^–1 (in 2003).

Model sensitivity analysis

Once, the model was calibrated using an initial set of parameters (Table 1) and tested against observations (healthy trees), we performed sensitivity analyses to assess the effects of defoliation on 3-PG predictions for decline areas. In order to assess the relative influence of different model parameters on model results related to forest decline (DBH and BAI), we selected seven of the 3-PG parameters, likely to be related to drought stress. These were (Table 1): stem partitioning ratio (pFS₂₀), maximum litterfall rate (γ_Fx), maximum canopy conductance (g_Cx), specific leaf area for mature aged stands (σ₁), age at which specific leaf area = ½ (σ₀ + σ₁), age at full canopy cover (t_c), and canopy boundary layer conductance (g_B). The initial values for these parameters was based on field observations (Figs. S2 and S3 [supplementaries]). These simulations were based on empirically changes (increase and decrease) of those seven 3-PG parameter values (Table 2) after a first drought period (1983; see Sánchez-Salguero et al. 2012a). For this analysis we used a 32 year (1973–2008) and a 29 year (1976–2008, Black pine) annual time series of DBH and BAI for Scots pine and Black pine, respectively. Finally, the 3-PG model structure enables to study management intervention. Therefore, it was used to evaluate the impact of different thinning intensities on Scots pineand Black pine growth. Specifically, simulations were conducted at two levels of thinning intensities (remove of 50% and 60% of trees).

Table 2. Sensitivity analysis results for Scots pineand Black pine in Sierra de los Filabres (Southern Spain) shown as DBH values at age 35 yr and 32 yr respectively , produced by the 3PG model assuming 20 and 40% decrease and increase of the value of eleven 3PG model parameters

The model version used in this study was the 3-PGpjs2.5 (Sands & Landsberg, 2001) which is implemented as a Microsoft Excel spreadsheet with a user-interface that facilitates data entry and interpretation of results (FS 599_ALL EX + ANS 2015RHW.xlsm).

Assessment of accuracy of model performance

The outputs of 3-PG include stand attributes (average stem DBH, basal area ha^–1, volume ha^–1, biomass ha^–1, stand density and mortality) as well as physiological attributes (leaf area index, NPP, gross primary production, and transpiration). 3-PG results were validated by comparing measured and modelled data of average diameter at breast height (DBH, cm), basal area increment (BAI, cm² yr^–1), and leaf area index (LAI, m² m^–2) for healthy and decline conditions. We used DBH and BAI measurements from our testing data set to compute the bias of 3-PG predictions, as well as the R-square of linear regressions of predicted (dependent variable) versus observed (independent variable) data. Model error was evaluated using root mean squared error (RMSE) (Willmott, 1981).

Estimation of diameter at breast height growth on healthy trees

Nominal values applied for the calibration of the 3-PG model are included in Table 1 such as: site specific functions, and climate data from the plantation under study. Based on this model, we predicted mean diameter at breast height (DBH, cm), basal area increment (BAI, cm² yr^–1) and leaf area index (LAI, m² m^–2) for Scots pineand Black pinetrees and we compared these results with field observations from permanent plots. This analysis was done to establish a performance baseline of 3-PG for plantations of Scots pineand Black pinegrowth predictions.

Figure 1 shows the comparison between observed and simulated values of DBH growth for healthy trees of Scots pineplantation during a 35-year period (1973–2008 years). Predicted DBH was in near perfect agreement with observed measurements. At age 32, the Scots pinestand had a mean DBH of 16.0 cm while the model predicted 16.2 cm. A good correlation between measured and estimated DBH data was obtained, with an R² higher than 0.99, and an RMSE of 0.33 cm (Fig. 3a).

Regarding BAI, there was a lower agreement between observed and simulated values (Fig. 1b). The 3-PG model tended to underestimate tree BAI at younger stand ages during the early development stage (15 years old) due to the higher variation of our BAI observations (Fig. 1b). However, the large differences between modelled and observed data during early stand development had largely disappeared by middle-age period (13 years old). Overall, there was poor correspondence between modelled and observed BAI with a lower regression R-square values (R² = 0.83, RMSE = 1.54) (Fig. 3b).

Figure 2 shows the comparison between observed and simulated values of DBH growth for Black pineduring a 32-year period (1976–2008). As with the results obtained for Scots pine, predicted DBH was in very good agreement with observed measurements in Black pine(means DBH=15.5 cm while the model predicted 15.9 cm). The regression between modelled and observed DBH, had an R² higher than 0.99 and an RMSE of 0.35 cm (Fig. 3c). There was a lower agreement between observed and simulated values for BAI (Fig. 2d). The 3-PG model tended to underestimate tree BAI at younger stand ages for the early development stage (7 years old) due to the higher variation of our BAI observations (Fig. 2d). However, the large differences between modelled and observed during early stand development were more marked by middle-age period (10 years old). There was a low correspondence between modelled and observed BAI (R² = 0.68, RMSE = 1.52) (Fig. 3d).

3-PG model sensitivity analysis to simulate forest decline trees

In order to identify the relative influence of different model variables on the ability to predict results for forest decline trees (DBH and BAI), field data were used to fit site specific empirical functions. These simulations were based on changes either positive or negative of seven 3-PG parameter values (Table 2) after a first drought period (1983).

Our sensitivity analysis results indicated that the most pronounced effects on DBH predictions for decline trees could be due to parameterizing the foliage: stem partitioning ratio (pFS₂₀), maximum litterfall rate (γ_Fx), maximum canopy conductance (g_Cx), specific leaf area for mature aged stands (σ₁), age at which specific leaf area = ½ (σ₀ + σ₁), age at full canopy cover (t_c), and canopy boundary layer conductance (g_B) (Table 2).

By using this new set of parameters related to decline process (Table 3), estimations were able to appropriately predict the DBH of decline trees during the 1983–2008 periods for both species (Figs. 1c and 2c). At age 32, the Scots pinestand had 14.1 cm while the model predicted 14.4 cm, and Black pinestand had 14.2 cm while the model predicted 14.2 cm. Model estimates were highly correlated with observed measurements with slopes near unity with an R² of 0.99, and an RMSE of 0.31 cm (Figs. 3c, 2 g, 3g). On decline trees, 3-PG predicted BAI values showed a lower correlation with field data (Figs. 1d, 2d, 3f and 3h) (R² = 0.69, RMSE = 1.84 cm² yr^–1, and R² = 0.36, RMSE = 2.57 cm² yr^–1 respectively).

Table 3. Values for the parameters used in 3-PG model for Scots pine (Pinus sylvestris L.)and Black pine (Pinus nigra Arnold) in Sierra de los Filabres (Southern Spain) related to decline process under stress conditions

Leaf area index (LAI) was predicted for both species and physiological status (Fig. 4). It was always much higher for the healthy trees than the decline trees. Predicted highest LAI value occurred between ages 16–17 (P. sylvestris in 1992, LAI = 1.14 m² m^–2; P. nigra in 1989, LAI = 1.47 m² m^–2). There was a substantial difference between healthy (P. sylvestris, LAI = 0.76 m² m^–2; P. nigra LAI = 0.74 m² m^–2)and declining trees (P. sylvestris LAI = 0.29 m² m^–2; P. nigra LAI = 0.37 m² m^–2), where at age 32-35 observed LAI was much lower than predicted LAI. (Fig. 4). Although we only had one year of measured LAI data for comparison purposes, the observed LAI showed a limited agreement with predicted LAI.

Sensitivity analysis of silvicultural treatments

Model predictions for DBH of Scots pineunder two different thinning intensities (50% and 60%) were consistently higher than predicted measurements, with an estimated DBH increment of 1.3 cm after 7 years (17.4/17.6 cm versus 16.3 cm) (Fig. 5). However, predictions after thinning were relatively closer on both silvicultural treatments, with a trend indicating a larger DBH overestimation as stands aged.

3PG-model to assess decline processes

In this study, we applied the 3-PG model to describe the physiological process behind forest decline in pine forests (McDowell et al., 2008; Williams et al., 2013). Supporting our hypothesis, a comparison of the model simulations with growth measurements indicated that a fixed physiological parameter set in 3-PG was able to predict DBH, BAI and to a lesser degree LAI with reasonable accuracy under healthy and decline conditions for plantations of Scots pineand Black pinein Southern Spain. We were then able to reproduce the observed growth reduction indicated by dendrochronology data by reducing leaf conductance and the available foliage area, which lead to less photosynthetic capacity and slower recovery after a drought stress.

Previous studies have shown that 3-PG could produce accurate estimates of Scots pineand Black pineDBH and were able to simulate growth patterns under different ecological conditions (Esprey et al., 2004; Landsberg et al., 2005; Patenaude et al., 2008; Xenakis et al., 2008). Nominal values applied for this study were based on previous studies (Sands & Landsberg, 2001; Patenaude et al., 2008). This paper presents a novel approach in the application of 3-PG model to accurately estimate the effects on growth of forest decline processes in pine stand based on the analysis of three parameters (DBH, BAI and LAI).

In the first approach, we modified sixteen of the 3-PG variables using field stand observations (Table 1) (Sands & Landsberg, 2001; Patenaude et al., 2008). All of these parameters are directly related to photosynthesis, stomata response and carbon accumulation, as they affect foliage biomass production and litter input (Landsberg & Sands, 2010). The model was able to predict DBH with reasonable accuracy for the base-line physiological set; although, there were some systematic differences between estimated and measured values of DBH. The model tended to under-predict DBH during early stand development and then over-predict later (Figs. 1 and 2). However, the good correlations between predicted and observed data may be explained for two reasons. First, we used current monthly weather data instead of monthly average, which allowed accounting for drought and heat discrete events as previously noted on other studies (Almeida et al., 2004; Bryars et al., 2013).

Second, under the particular study conditions (even-age plantations, regular tree distribution, no thermal growth limitations), pine growth appears to be much more dependent on leaf area development and the quantity of intercepted solar radiation than on changes in the rates or efficiencies of specific physiological processes (Beadle et al., 1985a,b,c; Samuelson et al., 2004).

The accuracy of the BAI estimates was lower but varied during different stages of stand development (Fig. 3). The low model accuracy for BAI was due to poor predictions during early and late stand development. It is noteworthy that in the Filabres area, P. sylvestris grew (increased biomass) well in the wet 1970s, when conditions were much more favourable for establishment, than in the drier 1990s. In the case of Scots pine, an introduced species there, the plantations are located well beyond the natural southern boundary for the species’ distribution area (Sánchez-Salguero et al., 2012a). These results are biologically reasonable and illustrated how the processes embedded in the 3-PG model controlled tree biomass production based on parameter values (Sands, 2004; Landsberg & Sands, 2010).

Once the base growth line was established, the physiological parameter set was modified to describe forest decline process (Hernández-Clemente et al., 2011; Sánchez-Salguero et al. 2012a,b). Previous studies have established physiological parameters related to forest decline (McDowell et al., 2008). On decline trees, values of stem partitioning ratio (pFS₂₀), maximum litterfall rate (γ_Fx), maximum canopy conductance (g_Cx), specific leaf area for mature aged stands (σ₁), age at which specific leaf area = ½ (σ₀ + σ₁), age at full canopy cover (t_c), and canopy boundary layer conductance (g_B) were different from those used on other simulations for pine species (Landsberg et al., 2005; Patenaude et al., 2008; van Oijen et al., 2013), and they were important for DBH and BAI predictions under drought stress (Esprey et al., 2004). However, the use of those parameters led to values which are in good agreement with field data (Figs. S2 and S3, [supplementaries]). The litterfall rate was fitted to match age-related decline of growth observed by Sánchez-Salguero et al. (2012a), whereas values used for the specific area were 4.6 m² kg^–1 and 4.3 m² kg^–1 for Scots pine mature leaves and values of 4.2 m² kg^–1 and 3.5 m² kg^–1 for Black pine were within the lower range of those measured in England (Mencuccini & Bonosi, 2001; Patenaude et al., 2008). More restrictive physiological values (i.e., maximum canopy conductance (g_Cx) and canopy boundary layer conductance (g_B) affecting water photosynthetic status) led to a different biomass allocation to the stem or to foliage, resulting in lower DBH and BAI estimates. On the other hand, the inaccuracy in predicting LAI in older stand ages seems to be more related to an inaccurate prediction of stem allometry or growth allocation (Fig. 4; Beadle et al., 1982); although, measured LAI data for the plots used in this study were limited. Again it is important to point out that this relationship was biologically reasonable and consisted with the processes are embedded in the 3-PG model (Landsberg et al., 2005). However, it is important to recognize that it was possible to obtain similar values with different combinations of values for the tuned parameters, suggesting that a deeper analysis, including sensitivity analysis, will be needed. The sensitivity analysis of those parameters on growth under forest decline process has demonstrated the importance of selecting an adequate physiological set value to obtain an accurate model fit under stress conditions.

Silvicultural treatments

On the adaptive silviculture treatments, the model tended to show a similar increase of DBH during post thinning stand development (Fig. 5) for both thinning levels. The average variation of DBH growth 10 years after thinning was estimated as 2.6 cm higher than the current stand. This inter-rotation analysis highlighted the strong influence of management practices on pine forest growth and drought adaptation (Fontes et al., 2010). The increase in the model DBH growth estimation could be due to water availability related to tree density (Martín-Benito et al., 2010), effects on stem C allocation (Oleksyn et al., 1999), and reduction of the ratio of foliage: stem biomass following severe thinning (Vanninen et al., 1996). However, it has been suggested that silviculture has little effect on stem allometry (Retzlaff et al., 2001). Since there are few empirical estimates of the impact of silvicultural treatments on forest decline response (Millar et al., 2007; Sánchez-Salguero et al., 2012b, 2013), field data of foliage and stem partitioning coefficients from component growth, litterfall, and stem mortality should be collected.