Combining hyperspectral remote sensing and eddy covariance data streams for estimation of vegetation functional traits

Remote Sensing (RS) has traditionally provided estimates of key biophysical properties controlling light interaction with the canopy (e.g., chlorophyll content (Cab) or leaf area index (LAI)). However, recent and upcoming developments in hyperspectral RS are expected to lead to a new generation of products such as vegetation functional traits that control leaf carbon and water gas exchange. This information is pivotal to improve our understanding and capability to 25 predict biosphere-atmosphere fluxes at global scale. Yet, the retrieval of key functional traits such as maximum carboxylation rate (Vcmax) or the Ball-Berry stomatal sensitivity parameter (m) remains challenging, as they only have a weak and indirect influence on optical reflectance factors. Recently, the assimilation of different observations in coupled soilvegetation-atmosphere transfer (SVAT) and radiative transfer models (RTM) is allowing Vcmax and m estimates; notably using the Soil Canopy Observation of Photosynthesis and Energy fluxes (SCOPE) model. In this work we assess the 30 potential of airborne and satellite emulated hyperspectral imagery jointly with eddy covariance (EC) data for the retrieval of functional traits. Specifically, we made use of time series of gross primary production (GPP) and thermal irradiance measured with net radiometers, together with 17 hyperspectral airborne images. The potential of satellite-borne sensors was tested with emulated EnMAP imagery from the airborne data. EnMAP was selected because of the availability of the emulator, and because is one of the foreseen hyperspectral satellite missions expected to contribute to a new generation of 35 RS products. We estimated ecosystem functional traits by inverting the senSCOPE model, a novel version of SCOPE https://doi.org/10.5194/bg-2019-501 Preprint. Discussion started: 10 February 2020 c © Author(s) 2020. CC BY 4.0 License.

RS of vegetation function faces the fact that optical and TIR signals that are mechanistically related with photosynthesis and transpiration originate inside the leaf. Therefore they vary within the canopy (light or wind speed gradients, among others) 85 and are modified in their transfer from leaves to top of the canopy (TOC), according to sun-view geometry and vegetation structure. This complicates the application of empirical methods to understand or generalize the connection between RS observations and plant function. For this reason, there is an increasing on the use of models physically describing both the radiative transfer and physiological process. Due to the complexity of these processes and the weak influence that some functional traits have on RS observations, some authors argued about the need to jointly use RS and EC data for their 90 retrieval. For instance, Pacheco-Labrador, et al., (2019) showed the successful retrieval of V cmax , m, and Cab, jointly constraining the Soil Canopy Observation of Photosynthesis and Energy fluxes (SCOPE) (van der Tol et al., 2009) model using with chamber-based fluxes and proximal sensing data. In fact, in the last years the interest of two different communities seems to converge. TBM have used RS data to constrain traits or to inform biophysical variables into models (Alton, 2017;Xie et al., 2018) whereas RTM have been also coupled to simplified vegetation dynamic (Koetz et al., 2005) or 95 physiological models (Xin et al., 2015), as well as to TBM (Migliavacca et al., 2009) or crop models (Thorp et al., 2012;Dorigo et al., 2007) in order to assimilate RS observations. However, to our knowledge an attempt to jointly retrieve functional traits using hyperspectral imagery combined with EC data is lacking in the literature.
Nowadays, some RTM simulate optical (e.g., PRI, SIF) and / or TIR signals related with plant photosynthesis taking place at leaf level (van der Tol et al., 2009;Vilfan et al., 2018;Yang et al., 2017;Hernández-Clemente et al., 2017). However, this 100 modeling is not informative of the physiological processes originating such signals. This connection can only be achieved coupling RTM to models representing such processes as for example the soil-vegetation-atmosphere transfer models https://doi.org/10.5194/bg-2019-501 Preprint. Discussion started: 10 February 2020 c Author(s) 2020. CC BY 4.0 License.
(SVAT). The state-of-the-art model coupling optical, SIF and TIR RTM with energy balance and photosynthesis models is SCOPE (van der Tol et al., 2009). Further improvements of this model are mSCOPE , which allows representing vertically heterogeneous canopies, and senSCOPE (Pacheco-Labrador et al., 2020), which improves the 105 representation of canopies featuring mixed green and senescent leaves (Pacheco-Labrador et al., 2020). Using satellite imagery, SCOPE has been used to obtain estimates of physiological parameters of vegetation such as V cmax and / or m exploiting reflectance factors (R λ , where λ denotes spectral) and SIF (Zhang et al., 2014), R λ , SIF and EC fluxes (Zhang et al., 2018), R λ and TIR data (Bayat et al., 2018), or R λ and EC fluxes (Dutta et al., 2019). Also, proximal sensing data have been used to constrain SCOPE providing estimates of functional parameters (Pacheco-Labrador et al., 2019;Hu et al., 2018); and 110 more recently, airborne imagery has been used to retrieve V cmax from R λ and SIF (Camino et al., 2019).
From all these works, only Camino et al, (2019) validated their retrievals against actual measurements of functional traits (V cmax ), as gas exchange leaf-level measurements usually needed to obtain validation data of functional parameters are time consuming and are not feasible for large areas and / or diverse natural ecosystems. In these circumstances, pattern-oriented model evaluation has been used to assess the suitability of different models and parameter estimates (Carvalhais et al., 115 2014;Reichstein et al., 2011;Migliavacca et al., 2013;Grimm and Railsback, 2012;Luo et al., 2012). Pacheco-Labrador et al, (2019) used pattern-oriented model evaluation approach to assess the suitability of functional and biophysical estimates obtained by SCOPE and senSCOPE model inversion against ground observations. On one hand, the relationships between variables measured at canopy scale (e.g., nitrogen concentration) with parameter estimates (e.g., V cmax and chlorophyll content (C ab )) were compared with relationships expected from the theory or the literature (e.g. Feng and Dietze (2013)). On 120 the other hand, variables derived from predicted fluxes (e.g., evaporative fraction (EF)) were compared with observations from a nearby EC system. Such evaluations provided relevant information about models structure and ill-posed solutions of their inversion.
The new generation of hyperspectral satellite-borne sensors (Rast and Painter, 2019) brings new possibilities for the spatiotemporal characterization functional traits of vegetation, even though the algorithms to convert these data into 125 information must be developed and evaluated. Some of these sensors are already operational (DESIS (Kerr et al., 2016), PRISMA (Galeazzi et al., 2008)); or will be in the next few years (e.g., EnMAP (Guanter et al., 2015), HyspIRI (Lee et al., 2015), etc. see Rast and Painter (2019) for a complete overview). Meanwhile up-coming missions also offer emulation capabilities (e.g., the EnMAP End-to-End simulator (Segl et al., 2012)); which allows the evaluation of their potential for different applications before they operate. 130 In this work we aim at developing a methodology to use of hyperspectral remote sensing in combination with EC stations to estimate functional and biophysical parameters of vegetation at the ecosystem scale. We constrain senSCOPE model with a combination of multiple hyperspectral R λ in different phenological periods, diel tower-based TIR observations and GPP estimates in a Mediterranean tree-grass ecosystem manipulated with nitrogen (N) and N plus phosphorus (P). The inversion is tested both on airborne as well as on synthetic EnMAP hyperspectral imagery, simulated to increase uncertainties and down-grade resolutions to levels expected in a space-borne sensor. Biophysical and functional parameters are assessed using ground observations and pattern-oriented model evaluation approach.

Study site
The study area is located in Majadas de Tiétar (39º56'24.68'' N,5º45'50.27'' W;Cáceres,Spain). This is a Mediterranean 140 tree-grass ecosystem where a large scale manipulation experiment with N and P was conducted at the ecosystem scale (MANIP, see El-Madany et al., (2018) and Nair et al., (2019)). Three EC towers monitor three areas under different manipulation regimes. The first one (CT) operates since 2003 and is used as control, the second and the third ones operate since 2014 and their footprints -~ 24 ha-, have been fertilized with N (NT) and N+P (NPT) respectively. Towers are sufficiently separated (~500 m) to prevent 80 % iso-lines of their respective footprint climatology to overlap (El-Madany et 145 al., 2018). These areas were fertilized twice: 50 and 10 kg P ha -1 in the form of triple superphosphate (Ca(H 2 PO 4 ) 2 ) were applied in November 2014 and March 2016, respectively in NPT. 100 and 20 kg N ha -1 were applied in March 2015 and March 2016 in the form of ammonium nitrate (NH 4 NO 3 ) in NPT, while in NT was in form of calcium ammonium nitrate (5Ca(NO 3 ) 2 NH 4 NO 3 ) to compensate the Ca added in NPT with the P fertilizer. Fertilization modified species, ecosystemlevel biophysical and functional properties and fluxes (Nair et al., 2019). 150 The climate is Continental Mediterranean, mean annual temperature and precipitations are 16.7 ºC and ~650 mm, respectively. There is a strong seasonality with hot and dry summers reaching ~40 ºC of maximum daily temperature, and with scarce precipitation -less than the 6 % of the annual rainfall (Casals et al., 2009). The ecosystem combines two vegetation layers: grass and trees. The grassland layer is composed of a large number of annual species comprehending three main PFT (grasses, forbs and legumes, see Migliavacca et al, (2017)), and supports low grazing pressure (<0.3 cows / ha). It 155 features a strong phenological behavior controlled by radiation and water availability in spring-summer, which produces a mayor and a minor greening peaks in spring and autumn, respectively; and a dry season in summer (Luo et al., 2018).
Replacement of species and mortality accumulates senescent and death standing material through the year, which can already represent up to the 30 % of the leaf area before the beginning of the dry season. The ecosystem is composed by scattered trees, mainly Quercus ilex L. subsp. ballota [Desf.] Samp (Evergreen Holm Oak). Trees fractional cover is ~20 %, average 160 tree distance, tree height, horizontal and vertical crown radius are 18.8 m (standard deviation, σ = 5.0 m) and 7.9 m (σ = 0.9 m), 4.2 m (σ = 0.9 m) and 2.7 m (σ = 0.9 m), respectively . Evergreen Holm oak leaves partially renovate every spring so leaves of up to three cohorts can coexist within the same crown. Some years, new leaves can sprout twice, though biophysical and spectral properties leaves from different cohorts are quite similar after the first 9 months of development (González-Cascón et al., 2016;Gonzalez-Cascon et al., 2018;Pacheco-Labrador et al., 2017b). 165 https://doi.org/10.5194/bg-2019-501 Preprint. Discussion started: 10 February 2020 c Author(s) 2020. CC BY 4.0 License.

Hyperspectral Airborne and EnMAP synthetic imagery
Airborne images were acquired in seven different campaigns that took place in the study site between May 2011 and May 2017 covering the spring growing and the summer dry seasons (Table 1). A total of 17 sites-images were used for the analysis. Imagery was acquired by the Compact Airborne Spectrographic Imager CASI-1500i (Itres Research Ltd., Calgary, AB, Canada), operated by the Instituto Nacional de Técnica Aeroespacial (INTA). CASI featured 144 spectral bands, 170 spectral range ~350-1050 nm, full width half maximum ~5.5 nm, and field of view ~40 º. Images were acquired around solar noon, at least always once in the main wind direction (east-northeast and west-southwest) along the longest axis of the EC tower footprints, with a pixel size ~0.90 x 1.66 m approximately. However, for all the campaigns after 2012 also overpasses with different azimuth orientations also were acquired. Images were atmospherically corrected by INTA using ATCOR-4 TM (ReSe Applications GmbH, Germany) and improved with empirical line correction (Smith and Milton, 1999) using ground 175 calibration targets measured during flight campaigns with an ASD Fieldspec3 TM spectroradiometer (Analytical Spectral Devices, Inc., Boulder, CO, USA). For most of the campaigns, ancillary measurements of water vapor and aerosol optic thickness were acquired using a CIMEL CE318-NE sun photometer (Cimel Electronique, Paris, France). Nearest neighbor resampling was used during geometric correction to prevent spectral mixture of different covers; output spatial resolution was 1 m. For the campaigns after 2012, several overpasses acquired in the main wind direction as well as in the solar 180 principal plane were combined in mosaics featuring pixels of 1.5 m, which were used instead of the individual overpasses to reduce directional effects of the sun-view geometry. The combination of single overpasses and mosaicked is not ideal, since the assimilated data might be subject to different directional effects (see Discussion section). Dates of each overpass can be seen in Table 1. In order to test the potential of satellite-borne hyperspectral imagery featuring higher uncertainties and coarser spectral and spatial resolutions than airborne data, we used the airborne images as a reference and simulated synthetic EnMAP imagery with the EnMAP End-to-End simulator (Segl et al., 2012), available for this work. Since no other satellite-borne hyperspectral data was available at the site, the simulation was a reasonable alternative. Synthetic imagery featured EnMAP 190 https://doi.org/10.5194/bg-2019-501 Preprint. Discussion started: 10 February 2020 c Author(s) 2020. CC BY 4.0 License. spectral configuration (FWHM ϵ [5.75,9.81] nm between 423 and 930 nm), spatial resolution (~30 m) and include geometric and atmospheric uncertainties (Guanter et al., 2015). However, in fact these uncertainties were added to those already present in the CASI datasets; which made these images a rougher test for the methodology under evaluation. Squared cut-outs of 1 km long centered on each EC tower were used for the later analyses both for CASI and EnMAP synthetic images.

Eddy Covariance and Footprint Climatology 195
Three EC towers monitored ecosystem carbon dioxide (CO 2 ), and water vapor concentrations with an enclosed-path infrared and Zonen, Delft, Netherlands) also at a height of 15 m. Quality assessment and quality checks were conducted according to Mauder and Foken (2005). Low quality data, and data under low turbulent mixing (Papale et al., 2006) or during rainy events were discarded. GPP was computed subtracting ecosystem respiration estimated according to Reichstein et al. (2005) from net ecosystem exchange using the REddyProc R Package (Wutzler et al., 2018). From the final dataset only daytime highquality GPP data (i.e. coming from half hours where NEE was measured with good quality) were used. Fluxes uncertainty 205 (σ) was computed from the standard deviation of the marginal distribution sampling of the gap-filling procedure (Reichstein et al., 2005). Footprint climatology was calculated every half-hour according to Kljun et al. (2015) and used later minimize the spatial mismatch between EC footprint and RS imagery. Further details can be found in Perez-Priego et al., (2017) and El-Madany et al., (2018).

Biophysical variables sampling and up-scaling 210
Biophysical variables were estimated at ecosystem scale for evaluation of parameter estimated via inversion of senSCOPE.
Destructive sampling of grass and tree leaves was always performed simultaneously to each airborne campaign, as well as (C w , g cm -2 ) and dry matter content (C dm , g cm -2 ) were determined in the laboratory with gravimetric methods and a scanner following protocols described in Mendiguren et al, (2015) and Melendo-Vega et al, (2018). For the grass layer, after 2012, these variables were separately measured for green and senescent plants; which allowed determining the fraction of green leaf area (f green ). Trees f green was assumed 1. Tree leaves from the current year and previous years were sampled from the north and south sides of between 5-18 crowns in each campaign; C w and C dm were determined using the same methods than 220 used for the grasses. Tree LAI was determined combining different methods: indirect measurements using the LICOR LAI 2200-C instrument (2018) and direct estimations using litter traps and leaf turnover rates (Melendo-Vega et al., 2018). Due https://doi.org/10.5194/bg-2019-501 Preprint. Discussion started: 10 February 2020 c Author(s) 2020. CC BY 4.0 License.
to lower temporal frequency of tree LAI measurements, a seasonal model was used to predict tree LAI as a function of DoY combining these data. The reduced tree fraction cover and tree LAI variability compared with grasses minimize the uncertainties introduced by this approach when integrated at ecosystem scale. Nitrogen concentration per total (N mass , %) and 225 green unit mass (N mass,green , %) was also determined. When missing, seasonal models were used to predict N mass and N mass,green as well as f green (the last two in grasses). Trees C ab and carotenoids content (C ca ) were determined from SPAD measurements and a model calibrated from samples of the site (Gonzalez-Cascon et al., 2017). From 2016 on, grass C ab and C ar were determined via destructive sampling (Gonzalez-Cascon and Martin, 2018). The relationship found between C ab per unit total leaf area and N mass (Pacheco-Labrador et al., 2020) was used to estimate grass C ab in the campaigns where it was missing. 230 N mass,green was also used to estimate grass V cmax using the relationship in Feng and Dietze (2013); since no other data were available, a constant V cmax = 45 μmol m -2 s -1 was estimated as an acceptable value for trees according to literature (Vaz et al., 2010;Vaz et al., 2011;Limousin et al., 2010).
Since 2014, the herbaceous layer was sampled in between 6 and 33 quadrants of 30 x 30 cm inside each tower footprint for determination of isotope signature of plant material (δ 13 C, ‰) during the airborne campaigns. Tree leaves were sampled 235 inside each tower footprint, a total of 8 trees for each footprint and 4 branches sampled along 4 different cardinal positions.
This sampling was done in winter, when the differences in terms of phenology between leave cohorts are minimal across trees and the traits more representative of the average conditions over the season. The carbon isotopic composition of dried samples was analyzed using a Delta Plus isotope ratio mass spectrometer (Thermo Fisher, Bremen, Germany) coupled via a ConFlowIII open-split to an elemental analyzer (Carlo Erba 1100 CE analyzer; Thermo Fisher Scientific, Rodano, Italy), and 240 measured according to Werner et al. (1999;2001) and Brooks et al. (2003). δ 13 C was calculated using Eq. 1 (Brand, 2013;Coplen, 2011), scaled to the δ 13 C IAEA-603 -LSVEC scale, based on calibrated in-house-standard (Acetanilide: -30.06 ± 0.05 where 13 R sample and 13 R standard are 13 C/ 12 C ratio of the sample and of the standard, respectively. 245 All the variables defined per leaf area were up-scaled to ecosystem level according to grass and tree LAI, fraction vegetation cover tress and grasses (f VC,grass = 100 %, f VC,trees = 20 %, respectively), and seasonal estimates of fractions of current and former years leaves for the trees (see Melendo-Vega et al, (2018)). Grass V cmax was scaled using LAI in the green fraction.
Variables relative to mass, such as N mass , N mass,green and δ 13 C were scaled considering also C dm . Scaled δ 13 C was used to compute 13 C photosynthetic discrimination (Δ 13 C) as a function of atmospheric CO 2 and plant-measured isotope signature 250 (δ 13 C atm and δ 13 C grass , respectively) following Eq. 2: https://doi.org/10.5194/bg-2019-501 Preprint. Discussion started: 10 February 2020 c Author(s) 2020. CC BY 4.0 License.

senSCOPE model parametrization
senSCOPE (Pacheco-Labrador et al., 2020) is an extended version of the model SCOPE that separates radiative transfer, energy balance and photosynthesis for green and senescent leaves randomly mixed within a homogeneous canopy. While 255 SCOPE parametrizes a single leaf type, senSCOPE requires leaf parameters for both senescent and green leaves; where senescent leaves only present senescent pigments and green ones feature any pigment but senescent. Like in Pacheco-Labrador et al., (2020), during inversion we assumed the same C dm for both leaf types, but that C w of green leaves is four times higher than in senescent leaves (Kidnie et al., 2015). These assumptions allowed us to later invert the model optimizing averaged leaf parameters rather than the individual parameters of each leaf type. senSCOPE includes f green , 260 representing the fraction of leaf area corresponding to green leaves; as in Pacheco-Labrador (2020) we used a neural network model to predict f green as a function of the averaged leaf parameters (X). Therefore, during inversion f green is predicted as a function of the averaged leaf parameters, and then the parameters corresponding to green and senescent leaves (X green and X senes , respectively) are calculated solving Eq. 3: where r is the ratio X senes /X green , which is different for each parameter according to the assumptions described above; and f senes = 1-f green is the fraction of leaf area corresponding to senescent leaves. The neural network predicting f green was trained over a look-up table (LUT) of 5000 samples generated via Latin Hypercube Sampling of different leaf parameters and f green . C ca was forced as a function of C ab according to Sims and Gamon (2002), with random normal noise of μ = 0 and σ = 4.5 μg cm -2 .
Parameters were limited according to the same bounds set for the inversion (see section 2.6); except C ab , C ca and senescent 270 pigments content (C s ) whose maximum values were 50 μg cm -2 , 20 μg cm -2 and 4 arbitrary units (a.u.), respectively. Different test showed that the emulator must be trained pigment contents bounded with values close to those of the vegetation represented when f green is close to 1 or 0. Averaged parameter values were simulated mixing the LUT pigments according the simulated f green , assuming that in green leaves C s = 0 a.u., and that in senescent leaves C ab = C ca = 0 μg cm -2 .
No assumptions were done about C w and C dm and the LUT values were used as the averaged value from both leaf types. The 275 neural network was then trained to predict f green from the average values using the 60 % of the LUT for training and the remaining 40 % for validation. Training R 2 , root mean squared error (RMSE) and mean error (ME) were 0.817, 0.124 and -.001, respectively. Analogously, validation statistics were 0.778, 1.75 and 0.106. senSCOPE uses the brightness-shape-moisture (BSM) to represent soil optical properties. In this work we used the parameters determining soil bright estimated in Pacheco-Labrador et al, (2019); and we used superficial (5 cm) soil moisture 280 content (SM p , vol vol -1 ) registered in the EC towers to force the effect of moisture on soil reflectance factors.
Other model parameters were also determined as a function of forcing variables. Roughness length for momentum of the canopy (z 0 , m) and canopy displacement height (d, m) were calculated as fractions of canopy height (h c , m) using the factors 0.15 (Plate and Quraishi, 1965) and 0.66 (Brutsaert, 1982), respectively. Also, leaf drag coefficient (C d,l ) was calculated according to Campbell (1977) as a function of wind speed (u, m) and d (Eq. 4) Soil boundary layer resistance (r bs , s m -1 ) was determined as a function of the friction velocity (u * , m s -1 ) internally calculated by the senSCOPE, according with Thom (1972) Soil resistance to evaporation from the pore space (r ss ) was prescribed as a function of SM p using the model fitted in 290 Pacheco-Labrador et al, (2019) in the same study site from lysimeters data (Perez-Priego et al., 2017); however, during the second step of the inversion (see Sect. 2.6), this variable was also estimated.  around the overpass (GPP ovp ). In Step#1, the biophysical parameters of the model (Table 2) and a first guess of V cmax (V cmax,S1 ) were estimated minimizing Eq. 6. 305

Model inversion and evaluation
where subscripts "obs" and "pred" stand for observed and predicted variables, λ is the wavelength in nm and σ GPP is the uncertainty of GPP. All the forcing variables of senSCOPE (meteorological conditions, SM p , etc.) were linearly interpolated to the time of the overpass. Then, parameter posterior uncertainties were calculated as proposed in Omlin and Reichter (1999); where the Jacobian matrix (J) provided correlated posterior distributions of the parameters that were truncated using the truncated Normal and Student's t-distribution toolbox (Botev, 2017;Botev and Ecuyer, 2015). The standard deviation of 200 realizations of V cmax,S1 was used as the posterior uncertainty of this variable ( cmax,S1 ) with a minimum threshold of 5 μmol m -2 s -1 . 315 In Step #2, diel cycles of GPP and up-welling longwave irradiance ( T ↑ ) ±1 day around the overpass were used to estimate the functional parameters V cmax and m; only daytime (down-welling short wave irradiance (R in ) > 60 W m 2 and sun zenith angle (θ s ) < 60º) data were used. Since Step#2 relies on a large number of observations, half-hourly GPP and T ↑ were not further smoothed. Also, the first guess of V cmax,S1 and its posterior uncertainty (with a minimum value) were used to regularize the solution. Unlike in Pacheco-Labrador et al, (2019), r ss was not prescribed in Step#2, but also estimated in this 320 step assuming a constant value for the three days where predictions and observations were compared. Attempts to prescribe https://doi.org/10.5194/bg-2019-501 Preprint. Discussion started: 10 February 2020 c Author(s) 2020. CC BY 4.0 License. this parameter suggested that the estimates might not be representative of the heterogeneity of soil conditions in the whole area. Step#2 cost function is presented in Eq. 7: where t represents each timestamp of the time series of GPP and T ↑ , T,obs ↑ is the uncertainty of T ↑ and and γ is a 325 regularization factor equal 100 and V cmax,S2 is the estimated V cmax in Step #2.
After Step#2, the J is numerically computed for all the estimated parameters in Step#1 and #2, and used to estimate their corresponding posterior uncertainties as formerly described. Uncertainties were propagated to the fluxes and R λ predicted by senSCOPE running 200 realizations of the model.    Fig. 3a,b present the observed and predicted NDVI for CASI and EnMAP, respectively. Observed NDVI is computed from the R λ in each tower used to constrain the inversion. Analogously, Fig. 3c,d present the ME of the fit of R λ (ME R ) in the visible and the NIR regions for both sensors. Observed NDVI shows the 345 phenological state of vegetation in each campaign, and reveals differences induced by fertilization since Apr-2015, being the NDVI of fertilized towers higher than CT. Campaigns in Oct-2012 and Jul-2015 coincided with the dry season, and NDVI values are low, whereas the campaign in May-2017 took place during the transition from the green peak to senescence.

NDVI differences induced by fertilization disappear during the dry season.
R λ shows a good fit after the inversion. ME R were similar for both sensors, but in general R λ was more accurately fit for 350 EnMAP. The largest errors were found during the dry season, when R λ was underestimated in the visible region and overestimated in the NIR. CASI RMSE in the visible and NIR regions were 0.0039 and 0.0093, respectively, whereas EnMAP RMSE in the same regions were 0.0028 and 0.0080, respectively.   (Fig. 4a, b, c) and EnMAP imagery (Fig. 4d, e, f) using Total Least Squares to account for similar error on the y and x axis (Golub and Loan, 1980). Green and orange colors indicate spring (NDVI ≥ 0.40) and dry season 360 (NDVI < 0.40), respectively. LAI (Fig. 4a,d) was slightly overestimated, especially for low values; however mid values were underestimated with CASI (Fig. 4a). f green (Fig. 4b,e) was underestimated for both sensors in the dry season, but more closely predicted in spring. C dm (Fig. 4c,f) Fig. 5a,b evaluate predicted and observed chlorophyll content per unit ground area (C ab,ground ) against N content in green vegetation per unit ground area (N ground,green ). Estimated parameters pretty well reproduced the C ab -N relationship followed by field observations and literature. Fig. 5d,c shows V cmax (per unit total leaf area, V cmax,total ) against N mass ; in all the cases, V cmax,total estimates positively correlated with N mass , but values are lower than the ones estimated from field data and literature. The relationships are not significant (CASI) or weak (EnMAP); this is in part due to the low range of variation of 375 N mass (and therefore of V cmax,total ) (e.g., see Pacheco-Labrador et (2019)), so that the uncertainty of the retrieval prevented significant relationships. Field V cmax,total was estimated using N mass,green and relationships from the literature. As expected, both estimates of C ab,ground and V cmax,total presented the lowest values during the dry season. Similarly, spring m (Fig. 5e,f) estimates and predicted underlying water use efficiency (uWUE) (Fig. 5g,h) were evaluated against measurements of Δ 13 C. Significant positive relationships were found between m and Δ 13 C for CASI, whereas significant and negative relationships were found 380 between uWUE simulated with senSCOPE after the inversion and Δ 13 C for both sensors.
https://doi.org/10.5194/bg-2019-501 Preprint. Discussion started: 10 February 2020 c Author(s) 2020. CC BY 4.0 License. Figure 5. Pattern-oriented model evaluation of: chlorophyll content per unit ground area vs. nitrogen content per unit ground area (a,b), maximum carboxylation rate per unit total leaf area vs. nitrogen concentration (c,d), Ball-Berry stomatal sensitivity vs. 13 C isotope discrimination (e,f), and predicted underlying water use efficiency vs. 13 C isotope discrimination (g,h). CASI estimates/predictions are 385 presented always on the left (a,c,e,g) and EnMAP estimates on the right (b,d,f,h) of each pair of subplots. * stands for p < 0.05 and • for 0.05 ≤ p < 0.10 Fig. 6 compares observed and predicted fluxes after inversion for CASI (Fig 6 a-e) and EnMAP (Fig 6 f-j). GPP (Fig. 6a,f) and T ↑ (Fig. 6b,g) are constraints of the inversion; the first was accurately fit, whereas the second was overestimated (shifted ~30 W m -2 ). λE (Fig. 6c,h) was overestimated in the dry season, whereas H (Fig. 6d,i) was biased (slopes ~0.45). 390 Consequently, EF resulted overestimated in the dry period, and underestimated for low values in spring. However, it was well predicted for values > 0.8, when evapotranspiration dominates. Differences between predicted and observed λE and H fluxes were strongly related with the energy balance closure error (ε EB ) of the EC data (Fig. A1); showing R 2 ~ 0.84 and RMSE ~ 32.5 W m -2 . Median ε EB was ~111 W m -2 (σ ~ 81 W m -2 ), which in relative terms represents the ~21.5 % (σ ~ 14.2 %) of the observed net radiation. However, during the dry season, predicted λE and H error are slightly biased respect to ε EB . 395 https://doi.org/10.5194/bg-2019-501 Preprint. Discussion started: 10 February 2020 c Author(s) 2020. CC BY 4.0 License.  Fig. 7 presents the estimates of r ss compared vs SM p . and the exponential decay model was fit. CASI estimates (Fig. 7a) were very low in spring for SM p > 12 %, and then increase below this threshold. Results for EnMAP are similar (Fig. 7b). In both sensors, some of the estimates in SM p ~ 15 % present high values, separating from the expected decreasing relationship.

Discussion
The new generation of hyperspectral imaging satellite missions is expected to improve the spatiotemporal characterization of vegetation function and functional traits Butler et al., 2017). To achieve this target, advances in uncertainty quantification of multiple data streams, algorithm development, as well as modelling and data integration are 410 essential . In this context, we prove the capability of hyperspectral imagery jointly with data from EC stations and coupled RTM-SVAT models to inform about vegetation functional traits. We show that an inversion of senSCOPE model in a Mediterranean tree-grass ecosystem using multiple constraint (GPP, T ↑ and remote hyperspectral reflectance factors) lead to robust estimates of key biophysical and functional traits of vegetation. The characterization of vegetation functional traits is necessary to improve our understanding and simulation of terrestrial water and carbon fluxes. 415 In this context, EC networks are a keystone for the development, evaluation and benchmarking of TBM (Bonan et al., 2011;Williams et al., 2009) as well as statistical models . These models are expected to predict later flux information 'everywhere, all of the time' (Baldocchi, 2014); partly thanks to the use of global information about vegetation biophysical parameters provided by RS. One of the limits of this approach is the lack of knowledge on the spatiotemporal distribution of key functional traits controlling these fluxes (Rogers, 2014;Rogers et al., 2016;Schaefer et al., 2012); which 420 cannot be directly inferred by RS since they have little influence on the outgoing radiance observed from space. The estimation of functional parameters such as V cmax or m is complicated even from field observations. These parameters must be estimated from leaf level measurements, which are time and labor demanding and not always feasible. Then, even if https://doi.org/10.5194/bg-2019-501 Preprint. Discussion started: 10 February 2020 c Author(s) 2020. CC BY 4.0 License.
feasible, leaf-level estimates need to be up-scaled to the ecosystem level considering the intra and inter-specific variability of the vegetation around the EC tower (Sprintsin et al., 2012), increasing their uncertainty. For this reason, process-based model 425 inversion has been used to estimate these parameters from fluxes in EC stations (Zheng et al., 2017;Reichstein et al., 2003). However, coherent estimation or scaling of these parameters strongly depends on the description of canopy structure (Sprintsin et al., 2012); so that lack or uncertain knowledge on the biophysical properties of vegetation around the tower can result in uncertain prediction of fluxes (Wang et al., 2019) and consequently in uncertain estimates of the functional traits.
Alternatively, hyperspectral imagery is an independent source of information on canopy structure, biochemistry (Ustin and 430 Gamon, 2010;Schaepman et al., 2009) , and partly to function (Serbin et al., 2015); its numerous and fine bands are sensitive to narrow spectral features and allow solving over-determined systems, which increases the accuracy of the retrievals (Goetz, 2009). The present manuscript demonstrates that the combination of hyperspectral and EC data in models accurately representing radiative transfer, energy balance and photosynthesis allows a coherent estimation of both biophysical and functional traits. This is possible because the multiple-constraint imposed by the combination of RS and EC observations is 435 more robust against model errors than the constraint imposed by RS data only (Pacheco-Labrador et al., 2019b). Therefore, this combination might be also more suitable than the mere ingestion of existing RS products of biophysical parameters such as LAI, which in some cases include ecosystem-dependent uncertainties (Wang et al., 2019). In a scenario where both EC flux stations networks are developing quickly (e.g. FLUXNET, ICOS, OzFLUX, and NEON) and hyperspectral satellite missions are increasing (e.g. DESIS, PRISMA, HyspIRI, EnMAP, etc.), our work calls for a combination of hyperspectral 440 and EC data together with models that accurately represent radiative transfer, energy balance and photosynthesis (RTM-SVAT). The characterization of functional (and biophysical) traits in these networks could be later used to inform, constrain and evaluate TBM, to inform TMB input uncertainties (Prichard et al., 2019), or to benchmark estimates provided by different models and methods (e.g., Luo et al., (2019), Croft et al., (2017), Zhou et al., (2014) or Xie et al., (2018)). Also, when comprehensive enough, they could be predicted at global scale using data-driven approaches  or RS 445 (Serbin et al., 2015).
The application of the multiple-constrained inversion of a coupled SVAT-RTM model at ecosystem scale implies facing sources of error and uncertainty that were not present in previous works carried out at close range over grassland plots (Pacheco-Labrador et al., 2020;Pacheco-Labrador et al., 2019). These can be summarized as: 1) the spatial mismatch between optical, thermal and EC footprints, 2) uncertainties in the estimation of GPP and EC energy balance closure, 3) the 450 atmospheric correction of airborne imagery and the addition of uncertainties in the EnMAP End-to-End simulator, and 4) errors induced by the representation of a Mediterranean tree-grass ecosystem featuring species with very different function (grasses and trees) when using a 1-D model to represent light interaction and photosynthesis. In some cases, we have applied measures to correct some of these uncertainties that could be generally used. Despite of these uncertainties, the estimated parameters are coherent with field data and pattern-oriented model evaluation demonstrated the robustness of the proposed 455 approach. This can be explained since the combination of multiple constraints carrying different and relevant information on the parameters to be estimated is able to counterweight the effect of uncertainties both in model and observations (Keenan et https://doi.org/10.5194/bg-2019-501 Preprint. Discussion started: 10 February 2020 c Author(s) 2020. CC BY 4.0 License. al., 2011;Wutzler and Carvalhais, 2014). An evaluation of different sets of constraints of the SCOPE showed the strong constrain of GPP on functional traits, and its capability to prevent deviations of biophysical parameters due to uncertainties (Pacheco-Labrador et al., 2019). The robustness of the method proposed is promising for its application in different EC 460 stations monitoring different ecosystems, which is still to be tested.
The first abovementioned source of uncertainty is the spatial mismatch between optical and an EC footprints, which is inherent to the combination of RS and EC data. To minimize this mismatch, we used footprint climatology approach (Kljun et al., 2015); which can be applied when pixel size is relatively small compared with EC footprint size. However, to simplify modeling and assimilation, we averaged R λ convolved with half-hourly footprint PDF to provide a representative spectral 465 signature of the area generating the fluxes. Pacheco-Labrador et al., (2017a) found that the spatial variability of the EC footprint induced limited changes in the spectral signal integrated from CASI imagery in this site; and that this variability had little impact in the modeling of GPP. For this reason the spectral signal integrated the week around the image acquisition was considered representative enough of the fluxes source; whereas the spectral variability induced by half-hourly EC footprints was used to weight the cost function in inversion. Nonetheless, this approach might not be advisable for sites 470 where the properties of the source area strongly vary with footprint displacement, and its suitability should be evaluated at each EC station. A second mismatch exists between the footprints of the T ↑ hemispherical sensor and the EC / RS footprints (Marcolla and Cescatti, 2018). Radiation sensors in EC towers are supposed to be representative of the ecosystem monitored, but this is a challenge in many EC sites (Leuning et al., 2012). This mismatch is especially problematic in ecosystems that are heterogeneous and / or that feature scattered 3D volumes; in these cases the contribution of sunlit and shaded vegetation 475 and different plant types to the hemispherical sensors can vary through the day as a function of sensor location and sun position. These contributions might not be coherent with the radiative regime of the EC footprint in all the cases (Marcolla and Cescatti, 2018), which induces part of the error in the closure of the energy balance (Leuning et al., 2012). In this work, we did not correct for this mismatch, which would require detailed 3D information of the tower environment and of the location of the sensor. Footprint differences might have induced the overestimation of T ↑ (Fig. 6b,g) and increased ε EB since 480 the contributions of colder tree crowns might be enhanced by the closer distance to the sensor, and because cooler shaded grass patches are not represented by the model.
A second source of uncertainty is related to the EC data. We attempted to reduce the inherent noise to EC GPP data by smoothing time series (Damm et al., 2010) around the RS observation time. This was relevant for the estimation of biophysical parameters in Step#1, which relies on a single GPP value. These parameters are strongly constrained by GPP via 485 APAR (Pacheco-Labrador et al., 2019), and unlike V cmax , are not refined in Step#2. For this reason, we also aggregated GPP data 1.5 h around each overpass. Results suggest that estimated biophysical parameters are robust or at least consistent within the dataset. Additionally, it must be considered that the energy balance closure gap in real EC datasets is relatively high (median ε EB ~21.5 % of R n in our case) compared to the gap achieved by models, in this case senSCOPE (median ε EB ~0.12 % of R n ). Large observational ε EB hinders the direct comparison of predicted and observed water and energy fluxes. 490 For example Fig.5 shows disagreement between these fluxes, resulting λE overestimated but in the dry season and H biased https://doi.org/10.5194/bg-2019-501 Preprint. Discussion started: 10 February 2020 c Author(s) 2020. CC BY 4.0 License.
(slope ~ 0.45), especially during spring. However, these differences strongly correlated with ε EB (Fig. A1) which suggest that this disagreement might be produced by observation uncertainties rather than by biased predictions. The evaluation of EF provides a more independent comparison of these fluxes and is closer to the 1:1 line than the original fluxes; but, still deviations from this reference are observed. So far, different inversions of SCOPE have relied either on remote thermal 495 imagery (Bayat et al., 2018), on EC TIR irradiance (this work) or on λE (Dutta et al., 2019) to constrain functional parameters related with transpiration such as the Ball-Berry sensitivity parameter; but no comparison of these constraints has been yet carried out. Some of the advantages and disadvantages of each variable are clear: TIR radiation measured in EC towers from hemispherical diffusors offer high temporal frequency, but the radiometric footprint does not match the EC footprint. The use of TIR imagery would allow the use of footprint climatology models minimizing the mismathc, however 500 diel series of TIR imagery feature low spatial resolution. On the other hand λE and H match the GPP footprint, but they present large uncertainties related with the lack of closure of the energy balance. The impact of these uncertainties in the inversion of SCOPE and similar models is still unknown.
RS data also include uncertainties related to the instrumentation, the atmospheric correction and directional effects; the last partly increased by the wide FOV of CASI (~40 º). We tried to minimize directional effects mosaicking different overpasses 505 when available. Space-borne hyperspectral imagery do not often provide multi-angular data which makes unlikely mosaicking or correction of angular effects; but at the same time it also presents cross-track field of views much narrower than our airborne sensor and therefore a lower range of observation angles (Guanter et al., 2015;Ungar et al., 2003;Lee et al., 2015). In the case of this study, directional effects are also expected to be minimized by the convolution of imagery with footprint climatology; which enhances the contribution of those pixels closer to the EC towers, where the flight tracks were 510 centered and therefore observed with lower zenith angles. Compared to the airborne imagery, different and even larger radiometric uncertainties (e.g., no empirical line correction is applied) and lower spatial and spectral resolutions are expected form space sensors. The applicability of the method had to be tested on a dataset comparable to space-borne hyperspectral data. Since no such imagery was available at the site; we used the EnMAP End-to-End simulator to produce synthetic images including spatial, instrumental and atmospheric uncertainties representative of this sensor (Segl et al., 2012). These 515 uncertainties were added to the ones already existing in the CASI dataset, producing an even rougher test for the method.
Despite of lower resolutions and larger uncertainties of the EnMAP images, results were consistent between both sensors. This suggests that space-borne hyperspectral sensors similar to EnMAP could be successfully used to estimate vegetation functional traits with the approach evaluated in this manuscript.
The fourth source of uncertainty is model error. The Mediterranean tree-grass ecosystem under study presents features not 520 well represented by senSCOPE and other models. Pacheco-Labrador et al, (2019) inverted SCOPE in a Mediterranean grassland of this site and found that the model was not suitable to characterize radiative transfer and physiological processes in such ecosystem due to the presence of large fractions of senescent leaves. This problem is partly corrected by senSCOPE, a modified version of SCOPE that separately represents green and senescent leaves. However, the problem is not completely solved. In that work senSCOPE improved the retrieval of C ab via inversion, but still overestimated NIR R λ and underestimated LAI suggesting that the optical properties of senescent material were not adequately represented. This conclusion was in agreement with those of of Melendo-Vega et al., (2018) in the same site. Moreover, we must consider that in the present work, senSCOPE (a 1-D model) is used to represent a Mediterranean tree-grass ecosystem featuring 1) a strong 3-D structure, 2) a mixture of species (grasses and trees) with very different function and ecological strategies, and 3) a grassland layer featuring large fractions of senescent leaves and death standing material, whose optical properties are 530 misrepresented. These conditions portray a harsh test for the inversion method; however, model parameters seem to be robustly estimated, are coherent with ancillary observations and between airborne and synthetic EnMAP imagery. Still, the effect of some of these uncertainties can be analyzed. For example, angular effects can be important in savanna-like ecosystems featuring a strong geometrical scattering component. Melendo-Vega et al, coupled the 1-D RTM PROSAIL (Jacquemoud et al., 1996;Verhoef, 1984) with the 3-D RTM FLIGHT (North, 1996) to represent tree-grass ecosystems and 535 compared this model with PROSAIL, showing that for a low tree cover, differences in predicted R λ minimized at nadir.
Consequently, the reduction of off-nadir view directional effects produced by footprint integration and mosaicking might have contributed to reduce the impact of 3-D structure during inversion. However, the 1-D model used does not reproduce tree crown shading on the grassland; which should affect R λ and APAR during the day course. Also, the misrepresentation of optical properties of the senescent material biases the estimates of some parameters in an attempt to compensate low NIR 540 scattering. At grass plot scale LAI, NIR R λ were underestimated and f green was overestimated (Pacheco-Labrador et al., 2020); whereas at ecosystem scale LAI and NIR R λ were not so strongly biased, and f green was underestimated in the dry period. In both cases C s and C dm were large, especially in the dry season, since both reduce NIR R λ . At ecosystem scale tree crowns increase NIR R λ in all the seasons (Pacheco-Labrador et al., 2017a) which might have prevented the underestimation of LAI. C ab and V cmax estimates and their relationship with nitrogen are consistent with field data. Retrieved V cmax is lower than 545 independent estimates based on nitrogen concentration and literature not specific of semi-arid Mediterranean ecosystems.
These estimates might present values too high for the study site. For example, the N mass -V cmax curve reported in Feng and Dietze (2013) and used to generate the evaluation dataset predict higher V cmax than values estimated in the grass plots in the same ecosystem (Pacheco-Labrador et al., 2019). An alternative hypothesis is related to the fact that tree crowns reduce down-welling irradiance impinging the grassland during the day, which is not adequately represented by the 1-D RTM. 550 Overestimation of APAR (and therefore GPP) might have been compensated decreasing the photosynthetic leaf surface (f green ), or V cmax to fit observed GPP. Pattern-oriented model evaluation of m and uWUE against Δ 13 C provided expected relationships when comparing spring m and uWUE estimates with Δ 13 C, a seasonally-integrated indicator of plant water use efficiency. In the dry season m estimates showed high values for both sensors (CASI and synthetic EnMAP), suggesting low uWUE. Wolf et al., (2006) also estimated large m values in a grassland during the dry season, and related them with 555 decoupled simulation of transpiration with other processes (which is not the case of the model used here), and to assumptions about leaf thermal properties. Contrarily, Bayat et al., (2018) did not found such large values during summer in a different grassland. However Bayat et al., (2018) inverted SCOPE model where all the leaf area can transpire; whereas in senSCOPE transpiration is limited to the green fraction. Looking at the fluxes can be noticed that in the dry season, predicted H is in https://doi.org/10.5194/bg-2019-501 Preprint. Discussion started: 10 February 2020 c Author(s) 2020. CC BY 4.0 License. agreement with observations, whereas λE is overestimated (Fig. 6); these discrepancies are slightly biased respect to EC ε EB 560 (Fig. A1). We hypothesize that m estimates in the dry season might not be realistic: Overestimation of T ↑ might have demanded large transpiration during summer time; when low GPP and f green (also underestimated) might have reduced the impact of m on the cost function, allowing large values with little effect on T ↑ . It must be also considered that the Ball-Berry model implemented in SCOPE and derived models relies on relative humidity; this approach does not well represent response leaf responses of both conductance and intracellular CO 2 concentration to stressing humidity conditions, and 565 alternative models have been proposed (Leuning, 1990(Leuning, , 1995. These models might also prevent the overestimation of m during stress conditions. Further work is needed to understand the challenges in the seasonal retrieval and validation of m estimates form RS. The previous analyses showed the need and the added value of evaluating functional trait estimates in the context of the development of new RS products of vegetation function and functional traits. Since direct observations of these traits would 570 not be as usual as other biophysical variables that are easier to measure in the field; new methods such as pattern oriented evaluation should be developed in parallel. Thorough evaluation contributes to understand the effect of different sources of uncertainty, model structure and to identify possibilities of improvement; and therefore it will be fundamental to support the development of the next generation of RS products. To the three challenges to overcome proposed by Schimel et al., (2019) we should also add the validation and evaluation of the new functional trait estimates as an additional challenge to overcome. 575

Conclusions
This work proves the potential of the new generation of hyperspectral space missions in combination with EC networks to characterize the spatiotemporal distribution of vegetation functional parameters relevant for the prediction of water and carbon fluxes. The method proposed has proven robust to different sources of uncertainty; some of which were minimized with different approaches (e.g., footprint climatology convolution, or the use of senSCOPE model). Repeated acquisition of 580 images over the same EC stations can provide not only seasonal variability of key functional and biophysical parameters, but also, it can characterize relationships between biophysical and functional parameters and other variables, such as those of C ab or V cmax with nitrogen, m or uWUE with Δ 13 C, or r ss with SM p . Spatiotemporal characterization of functional traits as well as these relationships could useful for validation, modeling or inversion of TBM and other models. Results also show that biases in the prediction of λE and H strongly correlate with observational energy balance closure error. Pros and cons of the 585 use of energy fluxes or TIR radiance from proximal or remote sensors to constrain leaf water exchange are discussed; further research is needed to determine which of these constraints is the most suitable for the estimation of functional parameters in

595
Code availability senSCOPE code and further developments, as well as the code for the multiple constraint inversion of the model are publicly available at https://github.com/JavierPachecoLabrador/senSCOPE.

Data availability
All data included in this study are available upon request by contact with the corresponding author.