Sun-induced chlorophyll fluorescence and photochemical reflectance index improve remote-sensing gross primary production estimates under varying nutrient availability in a typical Mediterranean savanna ecosystem

This study investigates the performances of different optical indices to estimate gross primary production (GPP) of herbaceous stratum in a Mediterranean savanna with different nitrogen (N) and phosphorous (P) availability. Sun-induced chlorophyll fluorescence yield computed at 760 nm (Fy760), scaled photochemical reflectance index (sPRI), MERIS terrestrial-chlorophyll index (MTCI) and normalized difference vegetation index (NDVI) were computed from near-surface field spectroscopy measurements collected using high spectral resolution spectrometers covering the visible near-infrared regions. GPP was measured using canopy chambers on the same locations sampled by the spectrometers. We tested whether light-use efficiency (LUE) models driven by remote-sensing quantities (RSMs) can better track changes in GPP caused by nutrient supplies compared to those driven exclusively by meteorological data (MM). Particularly, we compared the performances of different RSM formulations – relying on the use of Fy760 or sPRI as a proxy for LUE and NDVI or MTCI as a fraction of absorbed photosynthetically active radiation (fAPAR) – with those of classical MM. Results showed higher GPP in the N-fertilized experimental plots during the growing period. These differences in GPP disappeared in the drying period when senescence effects masked out potential differences due to plant N content. Consequently, although MTCI was closely related to the mean of plant N content across treatments (r = 0.86, p < 0.01), it was poorly related to GPP (r = 0.45, p < 0.05). On the contrary sPRI and Fy760 correlated well with GPP during the whole measurement period. Results revealed that the relationship between GPP and Fy760 is not unique across treatments, but it is affected by N availability. Results from a cross-validation analysis showed that MM (AICcv = 127, MEcv = 0.879) outperformed RSM (AICcv = 140, MEcv = 0.8737) when soil moisture was used to constrain the seasonal dynamic of LUE. However, residual analyses demonstrated that GPP predictions with MM are inaccurate whenever no climatic variable explicitly reveals nutrient-related changes in the LUE parameter. These results suggest that RSM is a valuable means to diagnose nutrient-induced effects on the photosynthetic activity. Published by Copernicus Publications on behalf of the European Geosciences Union. 6352 O. Perez-Priego et al.: Remote-sensing-based model of photosynthesis


Introduction
Human-induced nutrient imbalances are affecting essential processes that lead to important changes in ecosystem structure and functioning (Peñuelas et al., 2013). In spite of the crucial role of nutrients in regulating plant processes, efforts to describe and predict the response of photosynthesis to such changes with remote-sensing information have been limited. In the framework of the classical Monteith light-use efficiency (LUE) model (Monteith, 1972), estimates of photosynthesis (hereafter gross primary productivity, GPP) are based on three key quantities: (i) the fraction of photosynthetically active radiation (f APAR) absorbed by the vegetation; (ii) potential LUE (or maximum, LUE m ), normally taken from lookup tables and associated with plant functional types (Heinsch et al., 2006); and (iii) correction factors related to meteorological conditions that limit LUE m . Although nitrogen (N) deficiencies have been recognized as one of the main correction factors of LUE m (Madani et al., 2014), the predictive capability of LUE models is usually circumspect as they operate based on the general assumption that plants are under non-limiting nutrient conditions.
Very little attention has been given to nutrient-induced effects on f APAR and LUE in common formulations of LUE models. Light absorption by plants is given by chlorophyll pigments that enable photosynthetic processes. Assuming a correlation between leaf chlorophyll pigments and leaf N content, noting that N atoms are basic components of the chlorophyll molecular structure, several studies have demonstrated that leaf N content can be estimated through chlorophyll-related hyperspectral vegetation indices (Baret et al., 2007;Schlemmer et al., 2013). Among these indices, the MERIS (Medium Resolution Imaging Spectrometer) terrestrial-chlorophyll index (MTCI; Dash and Curran, 2004) has been used as a proxy for f APAR (Rossini et al., 2010;Wang et al., 2012). However, leaf N content is a functional trait that controls GPP not only because it scales with chlorophyll but also because it regulates enzyme kinetic processes driving photosynthesis and hence the physiological status of the plant (Huang et al., 2004;Walker et al., 2014). Then, prescribing biome-specific LUE parameters and correcting LUE m only for climatic and environmental conditions may hamper the accurate prediction of GPP (Yuan et al., 2014). For these reasons, recent literature has called for better physiological descriptors of the dynamic behavior of LUE .
The sun-induced chlorophyll fluorescence (SIF) or physiological-related reflectance indices such as the photochemical reflectance index (PRI) provide a new optical means to spatially infer LUE (Damm et al., 2010;Guanter et al., 2014;Rossini et al., 2015) and can provide diagnostic information regarding plant nutrient and water status Pérez-Priego et al., 2005;Suárez et al., 2008;Tremblay et al., 2012). From a physiological perspective, the efficiency of green plants to transform absorbed light into chemical energy during photosynthesis can be characterized by two main photoprotective mechanisms: (i) nonphotochemical quenching that can be detected using the Photochemical Reflectance Index (PRI), originally proposed by (Gamon et al., 1992) to track changes in the de-epoxidation state of the xanthophyll cycle pigments, and (ii) chlorophyll fluorescence, the dissipation of energy that exceeds photosynthetic demand (Krause and Weis, 1984). The PRI has been directly correlated with LUE Gamon et al., 1997;Nichol et al., 2000;Peñuelas et al., 2011;Rahman et al., 2004). However, such a relation may vary because of the sensitivity of the PRI to confounding factors like those associated with temporal changes in the relative fraction of chlorophyll : carotenoid pigment composition (Filella et al., 2009;Porcar-Castell et al., 2012), viewing angles and vegetation structure Grace et al., 2007;Hall et al., 2008;Hilker et al., 2008).
Alternatively, the estimation of SIF by passive remotesensing systems has been proven feasible in recent years from satellites Lee et al., 2013;Parazoo et al., 2014) to the field (Damm et al., 2010;Guanter et al., 2013;Meroni et al., 2011), and opens further possibilities to directly track the dynamics of LUE (Damm et al., 2010;Guanter et al., 2014). Although SIF correlates with LUE, such relations might not be conservative since chlorophyll fluorescence emission varies among species types (Campbell et al., 2008) or with stress conditions such as nutrient deficiencies (Huang et al., 2004;Mc-Murtrey et al., 2003) or drought (Flexas et al., 2002;Pérez-Priego et al., 2005). Likewise with the PRI, the retrieval of SIF from the apparent reflectance signal is not trivial as long as it is affected by the vegetation structure or canopy background components (Zarco-Tejada et al., 2013).
Comparable spatial and temporal resolutions of radiometric and ground-based GPP measurements are essential to accurately optimize LUE model parameters, particularly in heterogeneous ecosystems. Previous studies have related ecosystem-scale eddy covariance fluxes to radiometric measurements taken in single points to constrain LUE models. However, the explanatory power of LUE models might be greatly reduced by the spatial mismatch between radiometric and eddy covariance flux footprints (Gelybó et al., 2013;Porcar-Castell et al., 2015). Similar issues occur in smallscale factorial experiments where comparable measurements on an intermediate scale between leaf-scale cuvette measurements and ecosystem-scale eddy covariance measurements are required. Here, we tried to overcome such limitations by combining ground-based radiometric and CO 2 flux measurements with a similar extension of the measurement footprint using portable spectrometers and canopy chambers in a nutrient-manipulation experiment.
The main objective of this study was to evaluate whether traditional LUE models driven by meteorological and phenological data (MM) entail a limited assessment of the environmental controls on GPP. More particularly, we evaluated  (Table 1). Overall, leaf area measurements of the herbaceous stratum characterized the growing season phenology as peaking early in April and achieving senescence by the end of May (Table 1). The experiment consisted of four randomized blocks of about 20 m × 20 m. Each block was separated into four plots of 9 m × 9 m with a buffer of 2 m in between to avoid bound-  Table 1. Ancillary data resulting from the analysis. Green plant area index (PAI g ), fraction of PAI in different plant forms (f PAI ), and C, N, and P plant content. The N : P ratio also is shown. Data correspond to the mean value and standard deviation (SD) of the subsamples taken in each plot and treatment.  ary effects. In each block, four treatments were applied (see Fig. 1): 1. control treatment (C) with no fertilization 2. nitrogen-addition treatment (+N) with an application of 100 kg N ha −1 as potassium nitrate (KNO 3 ) and ammonium nitrate (NH 4 NO 3 ) 3. phosphorous-addition treatment (+P) with an application of 50 kg P ha −1 as monopotassium phosphate (KH 2 PO 4 ) and 4. N-and P-addition treatment (+NP), juxtaposing treatments 1 and 2.
Each fertilizer was dissolved in water and sprayed on foliage early in the growing season (21 March 2014). The same amount of water used in the fertilizer solutions (∼ 2 L m −2 ) was sprayed on the C treatment to avoid water imbalances among treatments. Within each plot, two permanent, non-disturbed parcels (32 in total, see black squares in Fig. 1) were dedicated to monitoring CO 2 fluxes (net ecosystem CO 2 exchange, NEE; and daytime ecosystem respiration, R eco ). While NEE measurements were performed over the course of the day (from early in the morning to late afternoon), spectral measurements were conducted simultaneously with flux measurements only around noon on half of the parcels (16 in total).
Flux and spectral measurements were carried out in four field campaigns: -campaign no. 1: before fertilization (20 March 2014), -campaign no. 2: 3 weeks after fertilization (15 April 2014) during the peak of the growing period, -campaigns no. 3 and no. 4: on 7 and 27 May 2014, respectively, concurring with the drying period were performed to evaluate joint effects related to physiological senescence processes.
Ancillary measurements were taken in every field campaign as follows: green plant area index (PAI g ) and aboveground biomass were directly measured by harvest in four parcels (0.25 m × 0.25 m) within each plot in the surrounding area where spectral and flux measurements were taken. All samples were refrigerated just after collection and transported for laboratory analyses. Fresh samples were separated into functional groups; the sample was scanned and green plant area was measured using image analysis (WinRHIZO, Regent Instruments Inc., Canada). Afterwards, fresh samples were dried in an oven at 65 • C for 48 h and weighed to determine dry biomass. To analyze the nutrient content in leaf mass, biomass subsamples were ground in a ball mill (RETSCH MM200, Retsch, Haan, Germany) and total C and N concentrations were determined with an elemental analyzer (Vario EL, Elementar, Hanau, Germany). P concentrations were also measured: 100 mg biomass subsamples were diluted in 3 mL of HNO 3 65 %, (Merck, Darmstadt, Germany) and microwave digested at high pressure (Multiwave, Anton Paar, Graz, Austria; Raessler et al., 2005). Afterwards, elemental analysis was conducted using inductively coupled plasmaoptical emission spectrometry (ICP-OES, Optima 3300 DV, Perkin Elmer, Norwalk, USA).

Flux measurements and meteorological data
Net CO 2 fluxes were measured with three transparent chambers of a closed dynamic system. The chambers consisted of a cubic (0.6 m × 0.6 m × 0.6 m) transparent low-density polyethylene structure connected to an infrared gas analyzer (IRGA LI-840, Lincoln, NE, USA), which measures CO 2 and water vapor mole fractions (W ) at 1 Hz. The chambers were equipped with different sensors to acquire environmental and soil variables, all installed at the chamber ceiling: photosynthetically active radiation (PAR) was measured with a quantum sensor (Li-190, Li-Cor, Lincoln, NE, USA) placed outside of the chamber to be handled and leveled; air and vegetation temperatures were measured with a thermistor probe (T a , type 107, Campbell Scientific, Logan, Utah, USA) and an infrared thermometer (T c , IRTS-P, Apogee, UT, USA); atmospheric pressure (P ) was measured inside the chamber using a barometric pressure sensor (CS100, Campbell Scientific, Logan, Utah, USA). The chambers were also equipped with soil temperature and humidity sensors; soil water content was determined with an impedance soil moisture probe (Theta Probe ML2x, Delta-T Devices, Cambridge, UK) at 5 cm depth and soil temperature (type 107, Campbell Scientific, Logan, Utah, USA) at 10 cm depth. Saturation vapor at surface temperature (i.e., air-to-leaf vapor pressure deficit, VPD) was computed using T c and relative humidity, which was derived from water vapor molar fraction measured with the IRGA (Perez-Priego et al., 2015). The chamber operated as a closed dynamic system. A small pump circulates an air flow of 1 L min −1 through the sample circuit: air is drawn from inside the chamber -through three porous hanging tubes spatially distributed through the chamber headspace -to the infrared gas analyzer; this air flow is then returned to the chamber. The hanging tubes allowed spatially distributed sampling, obviating the need to homogenize air during chamber deployment. Nevertheless, one small fan (12 V, 0.14 A) was fixed at 0.3 m on a floor corner of the chamber and angled 45 • upward.
A 0.6 × 0.6 m metal collar was installed in each permanent parcel of each plot. The collar provided a flat surface onto which the bottom of the chamber was placed. The chamber was open and ventilated for 1 min prior to measurement so that initial air composition and temperature in the confined environment of the chamber represented natural atmospheric conditions (as much NEE as Reco). For the NEE measurement, the transparent chamber was placed on the collar (closed position, lasted 3 min as a general rule), and fluxes were calculated from the rate of change in the CO 2 molar fraction (referenced to dry air) within the chamber. A similar procedure was carried out for R eco but using an opaque blanket that covered the entire chamber and kept it dark during the measurements (PAR values around 0). Fluxes were calculated according to Pérez-Priego et al. (2015).
Shortly, the flux calculation algorithm reduces flux uncertainties (i.e., NEE and R eco ) by including the change-point detection method to determine the stabilization time, which defines the initial slope of the regressions, and a bootstrap resampling-based method to improve confidence in regression parameters and to optimize the number of data points used for flux calculation. In addition, a statistical analysis of residuals was performed to automatically detect the best fit among alternative regressions (i.e., quadratic, hyperbolic tangent saturating function, exponential, linear). These analyses were implemented in a self-developed R Package (available upon authors request or at the following link http: //r-forge.r-project.org/projects/respchamberproc/). NEE and R eco measurements were taken over the course of the day (from sunrise to sunset) for each field campaign. Chamber disturbance effects and correction for systematic and random errors (i.e., leakage, water dilution and gas density correction, and light attenuation by the chamber wall) were applied according to Perez-Priego et al. (2015).

Field spectral measurements
Midday spectral measurements at canopy level were carried out under clear-sky conditions using two portable spectrometers (HR4000, OceanOptics, USA) characterized by different spectral resolutions. Spectrometer 1, characterized by a full width at half maximum (FWHM) of 0.1 nm and a 700-800 nm spectral range, was specifically designed for the estimation of sun-induced chlorophyll fluorescence at the O 2 -A band (760 nm). Spectrometer 2 (FWHM = 1 nm, 400-1000 nm spectral range) was used for the computation of reflectance and vegetation indices. Spectrometers were housed in a thermally regulated Peltier box, keeping the internal temperature at 25 • C in order to reduce dark current drift. The spectrometers were spectrally calibrated with a source of known characteristics (CAL-2000 mercury argon lamp, OceanOptics, USA), while the radiometric calibration was inferred from cross-calibration measurements performed with a calibrated FieldSpec FR Pro spectrometer (ASD, USA). This spectrometer was calibrated by the manufacturer with yearly frequency.
Incident solar irradiance was measured by nadir observations of a leveled calibrated standard reflectance panel (Spectralon, LabSphere, USA). Measurements were acquired using bare fiber optics with an angular field of view of 25 • . The average canopy plane was observed from nadir at a distance of 110 cm (43 cm diameter field of view) allowing for collecting measurements of 50 % of the surface area covered by the chamber measurements. The manual rotation of a mast mounted horizontally on the tripod allowed sequential observation of the vegetated target and the white reference calibrated panel. More in detail, every acquisition session consisted in the consecutive collection of the following spectra: instrument dark current, radiance of the white reference panel, canopy radiance and radiance of the white reference panel. The radiance of the reference panel at the time of the canopy measurement was then estimated by linear interpolation.
For every acquisition, 3 and 10 scans (for spectrometers 1 and 2, respectively) were averaged and stored as a single file. Five measurements were collected for each plot. Spectral data were acquired with dedicated software (Meroni and Colombo, 2009) and processed with a specifically developed IDL (ITTVIS IDL 7.1.1) application. This application allowed the basic processing steps of raw data necessary for the computation of the hemispherical conical reflectance factor described by Meroni et al. (2011).
The following indices were selected as suitable to investigate long-term nutrient-mediated effects on photosynthesis. The normalized difference vegetation index (NDVI; Rouse et al., 1974) was selected because it correlates well with plant area and among traditional spectral vegetation indices is used worldwide by classical LUE models as a surrogate for f APAR (Di Bella et al., 2004). The MTCI (Dash and Curran, 2004) was selected because it was specifically designed for canopy chlorophyll content estimation, and recently used as a proxy for f APAR as well as NDVI. In this study we used the PRI and SIF as surrogates for LUE. A scaled PRI (sPRI) calculated as (PRI+1)/2 was used. SIF was estimated by exploiting the spectral fitting method described in Meroni et al. (2010), assuming linear variation of the reflectance and fluorescence in the O 2 -A absorption band region. The spectral interval used for SIF estimation was set to 759.00-767.76 nm for a total of 439 spectral channels used. For methodological distinction among existing approaches, hereafter SIF is referred to as F760. Because F760 is affected by PAR we use the apparent chlorophyll fluorescence yield (Fy760; Rossini et al., 2010) computed as the ratio between F760 and the incident radiance in a nearby spectral region. A summary of the formulation to compute the vegetation indices and their corresponding target and proxy in the LUE model approach are presented in Table 2.

Relationship between GPP and remote-sensing data
Ecosystem-level GPP was computed as the difference between NEE and daytime R eco taken consecutively with the chambers. To assess how GPP is modulated by light among treatments and over the phenological cycle of the herbaceous stratum, we computed the parameters of photosynthetic light response curve (PLRC). Specifically, the Michaelis-Menten function was fitted to GPP and PAR data taken throughout the course of the day (from sunrise until sunset) for each field campaign and treatment as follows: where α is a parameter describing the photosynthetic quantum yield (µmol CO 2 µmol photons −1 ), and β is the parameter that extrapolates to GPP at saturating light condition (µmol CO 2 m −2 s −1 ). According to Ruimy et al. (1994), we used the optimized parameters of the PLRC as defined in Eq.
(1) to estimate the GPP at 2000 µmol quantum m −2 s −1 of PAR (hereafter referred to as GPP 2000 ). We evaluated direct relationships between those GPP measurements taken around noon (between 11:00 and 15:00 solar time) with the chamber (GPP noon ) and sequentially measurements of Fy760 and spectral indices (NDVI, sPRI, MTCI). In addition, to avoid confounding factors in the relationship between Fy760 and sPRI and photosynthesis, we also used GPP 2000 as a maximum photosynthetic capacity descriptor.

Monteith's light-use efficiency modeling approaches
Following Monteith's LUE framework (Eq. 2) two alternative modeling approaches were used:

Statistical analysis and model performance
All model formulations were optimized using GPP noon and spectral measurements taken at midday. Since the means of spectral measurements per treatment could have unequal variance, a Welch's t test was performed to evaluate significant differences between the mean values of the different vegetation indices for each treatment and over the four field campaigns. In addition, an analysis of covariance (AN-COVA) was used to test whether or not there was a significant interaction by the treatment effect between GPP noon and Fy760 and different spectral indices. Like vegetation indices, a t test was performed to the daily average of GPP taken over the course of the day (GPP daily ).

Cross-validation analyses and model evaluation
Different model formulations were evaluated in leave-oneout (LOO) cross-validation: from the whole data set com-  Meroni and Colombo (2009) posed by n observations, one data point at a time was removed. The model was fitted against the n−1 remaining data points (training set), while the excluded data (validation set) were used for model evaluation. The cross-validation process was then repeated n times, with each of the n observations used exactly once as the validation set. For each validation set of the cross-validated model, statistics were calculated. Model accuracy was evaluated by means of different statistics according to Janssen and Heuberger (1995): root mean square error (RMSE), relative root mean square error (rRMSE), determination coefficient (r 2 ), and model efficiency (ME). The model performances in LOO crossvalidation were also calculated and reported as RMSE cv , rRMSE cv , r 2 cv , and ME cv . The Akaike information criterion (AIC cv ) was used to evaluate the trade-off between model complexity (i.e., number of parameters) and explanatory power (i.e., goodness of fit) of the different model formulations proposed. The AIC cv is a method based on information theory that is useful for statistical and empirical model selection purposes (Akaike, 1998). Following Anderson et al. (2000), in this analysis we used the following definition of AIC cv : where n is the number of samples (i.e., observations), p is the number of model parameters, and RSS cv is the residual sum of squares divided by n.
The LUE model formulations proposed in Sect. 2.4 can be ranked according to AIC cv , where the model with lowest AIC cv is considered the best among the different model formulations.
All model parameters (MM, and RSM) were estimated by using a Gauss-Newton nonlinear least-square optimization method (Bates and Watts, 2008), and standard errors of parameters were estimated by bootstrapping (number of sampling, n = 500; Efron and Tibshirani, 1994), both implemented in the R standard package (R version 3.0.2, R Development Core Team, 2011).

Effects of fertilization on plant nutrient contents and GPP
Fertilization caused strong variations in leaf N and P content among treatments, plant forms and across field campaigns (Table 2), while total N content in plants ranged slightly between 13.8 ± 1.2 and 15.4 ± 1.7 mg g −1 for the C and +P treatments over the whole experiment. The largest increases in total N content were found in the peak of the growing season (no. 2, 20 March 2014), when +NP and +N treatments reached values of up to 23.7 ± 2.0 and 23.5 ± 4.1 mg g −1 , respectively. Although slightly lower, the differences in total N content between C and +P, and +NP and +N remained high over the drying period. Total P content was higher in +NP and +P treatments after fertilization, as compared to +N and C treatments. Consequently, the N : P ratio at the first campaign after fertilization (no. 2) achieved values of up to 14.2, 6.6, 6, and 3.7, in +N, C, +NP, and +P treatments, respectively. Similar differences in N : P between treatments were also observed during the drying period (no. 3 and no. 4, Table 2). On the other hand, PAI g ranged from 0.4 m 2 m −2 in campaign no. 4 to up to 2.5 m 2 m −2 in campaign no. 2. No differences were found in PAI g among treatments since grazing apparently offset any potential difference in the green aboveground production. Regarding variations in the fraction of plant forms, no significant differences were found between treatments. Fertilization caused significant differences in the GPP daily (p < 0.05) between N-addition treatments (mean values of 19.62 ± 4.15 and 18.19 ± 5.67 µmol CO 2 m −2 s −1 for +N and +NP, respectively) and C and +P treatments (14.31 ± 5.39 and 14.40 ± 4.09 µmol CO 2 m −2 s −1 , respectively) in the peak of the growing season (campaign no. 2); a relative difference of 37 % in GPP daily values was found between +N and +NP and C treatments. During the drying period, however, GPP was substantially downregulated (campaign nos. 3 and 4), and no significant differences were found in GPP daily , regardless of differences in plant N content observed among treatments. The potential photosynthetic capacity GPP 2000 (Fig. 2) derived from PLRC was similar in the four treatments in the pretreatment period (campaign no. 1; Fig. 2a). GPP 2000 varied throughout the season and peaked in the campaign no. 2 (15 April) in all treatments. At this time PLRC of the +N and +NP treatments diverged clearly from no N-addition treatments (C and +P; Fig. 2b). GPP 2000 was higher in +N and +NP treatments (18.6 and 20.1µmol CO 2 m −2 s −1 , respectively) compared to C and +P treatments (14.9 and 15.4 µmol CO 2 m −2 s −1 , respectively). After campaign no. 2, when the soil layer at 5 cm depth dried out appreciably (volumetric water content achieved values of 3 % vol., data not shown), vegetation progressively senesced and GPP 2000 in turn was downregulated and converged to similar values in all treatments, regardless of the higher N content observed in +N and +NP treatments as compared with C and +P treatments (Table 1). During the dry season, GPP 2000 decreased in all treatments ranging between 5.6 and 8 µmol CO 2 m −2 s −1 and no differences among treatments were observed ( Fig. 2c and d). These results indicate that the senescence of the herbaceous stratum, which is regulated by water availability, strongly modulated the photosynthetic capacity of the vegetation over the season.

Effects of fertilization on remote-sensing data
Optical properties of the analyzed plots were similar during campaign no. 1, before the nutrient application. A pronounced seasonal time course was observed for both Ph (sPRI and Fy760) and structural indices (St; NDVI and MTCI) with maximum values during the second campaign. It is interesting to note that, while for St indices the maximum values were reached in +N plots, +NP plots showed maximum Ph values. Vegetation indices and Fy760 then decreased in the drying period (Fig. 3). As for GPP, differences between treatments were more evident during campaign no. 2 when C plots showed statistically lower values for all the indices considered, while only MTCI was able to detect significant differences between N-fertilized plots (+N and +NP). Furthermore significant differences in Fy760 and MTCI between C and the other three treatments were found (p < 0.05) in the drying period (campaign no. 4). NDVI varied significantly with changes in PAI g with values of 0.4 in the campaign no. 4 up to 0.8 in the campaign no. 2 (p < 0.001, r 2 = 0.79).

Relationship between remote-sensing data and GPP
While Ph indices (Fy760 and sPRI) varied linearly with GPP noon in all treatments (p < 0.001, r 2 = 0.66 for Fy760 and p < 0.001, r 2 = 0.79 for sPRI, respectively; Fig. 4a and b), different patterns were observed for St: NDVI and GPP were best fitted by an exponential regression p < 0.001, r 2 = 0.77 (Fig. 4c), while a weak linear relationship between MTCI and GPP noon p < 0.05, r 2 = 0.45 (Fig. 4d) was found.
Although a weak relation between MTCI and GPP noon was found, MTCI was strongly correlated with plant N content (y = 14.17x − 2.49, p < 0.001, r 2 = 0.86). Note that these results are computed excluding data taken in the pretreatment campaign (no. 1) and differences in the relationship between remote-sensing data and GPP noon among treatments can be only attributed to nutrient-induced effects. The AN-COVA test did not show significant differences in slope or intercept of the relationship between GPP noon and sPRI, and NDVI across treatments. However, barely significant differences were found in the relationship between GPP noon and Fy760 (p < 0.1; Fig. 4b) and significant ones between GPP noon and MTCI (p < 0.01; Fig. 4d) between N-addition treatments (+N and +NP) and C treatments (C and +P).

Modeling GPP
Based on the AIC cv criterion, MM (VPD-SWC) outperformed MM-VPD, MM-SWC and RSMs. Although MM (VPD-SWC) showed high accuracy in the predictions (ME cv = 0.879, r 2 cv = 0.881), this model had a tendency to underestimate observation at high GPP noon values (see comparison between model predictions and observations, Fig. 6a-c). Note that the highest biases in modeled GPP noon values among MM models belong to +N and +NP treatments in field campaign no. 2. Since the four treatments experienced the same environmental conditions (i.e., comparable values of SWC, VPD, air temperature), this bias can be attributed to the higher N content (+N and +NP treatments) as compared to C and +P treatments. Remarkably, residuals of the MM (VPD-SWC) taken from periods with moist  Table 3. soil (SWC > 15) were significantly correlated with sPRI and Fy760 (p < 0.05; Fig. 7a and b, respectively). However, no biases between residuals and predictions were observed in RSM over the span of values and treatments (Fig. 8). Results from the evaluation of model performance indicated that RSM performs best when NDVI, rather than MTCI, is used as St in the Eq. (7) and, hence, as a proxy for f APAR (Table 3). Our results indicated that RSM performs best when either Ph (sPRI or Fy760) is combined with NDVI as St.

Effects of nutrients on GPP and remote-sensing data and their relationships
Nutrient fertilization, particularly N inputs, induced physiological changes manifested as an increase in photosynthetic capacity under high light conditions ( Fig. 2; Hirose and Werger 1994). As we expected, plant N content has shown to be a trait of photosynthesis that influences a variety of aspects of photosynthetic physiology (Ciompi et al., 1996;Sugiharto et al., 1990). These physiological changes were reflected on the optical properties, particularly on fluorescence and sPRI. The increase in fluorescence with N fertilization inputs was recently explained as the combined  effect that a higher N content has on (1) chlorophyll content, which magnifies APAR and enhances fluorescence signal, and on (2) the increased photosynthetic capacity that results in reduced non-photochemical quenching activity and consequently increases the fluorescence signal (Cendrero-Mateo et al., 2015). The relationships between GPP noon and Fy760 is not unique and may vary from optimal to nonoptimal environmental conditions (i.e., nutrient deficiencies, water stress), when other regulatory mechanisms might reduce the degree of coupling between fluorescence and photosynthesis (Cendrero-Mateo et al., 2015;Porcar-Castell et al., 2012). Although Fy760 was positively correlated with GPP noon , barely significant differences in the slope of this relationship were observed between treatments (Fig. 4b). Further studies are needed to fully explore the relationship between Fy760 and GPP noon under different stress conditions and over different ecosystems. However, if confirmed, the effect of nutrient availability on the relationship between Fy760 and GPP noon could have important implications in GPP modeling. This result suggests that the inclusion of a correction factor related to leaf N : P stoichiometry should be considered when modeling GPP assuming a linear relationship with fluorescence at plant functional type level Joiner et al., 2013).
In this study we also explored the capability of remote sensing to describe ecosystem functional properties defined as those quantities that summarize and integrate ecosystem processes and responses to environmental conditions and can be retrieved from ecosystem-level fluxes (e.g., GPP 2000 ) and structural measurements (Reichstein et al., 2014). GPP at light saturation (i.e., GPP 2000 ) is one example of an ecosystem functional property, shown here to be quite correlated to sPRI and Fy760 (Fig. 5). This result suggests that sPRI and Fy760 open also new opportunities for remote-sensing products to describe the spatiotemporal variability of essential descriptors of ecosystem functioning (Musavi et al., 2015). Inferring GPP 2000 using remote sensing has important implications both for monitoring global carbon cycle and for benchmarking terrestrial biosphere models.
MTCI was closely related to N content (r 2 = 0.86, p < 0.001), independent of other structural variables (i.e., PAI g ) and can be used as a good indicator of N availability. Although MTCI has been proven to be very sensitive to variations in chlorophyll contents (Dash and Curran, 2004) and hence linkable with light absorption processes, it was weakly correlated with GPP, particularly in plots added with N (+N and +NP; r 2 = 0.27, p < 0.01; Fig. 4d). A quite wide range of GPP noon values were found at high values of MTCI -high GPP noon values corresponding to the growing season and low ones to the drying period -which can be explained by two simultaneous mechanisms.
First, despite the high plant N content, physiological mechanisms including stomatal control or reduced carboxylation efficiency downregulate GPP (Huang et al., 2004) and ultimately might break the relationship between GPP noon and MTCI. Second, MTCI tracks changes in N content regardless of changes in canopy structure occurring during the dry season when grass achieved senescence (i.e., green to dry biomass ratio, PAI g ). More studies aimed at the separation of the combined effects of N and changes in green/dry biomass fractions on f APAR are essential. On the other hand, although NDVI followed the seasonal dynamic of PAI g , it saturated at high GPP noon values indicating the low ability of this index to detect spatial variations induced by N fertilization.
Although optical measurements were taken at high spatial resolution (< 0.36 m 2 ), the separation of confounding factors affecting sPRI or Fy760 is essential to elucidate the mechanistic association between sPRI or Fy760 and GPP. Like sPRI, the retrieval of Fy760 from the apparent reflectance signal can be also affected by vegetation structure or canopy background components (Zarco-Tejada et al., 2013). After optimization and selection of the best model parameters using NDVI and sPRI (or Fy760) as drivers, we analyzed the response of simulated GPP to variations in NDVI and sPRI ((or Fy760;Fig. 9). Results indicate that, at high GPP levels, it is Fy760 and sPRI but less NDVI that responded to GPP. However, at low GPP levels, either Fy760 or sPRI re-sponded to GPP on a small scale (Fig. 9b). Figure 9 suggests that the relationship between NDVI and sPRI or Fy760 is not unique and NDVI may play an important role in driving GPP in ecosystems characterized by marked seasonal variations. Our results highlight the complementarity between NDVI and Fy760 or sPRI. Particularly, NDVI assisted Fy760 or sPRI in predicting GPP under conditions with low biomass (i.e., low leaf area index), when confounding factors may affect Fy760 or sPRI. In semiarid ecosystems, the lack of sensitivity of sPRI or Fy760 to changes in GPP during dry conditions has been explained by the soil background effect on the reflectance signal (Barton and North, 2001;Mänd et al., 2010;Zarco-Tejada et al., 2013). Accordingly, Rahman et al. (2004) pointed out that conditions where sPRI performs best are in dense canopies with a low portion of bare soil.

Performances of different LUE modeling approaches
Here we aim at answering the question of how we can better simulate GPP using LUE modeling with varying nutrient availability and environmental conditions by drawing comparisons between the two model philosophies: RSM against MM approaches. There are an increasing number of studies focused on the development of LUE models driven by remotely sensed information to better explain spatiotemporal variations of GPP Rossini et al., 2012Rossini et al., , 2014. However, nutrient availability (and in particular N) greatly influences the spatial variability of LUE even within the same plant functional type (e.g., grasslands), and further studies are essential. The slightly better performance in cross-validation of the MM (VPD-SWC) against all model configurations, including RSM, supports the importance of a joint use of SWC and VPD as key parameters to constrain LUE in arid and semiarid ecosystems (Prince and Goward, 1995). However, residual analyses demonstrated that MM (VPD-SWC) was unable to track N-induced differences in GPP during the growing period, when both parameters are not limiting (Fig. 7). By contrast, accurate estimates of GPP were obtained with RSM both over the drying and the growing periods. These results also indicate the importance of physiological descriptors to constrain LUE, which prevails over structural factors controlling f APAR (i.e., green biomass) under given environmental conditions and encourage the use of hyperspectral remote sensing for diagnostic upscaling of GPP. With sPRI or Fy760 as a proxy for LUE, RSM is presented as a valuable means to diagnose N-induced effects on physiology. Our results show the limits of MM in predicting the spatial and temporal variability of GPP when LUE is not controlled by meteorological drivers alone (VPD, temperature, soil moisture). Accordingly, GPP is eventually biased whenever neither climatic nor structural state variables explicitly reveal spatial changes in the LUE parameter associated with plant nutrient availability; residuals showed a clear tendency to underestimate the highest modeled GPP values, significantly correlated to Fy760 and sPRI (Fig. 7). From a practical point of view, the forcing variables of RSM approaches may show a better observational coverage. In effect, the satellite-based retrievals of RSM forcing variables could additionally overcome representativeness limitations and potential regional or seasonal biases in meteorological fields (Dee et al., 2011). The uncertainties in forcing variables of MM (i.e., temperature, VPD and soil moisture) could propagate and affect the GPP estimates.

Concluding remarks
1. Fy760 and sPRI correlated well with GPP: both increased with N content and decreased with senescence.
2. MTCI can be used as a good descriptor of N content in plants, but the relationship with GPP breaks down under drought conditions.
3. Meteo-driven models were able to describe temporal variations in GPP, and soil moisture can be a key parameter to better track the seasonal dynamics of LUE in arid environments. However, meteo-driven models were unable to describe N-induced effects on GPP. Important implications can be derived from these results, and uncertainties in the prediction of global GPP still remain when meteo-driven models do not account for plant nutrient availability.
4. sPRI or Fy760 provides valuable means to diagnose nutrient-induced effects on the photosynthetic activity and, therefore, should be included in diagnostic GPP models.