Coupling Photosynthetic Measurements with Biometric Data to Estimate Gross Primary Productivity (GPP) in Mediterranean Pine Forests of Different Post-Fire Age

Abstract

Quantification of forest Gross Primary Productivity (GPP) is important for understanding ecosystem function and designing appropriate carbon mitigation strategies. Coupling forest biometric data with canopy photosynthesis models can provide a means to simulate GPP across different stand ages. In this study we developed a simple framework to integrate biometric and leaf gas-exchange measurements, and to estimate GPP across four Mediterranean pine forests of different post-fire age. We used three different methods to estimate the Leaf Area Index (LAI) of the stands, and monthly gas exchange data to calibrate the photosynthetic light response of the leaves. Upscaling of carbon sequestration at the canopy level was made by implementing a Big Leaf and a Sun/Shade model, using both average and variant (monthly) photosynthetic capacity values. The Big Leaf model simulations systematically underestimated GPP compared to the Sun/Shade model simulations. Our simulations suggest an increasing GPP with age up to a stand maturity stage. The shape of the GPP trend with stand age was not affected by the method used to parameterise the model. At the scale of our study, variability in stand and canopy structure among the study sites seems to be the key determinant of GPP.

1. Introduction

Forest ecosystems are vital for the survival of human kind as they provide material, economical and recreational resources [1], as well as fresh water, soil protection and other ecosystem services [2]. Currently many scholars focus on the role of carbon (C) cycling within forests, and the way this process is affected by global change [3] and interacts with the climate system [4]. To understand C dynamics, estimation of the rate of C sequestration and C stock change in forest ecosystems is crucial for ecological and global change research and policy.

Estimation of ecosystem C fluxes is usually made by accounting for the different components of ecosystem productivity [5]. The terrestrial Gross Primary Productivity (GPP) is the amount of C that plants enchain through photosynthesis in the form of carbon dioxide (CO₂). The annual and inter-annual GPP variation provides information about ecosystem dynamics and the rhythms of CO₂ transposition from the atmosphere into the biosphere [6]. GPP is controlled by factors that directly regulate photosynthesis as well as by stand properties that indirectly regulate plant growth [7]. Abiotic factors that influence photosynthesis operate on short-term (hourly), mid-term (daily to monthly) and long-term (years) time scales. At short time scales, solar radiation [8], temperature [9] and humidity [10] fluctuations are the main factors that regulate photosynthesis. At mid-term scales, water and temperature fluctuations [11,12] are considered the key determinants of photosynthetic C accumulation, while at longer time scales nutrient availability and disturbances can have more prominent effects [13,14]. In forest ecosystems, stand density directly affects the Leaf Area Index (LAI), the amount of foliage area per unit ground area [15,16]. LAI is determined from the stand morphology and characteristics [17], and it is a major biotic factor controlling within-stand light availability [13,18], microclimate [19] and water balance [20], and thus the overall GPP of a stand [21]. At sites with infrequent disturbances, as stands mature and LAI reaches a local maximum [22], nutrient and water limitation may occur due to competition among trees, leading to reductions of stand-level GPP. However, disturbances such as wildfires and their regimes (frequency, interval and severity) play a crucial role in GPP dynamics, by resetting various ecosystem features [14]. Thus, a sound understanding of how abiotic factors, such as solar radiation, temperature and precipitation, and biotic factors, such as stand and canopy structure, control GPP is important to understand how forest ecosystems function on seasonal and interannual timescales.

Climate change is expected to affect the dynamics of Mediterranean forests [23]. Projected temperature increases [24,25], the uncertain distribution of water availability [26,27] and the interplay between drought and wildfire [28,29] are some of the most important agents that could shift the rate of C assimilation in Mediterranean forests. In areas surrounding the Mediterranean basin, mean annual temperature is expected to increase by up to 3–4 °C and precipitation to decrease by 20%, leading to shifts in the forest fires regime [30]. As a result, fires could disturb the consistency of forest ecosystems, creating areas with different stages of succession and species synthesis [31]. Typically, in the post-fire period, the trajectory of ecosystem productivity follows a single peak curve with the maximum of productivity achieved at middle age stands [7,32]. Following such scenarios, it is important to accurately simulate C fluxes to elucidate the role of Mediterranean forests in the global C cycle and to develop effective reforestation strategies for C sequestration.

Estimation of the C exchange from forest ecosystems is challenging due to the complex relationships between biotic and abiotic factors [33]. Currently, two methodologies are implemented to infer forest productivity, the “top–down” and the “bottom–up” approaches [33,34]. The “top–down” or meteorological approach attains the estimation using eddy-covariance flux towers [35]. This method is based on the continuous measurement of CO₂ fluxes between ecosystems and the atmosphere, without segregating and quantifying the fluxes within the different C reservoirs of the ecosystem, for example the living biomass, the dead biomass and the soil biomass [36]. Although this method gives accurate estimations of net ecosystem productivity (NEP), i.e., the amount of C fixed by GPP after removing losses from autotrophic and heterotrophic respiration, it does so without the information of the components’ contribution. Therefore the estimation of GPP from eddy-flux data is made under a set of assumptions regarding the autotrophic and heterotrophic respiration rates [6,35]. Moreover, in cases that eddy-covariance flux towers operate in complex terrains and vegetation sites, additional approaches and data, such as mesoscale modelling, remote sensing and/or ecophysiological measurements, are used to accurately partition C fluxes [37].

The “bottom–up” or ecological approach consists in using finer spatial scale data, including biometric measurements, such as foliage, branch, stem and root growth and litter and root turnover rates, as well as gas exchanges measurements (photosynthesis and respiration) to understand C fluxes within an ecosystem [34,38]. Biometric data can additionally provide information about the structure and the properties of the canopy and the stand. Detailed biometric data can be used to estimate Net Primary Productivity (NPP), i.e., the amount of C fixed by GPP after excluding autotrophic respiration [39]. Such data can also be integrated with photosynthetic and respiration measurements, to ultimately infer GPP and NEP [40]. In this perspective the ecological approach is more complex but does provide the means for upscaling from the leaf to the ecosystem scale [36,41]. In contrast with the “top–down” approach, “bottom–up” studies necessitate the frequent physical presence of the researcher in the field and the associated difficulty to deliver continuous data [42]. By combining the “top–down” and the “bottom–up” approaches, more accurate estimations of ecosystem productivity can be obtained [43], although an agreement between the two methods is generally achieved after comparing multiple years of concurrent measurements [34].

Estimation of Mediterranean forests GPP has been primarily pursued through eddy-flux measurements and modelling studies [44,45,46,47,48,49]. Based on eddy-flux measurements, estimation of GPP ranges from 2.49 to 4.60 g C m⁻² day⁻¹ [45,47,49] throughout the year. On the other hand, model simulations yield relatively lower GPP from 1.15 to 2.66 g C m⁻² day⁻¹ [44,46,48] for typical Mediterranean forests. The range of aboveground NPP in Mediterranean forests has been estimated through biometric studies (0.13–0.55 g C m⁻² day⁻¹) [50] while annual NEP estimates from forest inventories can reach a rate of 0.25 g C m⁻² day⁻¹ [51]. Based on inferences from litterfall data, some studies estimate the annual NPP of Mediterranean pine forests at 1.64 Mg C·ha⁻¹·yr⁻¹ [50]. However, to our knowledge a systematic implementation of the “bottom up” approach integrated with ecophysiological measurements has not been implemented in Mediterranean pine forests.

To understand the interplay between tree morphological and physiological processes with the biotic and abiotic conditions of the stand and ultimately quantify C fluxes, process-based forest models can be used as a framework to integrate information on different ecosystem processes [52]. In this study we integrated the “bottom–up” approach with simple canopy photosynthesis models to estimate GPP variation in Mediterranean Pine forests (Pinus brutia Ten.) across a post-fire chronosequence on Lesvos Island, Greece. The “bottom–up” approach provides an advantage for studying Pine forests productivity on Lesvos, due to the topographic complexity and the forest structural heterogeneity arising from the fire history on the island. In regions of topographic and vegetation heterogeneity the eddy covariance flux tower approach could be biased [38]. Thus, even though the “bottom–up” approach is more labour-intensive, it can provide a sound understanding of C fluxes within an ecosystem [38], especially at finer spatial scales.

The aim of the study is to explore the key factors that control GPP variation in Mediterranean Pine forests of different post-fire age. Based on forest dynamics theory we expected an increase in GPP with post-fire age, until stands reached a stage of maturity, followed by a gradual decrease in productivity with increasing stand age. We anticipated seasonal GPP variation within stands to be primarily controlled by local weather conditions and phenological variation in foliar properties. Across different stand development stages, competition for light and the ability of canopy to capture radiative energy, expressed through LAI, was expected to be the main driver of GPP. We used three different methods to quantify LAI, based on direct and optical measurements. In order to systematically quantify variation in estimated GPP related to these sources of methodological variation, we used a factorial simulation protocol, which quantified variation in GPP across different methods of LAI and photosynthetic capacity representation.

2. Materials and Methods

2.1. Study Sites

Four permanent measurement plots of different post-fire age have been established and monitored on a monthly basis since July 2019 across the Pinus brutia forest ecosystem on the island of Lesvos, Greece (Figure 1). Lesvos lies in the Northeast Aegean Sea, with a typical Mediterranean climate, characterised by long xero-thermic periods and mild wet winters. Dry season begins at mid-April and ends at early October. The average annual temperature on the island is approximately 18 °C, with the minimum monthly temperature in January at 6.8 °C and the maximum monthly temperature in mid-July at 31 °C. Mean total annual precipitation is 645 mm with most rainfall falls in December and January [53].

The four plots (

Table 1. Characteristics of the study plots. AMAL is the 13 year-old plot, PEV is the 40 year-old plot, LML is the 72 year-old plot and ACHL is the 92 year-old plot.

Plot	Lat 1	Lon 1	Stand Age (y)	Inclination %	Orientation
AMAL	39.02	26.59	13	10	SE
PEV	39.16	26.37	40	5	SW
LML	39.16	26.38	72	10	SW
ACHL	39.13	26.30	92	0	-

¹ Coordinates are given in decimal degrees.

) were established across a post-fire chronosequence (13 to 92 years old) on similar elevation (166 to 306 m asl) and soil parental material (Ophilithic bedrock), using regional forest wildfire and geological maps. Soil texture and depth are presented in Table A1. During plot establishment we tried to minimise effects from recent disturbances by avoiding areas with obvious grazing or fire signs. All plots have mild slope (under 10%) and are dominated by P. brutia with several woody species, such as Quercus coccifera L., Pistacia lentiscus L., Erica spp. and Cistus spp., found in the understory, rarely above 1 m high, and mainly present in the two younger plots.

To determine stand age in October 2019, we collected increment cores from at least 12 trees per site. From each tree, we extracted two cores to minimise the possible impact of reaction wood. Samples were collected at breast-high with the use of a 5 mm increment borer. After surface preparation, tree-ring width was measured to 0.01 mm using the Time Series Analysis and Presentation software (TSAP-Win) [54] and a LINTAB (Rinntech Inc. Heidelberg, Germany) measuring table. All samples were visually and statistically cross-dated using TSAP-Win to identify the possible presence of missing or false rings [55]. Samples from each site were synchronised in one single mean chronology per site, the age of which determined the mean age of each stand (Table 1). Detailed information on dendrochronological analysis is provided in ([56], ‘under review’).

Square plots of 900 m² were deployed and each tree above 1.3 m received a unique code for the facilitation of biometric measurements. For each tree diameter at breast height (dbh—cm) and height (H—m) were measured in July 2019 and monitored on an annual basis thereafter. The status (alive, dead) of each tree has also been recorded since the start of plot monitoring period.

2.2. Gas Exchange and Functional Traits Measurements

From July 2019 to December 2020, we measured leaf-level CO₂ exchange using a LICOR 6400 (LICOR Inc., Lincoln, NE, USA) infrared gas analyser, on a monthly basis, except for August 2019/2020, October 2019 and November 2020. All measurements were made on days with typical weather conditions based on previous days’ weather forecasts. We randomly selected trees within the stand and measured light response curves (Lc) for six to ten trees each month. Measurements were made for three pairs of fully developed and healthy needles, from both sunlit (four to eight Lc measurements) and shaded (two to four Lc measurements) branches. A fully sunlit or shaded branch from the selected tree was cut by climbing on the tree and/or using telescopic scissors, and immediately placed in a water bucket where it was recut prior to leaf gas-exchange measurements [57]. Only current year’s needles were selected, or needles of the year before, in cases that the current year’s needles were not fully developed, by identifying the small and clustered sterile scales that indicate the end of the previous and the beginning of the new growth unit [58]. Light response curves were made using the 6400–02B LED Light Source chamber with red and blue LEDs (665 and 470 nm, respectively). The chamber was set up at ambient conditions of humidity and temperature (provided by the L6400 external sensors), at 400 ppm CO₂ concentration, by using CO₂ cartridges, and systematically changing the rate of photosynthetic active radiation (PAR) between 2000, 1800, 1600, 1400, 1200, 1000, 800, 600, 400, 200, 100, 80, 60, 40, 20 and 0 μmole quanta m⁻² s⁻¹ to measure net photosynthesis (A_net (μmole m⁻² s⁻¹)). At each PAR level, three logs of A_net were taken after at least a two-minute period during which the needles were left at the specific PAR flux to reach an equilibrium. At the end of each light curve, needles were acclimated under dark conditions for five minutes and an additional reading of leaf dark respiration (Rd (μmole m⁻² s⁻¹)) was taken. All sets of Lc measurements began 30 min to 1 h after dawn and ended at least 2 h after noon. Because the area of the needles inserted in the chamber did not fully cover the chamber area, we marked the edge points and corrected our measurements for the leaf area enclosed.

Each day, the needles used in the gas-exchange measurements were transferred to the laboratory for additional leaf traits measurements. In the lab, needles were treated in wet and dark conditions. One day later we measured wet mass (LWW—g) and thickness (Lth—mm) before scanning the needles with a portable scanner (iScan 900 dpi). Leaf area (LA—mm²) of the whole needles as well the marked area of the needles enclosed in the chamber were estimated using the ImageJ image analysis software [59]. We then oven-dried the needles for 72 h at 70 °C and measured their dry weight (LDW—g). Leaf dry mass per area (LMA—g m⁻²) and Leaf Dry Matter Content (LDMC) were estimated using average LWW, LDW and needle area [57].

2.3. Stand Level Measurements

At the stand level we quantified the monthly litterfall and branch fall rates. Five traps were made and installed in each plot. Every trap had a square frame (0.5 m × 0.5 m) and four legs (0.3 m height) made by PVC pipes (φ 16). Legs mounted the structure 0.3 m above the ground to prevent needle sepsis during the wet season. A fiberglass window net was installed inside each frame so as to create an inner curve with a maximum depth of 0.2 m, in order to prevent wind from removing needles and needle sepsis during the wet season. Litterfall accumulated in the traps was collected monthly and then oven-dried for 48 h at 70 °C. Needles were separated from other forms of litter, such as cones, twigs, and bark specimens. The mass of the needles was recorded and divided by the total area of the litterfall traps, for estimation of litterfall mass per m².

The LAI of each plot was estimated using three different methods. In the first approach, namely the Litterfall Method (LFM), we used data from the litterfall traps [60,61] to estimate LAI by dividing the annual sum of litterfall mass per m² with the annual average LMA and multiplying by the average needle lifespan (1.5 years, [62], C. Sazeides observations). In the second method, namely the Hemispherical Photography Method (HPM), we divided each plot into nine subplots where we took two fisheye photographs with a Canon EOS-60D and a fisheye lens (sigma 4.5 mm circular fisheye), mounted on a tripod and placed 1 m above the forest floor. The camera was always levelled, and its top side was oriented northward. All pictures were taken under overcast conditions, early in the morning during the spring of 2020, with a fixed aperture (f = 7.9) and were slightly underexposed (2/3 f-stops) [63]. Images were analysed using the HemiView Software (HemiView version 2.1 Delta-T Devices Ltd. Burwell, UK) and an average LAI per plot was estimated. In the third approach, namely the Ceptometer Method (CMM), an average plot-level LAI was estimated by systematically measuring incoming PAR within the plot using the AccuPAR linear PAR ceptometer (Decagon Devices, Inc., Pullman, WA, USA) and comparing with measurements outside the canopy at full light conditions. CMM LAI estimation was made with measurements taken around solar noon and the leaf distribution parameter set to X = 1. In each plot, six parallel transects with 5 m distance from each other were identified, and PAR measurements were taken every 5 m, at 1 m above the forest floor, counting a total of 36 LAI values for each plot, which were averaged to estimate a mean plot-level LAI.

2.4. Environmental and Remote Sensing Data

Ambient air humidity (RH) and temperature (T) were monitored in each site using iButton Hygrochron DS1923-F5# (iButtonLink, Whitewater, WI, USA) at a 5 h time step. The iButtons were placed at 2 m height on the north-facing side of a single tree in each plot. The period covered was between July 2020 and July 2021. For each month (i) and study site the average monthly temperature (Ti) and relative humidity (RHi) and standard deviation were calculated and compared between stands (see Figure A1). There were non-significant differences at ambient conditions among the sites in the same month. As the study plots are located within a small distance from each other, variation in Ti’s and RHi’s between plots are probably attributed to topographic variation and differences in stand structure. Since Ti and RHi data only partly overlapped with the simulation period, they were only used to describe site-specific weather conditions.

Hourly data for direct normal irradiance and diffuse horizontal irradiance were downloaded from Solcast (Solcast, 2019, Global solar irradiance data and PV system power output data, URL https://solcast.com/ (accessed on 21 January 2021)), a global solar forecasting and historical solar irradiance data company. The two variables were summed up to calculate the total normal irradiance, as the maximum irradiance that canopy received. Total irradiance (Wm⁻²) was converted to PAR (μmole quanta m⁻² s⁻¹) by multiplying with a constant equal to 2.1479 [64]. In addition, absorbed canopy radiation (I_c) was estimated from total irradiance using a constant absorbance coefficient (a = 0.75) to account for the fraction of incoming solar radiation that is absorbed by the foliage [65].

In order to compare our simulations with an independent dataset, we extracted for each study site the global GPP from the moderate resolution (500 m) dataset of Zhang et al. [66] for the 2000–2016 period. These estimates are based on an improved light use efficiency model and are driven by satellite data from MODIS and climate data from NCEP Reanalysis II. The mean value and the standard deviation (between years) across the 17 year-long period were used as a baseline to compare with our simulated GPP.

2.5. Statistical Analysis

For each light response curve (per branch, month and plot) we first corrected the needles’ area enclosed in the chamber and then fitted the Michaelis–Menten (MM) model across the PAR range:

$$A_{net}=A_{sat}PAR/\left(Km+PAR\right)-R_d$$

(1)

where A_sat is the light saturated net photosynthetic rate (μmole m⁻² s⁻¹), K_m is the half saturation coefficient (μmole quanta m⁻² s⁻¹), and R_d is the dark respiration rate (μmole m⁻² s⁻¹). The model was fit to data from the light response curves using the Global Optimization by Differential Evolution algorithm (package DEoptim) [67] and minimising the root sum of squares [52].

To test for seasonal differences in LMA and the parameters of the MM model we used a mixed effect (package lme4) [68] model, using plot as a fixed effect and month as a random effect term, due to the fact that our interest regarding months lies in the variation among them rather than the specific effect of each level [69]. The importance of accounting for seasonal variation was tested after comparing the Akaike’s Information Criterion (AIC) of the mixed effect model with a simple linear model without the monthly random component term (ΔAIC < 2). A Wilcoxon test was used to explore for differences between the parameters of the MM model between sunlit and shaded branches. An asymptotic test (package cvequality) was used to evaluate the coefficient of variation (CV) equality among the average and variant model estimations.

2.6. Coupling Biometric and Gas Exchange Data to Simulate Stand Level GPP

To estimate the GPP in each study plot we integrated the monthly gas exchange and biometric data with two simple canopy photosynthesis models. In particular, we ran the Big Leaf [70] and the Sun/Shade [71] canopy photosynthesis models in six different setups for the site-specific abiotic and biotic condition of each study plot. The twelve different setups are summarised in

Table 2. Overview of the six different model setups. See Equation (1) for description of the light response model parameters. Leaf Area Index (LAI) estimated using the following: Hemispherical Photography method (HPM), Litterfall method (LFM) and Ceptometer Method (CMM).

Model	Photosynthetic Capacity Setups	LAI Setups
Big Leaf (BL) Model	Monthly (variant) photosynthetic light response parameters	HPM
		LFM
		CMM
	Annual (average) photosynthetic light response parameters	HPM
		LFM
		CMM
Sun/Shade (SS) Model	Monthly (variant) photosynthetic light response parameters	HPM
		LFM
		CMM
	Annual (average) photosynthetic light	HPM
		LFM
		CMM

.

The “Big Leaf” approach (BL) is a simple method for estimating canopy photosynthesis, without any canopy profile information required [72]. It is based on the assumption that the scaling procedure from a single leaf to the canopy is governed by a linear relationship, assuming that canopy response to environmental stimulations is the same as that of a leaf [73]. Thus, canopy is considered as a single layer, with all leaves receiving the same radiation [71]. Requirements for this model are few, in terms of field-measured parameters as well as computational cost [74]. However, the assumptions of the BL model do not take into accounts the decrease of irradiance throughout the canopy and the relative position of leaves within the canopy. De Pury and Farquhar [71] descripted the Sun/Shade model (SS) that separates the sunlit and the shaded part of canopy. The SS model estimates the fraction of canopy that is sunlit or shaded, as well as the irradiance those two parts receive, providing a more realistic estimation of canopy photosynthesis.

Below we describe the coupling of the BL and SS models with the biometric and solar radiation data of the study plots. All algorithm developments, graphs and statistical analyses were made in R [75] using the packages lubridate [76], dplyr [77], fishmethods [78], ggplot2 [79], ggpubr [80] and tidyverse [81].

For each month and in each study plot, a discrete light response curve was available, which accounted for potential acclimation of photosynthesis to environmental variation [82,83]. The monthly Lc curves were used in the varying photosynthetic capacity model setups, while a mean Lc fitted with all available data (for measurements across all months Figure A2) was used in the constant photosynthetic capacity model setups.

In all cases, hourly radiation data were used to force the models and simulate hourly A_net under the constant and varying photosynthetic capacity model setups, with R_d set to 0 during daylight hours. Simulations were made for the period between July 2019 and December 2020. The hourly time-step simulations were summarised to monthly and annual GPP estimates for model intercomparison and validation. The estimation of the annual GPP was made for the period between July 2019 to June 2020.

2.6.1. Big Leaf Model

Scaling up productivity from the leaf to the canopy level following the BL model setup was made using:

$$GPP=A_{net}\left[1-\exp \left(-k_b LAI\right)\right]/k_b$$

(2)

where exp is the exponential function and k_b is the beam radiation extinction coefficient of the canopy:

$$k_b=0.5/\sin \beta$$

(3)

and β is the solar elevation angle (in radians) and 0.5 is the ratio of projected area to surface area of a hemisphere. Hourly values of β were estimated from the astrocalc4r function (package fishmethods).

2.6.2. Sun/Shade Model

The application of the SS model required additional calculations. We initially simulated the fraction of sunlit canopy according to:

$$L_{sun}=[1-\exp \left(-k_b LAI\right)]/k_b$$

(4)

where L_sun is the fraction of sunlit LAI.

The shade fraction of the canopy (L_sh) was estimated from the difference between LAI and L_sun, on an hourly basis. The fraction of PAR absorbed from the sunlit part of the canopy was calculated using:

$$I_{csun}=I_{lb}+I_{ld}+I_{lbs}$$

(5)

The absorbed irradiance from the sunlit canopy (I_csun) was estimated as the sum of the direct beam (I_lb), diffuse irradiance (I_ld) and scattered beam irradiance (I_lbs) absorption (Equations (A1)–(A3)). The total absorbed irradiance from the canopy (I_c) was calculated from:

$$I_c=\left(1-rcb\right)I_{lb}\left[1-\exp \left(-k_{bs}LAI\right)\right]+\left(1-rcd\right)I_d\left[1-\exp \left(-k_{ds}LAI\right)\right]$$

(6)

with the difference between I_c and I_csun representing the irradiance absorbed by the shaded part of the canopy (I_csh). Equations for calculating rcb, rcd, k_ds and k_bs, are summarized in Table A2. Scaling up to the canopy level following the SS model setup was made using:

$$GPP=[A_{sat,sun}{L}_{sun}{I}_{csun}/\left(K_{m,sun}+L_{sun}{I}_{csun}\right)]+[A_{sat,sh}{L}_{sh}{I}_{csh}/\left(K_{m,sh}+L_{sh}{I}_{csh}\right)]$$

(7)

3. Results

3.1. Variation in Stand Structure across the Post-Fire Chronosequence

Variation of mean diameter at breast height among the study sites (

Table 3. Variation in stand structure across the study plots. Leaf Area Index (LAI) was estimated with three different methods (see text for details). Diameter at breast height (dbh) measured for trees with a diameter > 1 cm. Hemispherical Photography method (HPM), Litterfall method (LFM) and Ceptometer Method (CMM).

Stand Age (Years)	Mean Tree Height (m)	Trees/Plot (no/900 m2)	Mean dbh (cm)	Basal Area (m2 ha−1)	LAI Method
					HPM	LFM	CMM
13	1.93	99	1.98	0.59	0.90	1.47	0.66
40	5.46	300	7.46	22.87	1.48	2.66	1.77
72	9.28	62	20.72	32.92	1.45	3.88	2.13
92	12.87	21	43.36	35.79	1.47	3.36	1.79

) indicated an increasing average tree size with time since fire, with a decreasing stem density after a peak in the 40 year-old stand. Overall, LAI followed an increasing trend with time since fire, peaking at the 72 year-old plot, and slightly decreasing in the 90 year-old plot. As the plots were chosen to be undisturbed after the fire event, the slightly decrease in stem density in the 90 year-old plot was due to the self-thinning process. The LFM estimated higher values for LAI compared to HPM and CMM.

3.2. Foliage Properties and Their Seasonal Variation

Outputs from the mixed effect model indicated that the average values of both LMA and the photosynthetic MM model parameters were not different between the study plots, except for higher A_sat in the 13 year-old plot, suggesting the potential for using a common (average) light response curve across all study sites (

Table 4. Plot-specific fixed effect estimates (±1 standard error) for Leaf Mass per Area (LMA, g m⁻²), saturated net photosynthetic rate (A_sat, μmole m⁻² s⁻¹), half saturation coefficient (K_m, μmole quanta m⁻² s⁻¹) and dark respiration rate (R_d, μmole m⁻² s⁻¹). AIC_plot indicates the AIC for a simple linear regression with plot, while AIC_plot+month indicates the AIC when seasonal variation was taken into account by including a monthly random component. ΔAIC indicates the difference between AIC_plot+month and AIC_plot.

Parameter	13 Year-Old	40 Year-Old	72 Year-Old	92 Year-Old	AICplot	AICplot+month	ΔAIC
LMA	139.3 (±4.5)	146.4 (±4.4)	142.5 (±4.4)	139.0 (±3.3)	2253	2231	−22
Asat	10.9 (±1.3)	5.9 (±1.3)	8 (±1.4)	7.9 (±1.0)	1395	1207	−188
Km	303.4 (±37.4)	223.3 (±36.9)	235.6 (±38.2)	239.4 (±27.9)	3235	3121	−114
Rd	0.53 (±0.1)	0.48 (±0.1)	0.45 (±0.1)	0.49 (±0.07)	158	111	−47

). However, significant variation in the random (month) component was also identified (Table 4), suggesting that seasonal variation should be considered (monthly photosynthetic light response model setup).

Among all foliar properties studied, LMA showed the lowest seasonal variation (Figure 2A–D, Table A3). Monthly A_sat variability followed the seasonal patterns of temperature and precipitation variation (Figure 2E–H) with the highest values observed in spring. Between plots the highest seasonal variation was found in the youngest and less dense plot (Table A3). For all plots, lower A_sat values were observed in the summer of 2020 (Figure 2E–H), while K_m covaried with A_sat (Figure 2I–L). Leaf dark respiration showed the highest seasonal variation, with middle age stands having the higher CV.

The saturated net photosynthetic rate and the dark respiration were not statistically different between sunlit and shaded needles (A_sat,sun = A_sat,sh, R_d,sun = R_d,sh) as inferred from the Wilcoxon tests applied within and across all plots (Table A4). However, the half saturation coefficient was different between sunlit (K_m,sun) and shaded needles (K_m,sh). In order to simplify our model parameterisation, a correction for K_m,sh was implemented by multiplying the value of K_m,sun by 0.76, estimated from the ratio of K_m,sh/K_m,sun across all monthly estimates.

3.3. GPP Simulations

Simulated monthly GPP across the four study plots and the twelve model setups is presented in Figure 3. All model setups illustrated a characteristic peak of GPP during the growing season with the highest values in June, following seasonal variation in PAR flux and temperature. In all model setups, simulated annual GPP increased with stand age (Figure 4) until maturity (72 year-old) and then slightly decreased, except for the HPM model setups, where GPP increased continuously with post-fire time.

In general, the BL model simulated lower monthly GPP fluxes compared to the SS model, yielding a lower annual GPP across all study plots (

Table 5. Annual GPP (g C m⁻² year⁻¹) computed from the sum of simulated daily GPP for the period of 07/2019 to 06/2020. For each site, the estimated long-term (2000–2016) GPP (±one SD) from an independent dataset [66] is also provided as a baseline. For the 40 year-old and the 72 year-old plots the independent dataset provided the same GPP, as these two plots are within a 500 m radius.

GPP Model	LAI Method	Photo Capacity	13 Year-Old	40 Year-Old	72 Year-Old	92 Year-Old
Big Leaf Model	HPM	variant	630	724	840	895
	LFM	variant	830	910	1172	1204
	CMM	variant	511	785	996	982
	HPM	average	616	674	751	778
	LFM	average	816	846	1043	1042
	CMM	average	499	730	890	852
Sun/Shade Model	HPM	variant	767	945	1075	1125
	LFM	variant	999	1255	1670	1667
	CMM	variant	658	1036	1315	1252
	HPM	average	779	903	987	1012
	LFM	average	1004	1195	1519	1483
	CMM	average	673	989	1205	1124
Zhang et al. (2017)	-	-	816 (±98.7)	1221 (±126.5)	1221 (±126.5)	1060 (±102.1)

). Compared with the independent GPP dataset predictions, the BL model estimates were within the mean ± one SD range 5 out of 24 times and the SS model estimates were within the baseline range 10 out of 24 times.

In terms of the variation related to the method used to estimate LAI, the LFM model setups yielded the highest GPP, as expected from the higher LAI input values. Under the LFM parameterisation the simulated annual GPP was 28% and 38% higher compared to the CMM and HPM methods, respectively. Compared with the independent GPP dataset predictions, the CMM model setup annual estimates were within the mean ± one SD range 4 out of 16 times, the HPM 5 out of 16 and the LFM 6 out of 16.

Models that run with the average photosynthetic capacity parameters illustrated a lower variability in simulated monthly GPP compared to the variant model setups (CVa = 42.3% vs. CVv = 52.4%, D_AD = 13.60, p < 0.001). During most of the simulation period the variant model setups simulated higher GPP than the average photosynthetic capacity models, except for the period between 07/2020 and 09/2020 (Figure 3). Compared with the independent GPP dataset predictions, the average photosynthetic capacity models were within the baseline range 7 out of 24 times, and the variant photosynthetic capacity models were so 8 out of 24 times.

4. Discussion

In this study we provided a semi-empirical framework to quantify GPP variation across Mediterranean pine forests of different post-fire age. In accordance with forest dynamics theory, we simulated an increase in GPP with stand age until a maturity stage, followed by a small decrease in older stands. The shape of the GPP curve with time did not significantly change with the method used to upscale CO₂ fluxes from the leaf to the stand level. The Big Leaf model setups systematically underestimated GPP compared to the Sun/Shade model setups. Differences in the method used to parametrise our framework yielded an average plot GPP between 844 and 1166 gC m⁻² y⁻¹ across different LAI estimation methods, and an average GPP range from 934 to 1010 gC m⁻² y⁻¹ between constant and variant photosynthetic capacity setups. Overall, our simulations are within the range of GPP suggested for Mediterranean pine forests, namely from 816 to 1600 gC m⁻² y⁻¹ [48,66,84].

The trajectory of gross primary productivity with time since disturbance in unmanaged forests is frequently reported to present a parabolic shape [7] following the development of stand structure and LAI. After a wildfire, GPP immediately decreases, due to reduction in the available photosynthetic tissue biomass. Recovery of GPP can begin from the first year after the fire [85] following a steadily increasing trajectory. The degree of GPP rebound during the first few years is affected by the severity of the disturbance [14,86] as well as by water and nutrient availability [85,87]. A decade after the fire, when newly established trees take over, strong associations between GPP and canopy structure have been observed [88]. From this stage onwards, the recovery of GPP can be mostly explained by the recovery of LAI [89]. Across our post-fire chronosequence, LAI increased from younger to mature stands and then decreased at the 92 year-old stand (Table 3), probably due to competition between trees for light, water and nutrients [20,90]. Optical methods, such as CMM and HPM, generally underestimate LAI values when compared to direct and/or semi-direct methods, such as LFM, especially at forest stands [17,91]. The trend of LAI with time since the last stand-replacing fire controlled variation in GPP and shaped the typical parabolic curve with a peak in productivity at mature stands (Figure 3). This GPP trend was found in CMM and LFM model parameterisations, where stand-level LAI estimates were higher. In the HPM parameterisation, the lower LAI field estimates used to force the models yielded a continuously increasing GPP with stand age, as a result of carbon assimilation rate being away from the asymptote expected for relatively dense forest canopies with high LAI values. The results from the comparison of the variation related to the method used to estimate LAI illustrated that LAI measured using the LFM has the best agreement with the independent GPP dataset predictions. Overall, the selection of the LAI parameterisation method can have important effects on the simulation of GPP in studies that couple biometric and gas exchange data to estimate forest productivity, as has also been highlighted in satellite-data driven GPP simulation studies [92].

The control of climate, canopy structure and leaf physiology in controlling GPP has gained a lot of attention over the last decades. Many studies have explored the effect of climate [93,94], leaf traits [95,96] and canopy structure [97,98] in determining GPP across different ecosystem types. For many ecosystems, GPP increases until a peak during the growing period when solar radiation, temperature and water availability are optimum and then gradually decreases, with the contribution of the individual climate factors shifting with latitude [99,100]. This pattern has also been observed for Mediterranean forests, with GPP peaking during solar radiation and temperature optima [28,101], although for pine forest in semi-arid conditions GPP seems to maximise earlier in the year, representing a shift from radiation-controlled to moisture-controlled conditions [102]. Overall, Mediterranean forests seem to achieve a lower maximum and annual GPP, compared to tropical, subtropical and temperate forests [100]. However, the relative control of climate, canopy structure and leaf physiology on forest GPP has been studied less within similar ecosystem types (see, however, [103,104]). Among our study plots, found within relatively short distances from each other (max distance 28.2 km), similar climatic conditions prevail (Figure A1). Furthermore, no differences in the photosynthetic capacity parameters have been identified between the study sites (Table 4) and thus variability in stand and canopy structure can be considered as the main source of simulated GPP variation among the sites. Across our study sites, canopy structure variability is related to the time since the last stand-replacing fire and thus the four study plots can be used to represent the expected GPP pattern with post-fire time. In accordance with ecosystem theory, GPP seems to peak at mature stands, and thus maintenance of this forest state can be important for carbon management purposes [7].

In our simulations, seasonal variation of GPP (Figure 3) seems to be mainly driven by seasonal patterns of PAR [100]. Simulation of GPP based on satellite-driven models is frequently made through light use efficiency models [105], where primary productivity is estimated as the product of incoming PAR, the fraction of absorbed PAR and various environmental scalars (functions that do or do not limit GPP based on prevailing weather conditions and nutrient availability). These models frequently assume a constant (throughout the year) maximum efficiency of the foliage to convert the absorbed energy to carbon (like our average photosynthetic capacity models), which is being reduced by the environmental scalars [66]. Field observations suggest, however, that the photosynthetic capacity of leaves varies substantially throughout the year [84,106], as our results also show (Figure 2E–P and Table 4), and thus the use of constant photosynthetic parameter values may be inappropriate when simulating GPP [84]. In our varying photosynthetic capacity model parameterisation, we used monthly estimates of the parameters of the photosynthetic light response curve to take into account non-stomatal limitation shifts in photosynthetic capacity, such as reduced mesophyll conductance and photochemical and enzymatic limitations [107]. We further assumed that short-term effects of factors such as temperature and water availability would be less important in determining daily carbon assimilation rates, and thus avoided the use of scalars or functions that regulate canopy conductance, for example through the implementation of stomatal conductance regulation algorithms [108]. Incorporating such short-term responses of photosynthesis to weather variability would probably increase the accuracy of our simulations and reduce simulated GPP fluxes. However, in this study our purpose was to describe a simplified scheme for upscaling measurements of leaf gas-exchange to the stand level, by avoiding the use of detailed physiological algorithms, and assuming that once variation in some key foliage properties has been taken into account, patterns of incoming solar radiation might be the key determinants of forest primary productivity [52]. We however acknowledge that our simplified framework depends strongly on the availability of monthly photosynthetic data that require a significant time spent in the field.

During our study period, a strong reduction in the photosynthetic rate in July 2020 has been observed, with A_sat = 5.365 (μmole m⁻²s⁻¹) being 22.3% lower than the respective July 2019 A_sat. Precipitation during June and July 2020 (36.7 mm) was more than three times larger than that of June and July 2019 (11 mm), suggesting that short-term effects of photosynthetic stomatal limitation might not be as important as other process such as the drought memory phenomenon [11]. In our case this would mean that in 2019, trees were adapted to dry conditions from the previous dry months, therefore creating the drought memory [109], reflecting a higher photosynthetic rate in July 2019 in contrast to July 2020, when trees did not have the pre-adaption conditions. Experiments have shown that pre-stressed plants maintain higher photosynthetic rates than non-pre-stressed plants at a drought event, even if both show similar stomatal conductance and respiration [110,111], supporting the use of seasonal shifts in photosynthetic capacity rather than regulations of stomatal conductance based on short-term weather variation.

The 13 year-old plot, which is the youngest in our plot network, illustrated the lower stand basal area and LAI. Due to small tree size the stand has not yet achieved canopy closure and is characterised by the most heterogenous microclimatic conditions. The low LAI in the 13 year-old stand results in a relatively higher percentage of tree foliage being sunlit compared to other plots (13 year-old, 56.6%, 40 year-old, 40.0%, 72 year-old, 35.9% and 92 year-old, 38.2%) [112], leading to higher A_sat due to photo-acclimation. Murchie and Horton [113] supported the view that light-demanding tree species (such as P. brutia) can acclimate by changing their leaf chlorophyll content and subsequently their photosynthetic capacity. Our data revealed that photosynthetic capacity illustrated a higher seasonal variation in plots with lower LAI and canopy closure. In such stands, parts of the canopy are sunlit for longer, or if shaded, receive more radiation than denser stands. This could lead to morphological and physiological adaptation of leaves to optimise the rate of photosynthesis [114,115], as has also been found across our study sites where needles at younger stands are optimized along the LMA–photosynthetic capacity spectrum [116]. In addition, the heterogenous within-stand microclimatic conditions also probably lead to the increased variability that has been observed in gas exchange (13 year-old A_sat and K_m CV vs. other plots, Table A3), because of limited buffering of abiotic variability from neighbouring trees. A negative trend in A_sat variability with LAI was also observed across plots, further supporting the potential effects of canopy structure-dependent photo-acclimation. Thus, the effects of canopy structure on GPP variation can be both direct (LAI amount of available photosynthetic tissue) and indirect (buffering effects on photosynthetic variability).

5. Conclusions

In this study we presented a simplified framework to estimate GPP in Mediterranean pine forests of different post-fire age. Overall, our simulations seem to be in agreement with GPP estimates from eddy-flux [84], process-based [101,117] and satellite-driven [66] models. Compared to the GPP values reported for temperate conifer forests [118], Mediterranean pine forest are probably achieving relatively lower C sequestration, due to water limitation during the dry summer period [119,120]. Mediterranean forests are currently positively contributing to the global carbon balance [121], but with the expected increase in fire frequency and severity, estimates of carbon sequestration among different stand ages are crucial to quantify their climate mitigation potential and design appropriate management strategies [118].