Sun-induced Chlorophyll fluorescence and PRI improve remote sensing GPP 1 estimates under varying nutrient availability in a typical Mediterranean savanna 2 ecosystem 3 4

Oscar Perez-Priego*, Jinhong Guan, Micol Rossini, Francesco Fava, Thomas Wutzler, Gerardo 5 Moreno, Nuno Carvalhais, Arnaud Carrara, Olaf Kolle, Tommaso Julitta, Marion Schrumpf, 6 Markus Reichstein and Mirco Migliavacca 7 1 Max Planck Institute for Biogeochemistry, Jena, Germany 8 2 Chinese Academy of Sciences and Ministry of Water Resources, State Key Laboratory of 9 Soil Erosion and Dryland Farming on Loess Plateau, Institute of Soil and Water 10 Conservation, Yangling, Shaanxi, China 11 3 Università degli Studi Milano-Bicocca, Remote Sensing of Environmental Dynamics 12 Laboratory, DISAT, Milan, Italy 13 4 Universidad de Extremadura, Forest Research Group, Plasencia, 10600, Spain 14 Departamento de Ciencias e Engenharia do Ambiente, DCEA, Faculdade de Ciencias e 15 Tecnologia, FCT, Universidade Nova de Lisboa, 2829-516 Caparica, Portugal 16 6 Fundación Centro de Estudios Ambientales del Mediterráneo (CEAM), Valencia, Spain 17 18


Introduction 98
Human-induced nutrient imbalances are affecting essential processes that lead to 99 important changes in ecosystem structure and functioning (Peñuelas et al., 2013). In spite of 100 the crucial role of nutrients in regulating plant processes, efforts to describe and predict the 101 response of photosynthesis to such changes with remote sensing information have been 102 limited. In the framework of the classical Monteith Light Use Efficiency (LUE) model 103 (Monteith, 1972), estimates of photosynthesis (hereafter gross primary productivity, GPP) are 104 based on three key quantities: i) the fraction of photosynthetically active radiation (fAPAR) 105 absorbed by the vegetation, ii) potential LUE (or maximum, LUE m ), normally taken from 106 look-up tables and associated with plant functional types (Heinsch et al., 2006) and iii) 107 correction factors related to meteorological conditions that limit LUE m . Although Nitrogen 108 (N) deficiencies have been recognized one of the main correction factors of LUE m (Madani et 109 al., 2014), the predictive capability of LUE models is usually circumspect as they operate 110 based on the general assumption that plants are under non-limiting nutrient conditions. 111 Very little attention has been given to nutrient-induced effects on fAPAR and LUE in 112 common formulations of LUE models. Light absorption by plant is given by chlorophyll 113 pigments that enable photosynthetic processes. Assuming a correlation between leaf 114 chlorophyll pigments and leaf N content, note that N atoms are basic components of the 115 chlorophylls molecular structure, several studies have demonstrated that leaf N content can be 116 estimated through chlorophyll-related hyperspectral vegetation indices (Baret et al., 2007;117 Schlemmer et al., 2013). Among these indices, the MERIS Terrestrial Chlorophyll Index 118 (MTCI, Dash and Curran, 2004) has been used as a proxy for fAPAR (Rossini et  during photosynthesis can be characterized by two main photo-protective mechanisms: i) non-133 photochemical quenching that can be detected using the Photochemical Reflectance Index 134 (PRI), originally proposed by (Gamon et al., 1992) to track changes in the de-epoxidation 135 state of the xanthophyll cycle pigments, and ii) Chlorophyll fluorescence, the dissipation of 136 energy that exceeds photosynthetic demand (Krause and Weis, 1984). The PRI has been 137  April and achieving senescence by the end of May (Table 1). 191 The experiment consisted of four randomized blocks of about 20 m x 20 m. Each block 192 was separated into four plots of 9 m x 9 m with a buffer of 2 m in between to avoid boundary 193 effects. In each block, four treatments were applied (see Fig. 1 Ancillary measurements were taken in every field campaign as follows: green plant area index 217 (PAI g ) and aboveground biomass were directly measured by harvest in four parcels (0.25m x 218 0.25m) within each plot in the area surrounding that where spectral and flux measurements 219 were taken. All samples were refrigerated just after collection, and transported for laboratory 220 analyses. Fresh samples were separated into functional groups, the sample was scanned and 221 green plant area was measured using image analysis (WinRHIZO, Regent Instruments Inc., 222 Canada). Afterwards, fresh samples were dried in an oven at 65 °C for 48 hours and weighed 223 to determine dry biomass. To analyze the nutrient content in leaf mass, biomass subsamples 224 were ground in a ball mill (RETSCH MM200, Retsch, Haan, Germany) and total C and N 225 concentrations were determined with an elemental analyzer (Vario EL, Elementar, Hanau, 226 Germany). P concentrations were also measured: 100-mg biomass subsamples were diluted in 227 3 ml of HNO 3 65%, (Merck, Darmstadt, Germany) and microwave digested at high pressure 228 For every acquisition, 3 and 10 scans (for Spectrometers 1 and 2, respectively) were averaged 306 and stored as a single file. Five measurements were collected for each plot. Spectral data were 307 acquired with dedicated software (Meroni and Colombo, 2009) and processed with a 308 specifically developed IDL (ITTVIS IDL 7.1.1) application. This application allowed the 309 basic processing steps of raw data necessary for the computation of the hemispherical conical 310 reflectance factor described by Meroni et al. (2011). 311 The following indices were selected as suitable to investigate long term nutrient-mediated 312 effects on photosynthesis. The NDVI (Rouse et al., 1974) Table 2. 327

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

Monteith's light-use efficiency modelling approaches 346
Following Monteith's LUE framework (Eq. 2) two alternative modeling approaches were 347 used: 348 i. Meteo-driven methods (MM); based on the MOD17 formulation, ƒAPAR is 350 approached through the relationship with NDVI and includes limiting functions 351 ƒ(meteo), which are based on climatic driving parameters to limit maximum LUE 352 (LUE max ). Alternatively, Eq.
(2) was reformulated as follows: 353  and 23.5±4.1 mg g -1 , respectively. Although slightly lower, the differences in total N content 428 between C and +P, and +NP and +N remained high over the drying period. Total P content 429 was higher in +NP and +P treatments after fertilization, as compared to +N and C treatments. 430 Consequently, the N:P ratio at the first campaign after fertilization (#2) achieved values of up 431 to 14.2, 6.6, 6, and 3.7, in +N, C, +NP, and +P treatments, respectively. Similar differences in 432 N:P between treatments were also observed during the drying period (#3 and #4, Table 2). On 433 the other hand, PAI g ranged from 0.4 m 2 m -2 in campaign #4 to up to 2.5 m 2 m -2 in campaign 434 #2. No differences were found in PAI g among treatments since grazing apparently offset any 435 potential difference in the green aboveground production. Regarding variations in the fraction 436 of plant forms, no significant differences were found between treatments. 437 Fertilization caused significant differences in the GPP daily (p<0.05) between N-addition 438 treatments (mean values of 19.62±4.15 and 18.19±5.67 μmolCO 2 m -2 s -1 for +N and +NP, 439 respectively) and C and +P treatments (14.31±5.39 and 14.40±4.09 μmolCO 2 m -2 s -1 , 440 respectively) in the peak of the growing season (campaign #2); a relative difference of 37% in 441 GPP daily values was found between +N and +NP and C treatments. During the drying period, 442 however, GPP was substantially down regulated (campaigns #3 and #4) and no significant 443 differences were found in GPP daily , regardless of differences in plant N content observed 444 among treatments. The potential photosynthetic capacity GPP 2000 (Fig 2)  as compared with C and +P treatments (Table 1). During the drying season, GPP 2000 decreased 455 in all treatments ranging between 5.6 and 8 μmolCO 2 m -2 s -1 and no differences among 456 treatments was observed (Fig 2 c and d). These results indicate that the senescence of the 457 herbaceous stratum, which is regulated by water availability, strongly modulated the 458 photosynthetic capacity of the vegetation over the season. 459  (Figure 3). As for GPP, differences between treatments were 467 more evident during campaign #2 when C plots showed statistically lower values for all the 468 indices considered, while only MTCI was able to detect significant differences between N 469 fertilized plots (+N and +NP). Furthermore significant differences in Fy760 and MTCI 470 between C and the other three treatments were found (p<0.05) in the drying period (campaign 471 #4,). NDVI varied significantly with changes in PAI g with values of 0.4 in the campaign #4 472 up to 0.8 in the campaign #2 (p<0.001, r 2 =0.79). 473 474

Relationship between remote sensing data and GPP 475
While Ph indices (Fy760 and sPRI) varied linearly with GPP noon in all treatments 476 (p<0.001, r 2 =0.66 for Fy760 and p<0.001, r 2 =0.79 for sPRI, respectively, Fig 4 a and b,), 477 different patterns were observed for St: NDVI and GPP were best fitted by an exponential 478 regression (p<0.001, r 2 =0.77 Fig 4 c), while a weak linear relationship between MTCI and 479 GPP noon (p<0.05, r 2 =0.45, Fig 4 d) was found. Although a weak relation between MTCI and 480 GPP noon was found, MTCI was strongly correlated with plant N content (y=14.17x-2.49, 481 p<0.001, r 2 =0.86). Note that these results are computed excluding data taken in the pre-482 treatment campaign (#1) and differences in the relationship between remote sensing data and 483 GPP noon among treatments can be only attributed to nutrient-induced effects. The ANCOVA 484 test did not show significant differences neither in slope nor intercept of the relationship 485 between GPP noon and sPRI, and NDVI across treatments. However, barely significant 486 differences were found in the relationship between GPP noon and Fy760 (p<0.1, Fig 4b) and 487 significant between GPP noon and MTCI (p<0.01, Fig 4d) between N addition treatments (+N 488 and +NP) and C treatments (C and +P). 489 Similar to GPP noon , GPP 2000 was also significantly related to mean midday sPRI 490 (r 2 =0.76, p<0.001, Fig. 5a) and Fy760 (r 2 =0.76, p<0.001, Fig. 5b). As expected, an 491 exponential regression fitted best for NDVI, while a poor relationship with MTCI was found 492 (data not shown).  Figures 6a-6c). 499 Note that the highest biases in modeled GPP noon values among MM models belong to +N and 500 +NP treatments in field campaign #2. Since the four treatments experienced the same 501 environmental conditions (i.e. comparable values of SWC, VPD, air temperature), this bias 502 can be attributed to the higher N content (+N and +NP treatments) as compared to C and +P 503 treatments. Remarkably, residuals of the MM (VPD-SWC) taken from periods with moist soil 504 (SWC>15) were significantly correlated with sPRI and Fy760 (p<0.05, Fig. 7 a and b,  505 respectively). However, no biases between residuals and predictions were observed in RSM 506 over the span of values and treatments (Fig. 8). Results from the evaluation of model 507 performance indicated that RSM performs best when NDVI rather than MTCI, is used as St in 508 the Eq.7 and, hence, as a proxy for fAPAR (Table 3). Our results indicated that RSM 509 performs best when either Ph (sPRI or Fy760) is combined with NDVI as St. In this study we also explored the capability of remote sensing to describe ecosystem 536 functional properties defined as those quantities that summarize and integrate ecosystem 537 processes and responses to environmental conditions and can be retrieved from ecosystem 538 level fluxes (e.g. GPP 2000 ) and structural measurements (Reichstein et al., 2014). GPP at light 539 saturation (i.e. GPP 2000 ) is one example of an ecosystem functional property, shown here to be 540 quite correlated to sPRI and Fy760 (Fig. 5). This result suggests that sPRI and Fy760 open 541 also new opportunities for remote sensing products to describe the spatiotemporal variability 542 of essential descriptors of ecosystem functioning (Musavi et al., 2015). Inferring GPP 2000 543 using remote-sensing has important implication both for monitoring global carbon cycle and 544 for benchmarking terrestrial biosphere models. 545 MTCI was tightly related with N content (r 2 =0.86, p<0.001), independent of other structural 546 variables (i.e. PAI g ), and can be used as a good indicator of N availability. Although MTCI 547 has been proven to be very sensitive to variations in chlorophyll contents (Dash and Curran, 548 2004) and hence linkable with light absorption processes, it was weakly correlated with GPP, 549 particularly in plots added with N (+N and +NP; r 2 =0.27, p<0.01 , Fig 4 d). A quite wide range 550 of GPP noon values were found at high values of MTCI -high GPP noon values corresponding to 551 the growing season and low ones to the drying period -which can be explained by two 552 simultaneous mechanisms. 553 First, despite the high plant N content, physiological mechanisms including stomatal control 554 or reduced carboxylation efficiency down-regulate GPP (Huang et al., 2004) and ultimately 555 might break the relationship between GPP noon and MTCI. Second, MTCI tracks changes in N 556 content regardless changes in canopy structure occuring during the dry season when grass 557 achieved senescence (i.e. green to dry biomass ratio, PAI g ). More studies aimed at the 558 separation of the combined effects of N and changes in green/dry biomass fractions on 559 fAPAR are essential. On the other hand, although NDVI followed the seasonal dynamic of 560 PAIg, it saturated at high GPP noon values indicating the low ability of this index to detect 561 spatial variations induced by N fertilization. 562 Although optical measurements were taken at high spatial resolution (<0.36 m 2 ), the 563 separation of confounding factors affecting sPRI or Fy760 is essential to elucidate the 564 mechanistic association between sPRI or Fy760 and GPP. Like sPRI, the retrieval of Fy760 565 from the apparent reflectance signal can be also affected by vegetation structure or canopy 566 background components (Zarco-Tejada et al., 2013). After optimization and selection of the 567 best model parameters using NDVI and sPRI (or Fy760) as driver, we analyzed the response 568 of simulated GPP to variations in NDVI and sPRI (or Fy760 , Fig 9). Results indicate that at 569 high GPP levels, Fy760 and sPRI but less NDVI shaped GPP. However, at low GPP levels, 570 either Fy760 or sPRI responded to GPP on a small scale (Fig 9b). and semi-arid ecosystems (Prince and Goward, 1995). However, residual analyses 593 demonstrated that MM (VPD-SWC) was unable to track N-induced differences in GPP during 594 the growing period, when both parameters are not limiting (Fig. 7). By contrast, accurate 595 estimates of GPP were obtained with RSM both over the drying and the growing periods. 596 These results also indicate the importance of physiological descriptors to constrain LUE, 597 which prevails over structural factors controlling fAPAR (i.e. green biomass) under given 598 environmental conditions and encourage the use of hyperspectral remote sensing for 599 diagnostic upscaling of GPP. 600 With sPRI or Fy760 as a proxy for LUE, RSM is presented as a valuable means to diagnose 601 N-induced effects on physiology. Our results show the limits of MM in predicting the spatial 602 and temporal variability of GPP when LUE is not controlled by meteorological drivers alone 603 (VPD, temperature, soil moisture). Accordingly, GPP is eventually biased whenever neither 604 climatic nor structural state variables explicitly reveal spatial changes in the LUE parameter 605 associated with plant nutrient availability; residuals showed a clear tendency to underestimate 606 the highest modeled GPP values, significantly correlated to Fy760 and sPRI (Fig.7). Fy760 are indistinctly used as a proxy of LUE, the NDVI is taken as fPAR. 696 697 698  Table 1. Ancillary data resulting from the analysis. Green Plant Area Index (PAIg), fraction of PAI in different plant forms (fPAI), 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. physiologically-related parameter of RSM referring to either sPRI or Fy760 as a proxy for LUE; PLRC, photosynthetic light response curve; PRI, photochemical reflectance index; R eco , daytime ecosystem respiration; RSM, remote sensing based models; SIF, sun-induced chlorophyll fluorescence; sPRI, scaled-photochemical reflectance index; st, structurally-related parameter of RSM referring to either NDVI or MTCI as a proxy for fAPAR; SWC, soil water content; SWC max parameter of the f(meteo) term; VPD, air-to-leaf vapor pressure deficit; VPD max and VPD min are fitting parameters of the f(meteo) term; α is a parameter describing the photosynthetic quantum yield; β is the parameter that extrapolates to GPP at saturating light condition.