Regional detection of canopy nitrogen in Mediterranean forests using the spaceborne MERIS Terrestrial Chlorophyll Index

Canopy nitrogen (N) concentration and content are linked to several vegetation processes and canopy N 10 concentration is a state variable in global vegetation models with coupled carbon (C) and N cycles. While there is ample C data available to constrain the models, widespread N data are lacking. Remote sensing and vegetation indices have been used to detect canopy N concentration and canopy N content at the local scale in grasslands and forests. In this paper we conducted a regional case-study analysis investigating the relationship between the Medium Resolution Imaging Spectrometer (MERIS) Terrestrial Chlorophyll Index (MTCI) time series from ESA ENVISAT at 1 km spatial resolution and both canopy N 15 concentration (%N) and canopy N content (g m) from a Mediterranean forests inventory in the region of Catalonia, NE of Spain. The relationships between the datasets were studied after resampling both datasets to lower spatial resolutions (20 km, 15 km, 10 km and 5 km) and at the initial higher spatial resolution of 1 km. The results at the higher spatial resolution yielded significant relationships between MTCI and both canopy N concentration and content, r = 0.32 and r = 0.17, respectively. We also investigated these relationships per plant functional type. While the relationship between MTCI and canopy N 20 concentration was strongest for deciduous broadleaf and mixed plots (r = 0.25 and r = 0.47, respectively), the relationship between MTCI and canopy N content was strongest for evergreen needleleaf trees (r = 0.20). At the species level, canopy N concentration was strongly related to MTCI for European Beech plots (r = 0.71). These results present a new perspective on the application of MTCI time series for canopy N detection, ultimately leading towards the generation of canopy N maps that can be used to constrain global vegetation models. 25


Introduction
Canopy nitrogen (N) concentration is an essential state variable in regional (Ollinger and Smith, 2005) and global vegetation models including both the carbon (C) and the N cycles (Zaehle and Friend, 2010;Smith et al., 2014).This variable has been linked to several vegetation traits and processes, at the leaf and canopy levels.At the leaf level, leaf N concentration, which represents the leaf N status expressed as a percentage of leaf dry mass (%N, N g 100g -1 DM), has been related to photosynthetic capacity (Evans, 1989;Reich et al., 1995;Reich et al., 1997;Reich et al., 1999;Wright et al., 2004), specific leaf area, leaf life span (Reich et al., 1999;Wright et al., 2004) and light use efficiency (Kergoat et al., 2008).Leaf N concentration expressed on a leaf area basis, also called leaf N content (N g m -2 ) has also been linked with chlorophyll content, Rubisco content (Evans, 1989) and photosynthetic capacity (Evans, 1989;Reich et al., 1995).At stand scale, canopy nitrogen concentration, which represents the leaf N concentration averaged over the stand canopy, has also been found to correlate with above ground Net Primary Productivity (NPP) (Reich, 2012), while canopy N content has been linked with the canopy light use efficiency (Green et al., 2003).
Given their links to many vegetation processes, leaf and canopy N variables could be used to constrain N cycle modules in global vegetation models.At the global scale, ample data is available to constrain models for the C cycle; however, data to constrain the N cycle are limited.Currently, canopy N data is not widely available and canopy N sampling campaigns are timeconsuming and thus expensive tasks.Moreover, upscaling from local sampling campaign measurements represents an additional limitation.In this perspective, local, regional or even global remotely sensed canopy N estimates will be a valuable addition, enabling us to collect information in a less time intensive and expensive manner than traditional on-field sampling campaigns.Such near global canopy N estimates will be beneficial as input in global vegetation models or to calibrate and validate these models.
Currently, different remote sensing techniques have been applied to detect canopy N in terrestrial vegetation.Imaging spectrometry from either airborne or spaceborne sensors coupled with different analysis methods, including partial least squares regression (PLS), continuum removal, spectral unmixing or vegetation indices, has proven efficient in improving N sensing capabilities at the local scale (Smith et al., 2003;Ollinger et al., 2008;Huber et al., 2008;Martin et al., 2008;Schlerf et al., 2010;Wang et al., 2016).
Among other techniques, ratios or normalized differences of reflectance bands in the Red and Near Infrared (NIR) regions of the spectrum, the so called vegetation indices (VI) (Glenn et al., 2008), are one of the most straightforward methods for canopy N detection.Combined with in situ hyperspectral devices, vegetation indices have been extensively used for leaf or canopy N detection in agricultural systems (Peñuelas et al., 1994;Filella et al., 1995;Hansen and Schjoerring, 2003;Tian et al., 2011;Schlemmer et al., 2013;Li et al., 2014).Vegetation indices have also been applied to airborne or spaceborne acquired imagery in natural environments (Ramoelo et al., 2012;Wang et al., 2016).
A particular vegetation index, the MERIS Terrestrial Chlorophyll Index (MTCI) has been proposed for detecting canopy N (Clevers and Gitelson, 2013).MTCI was originally computed from three reflectance bands from the Medium Resolution Imaging Spectrometer (MERIS) aboard the European Space Agency (ESA) ENVISAT satellite at a spatial resolution of 1 km.
However, it can also be obtained from other sensors' reflectance data and a similar product will be available from the ESA Sentinel-2 satellite mission (Drusch et al., 2012).It was first developed to estimate chlorophyll content (Dash and Curran, 2004Curran, , 2007)).Since then, other applications of this index have been described, among which the possibility to estimate Gross Primary Productivity (GPP) from natural (Harris andDash, 2010, 2011;Boyd et al., 2012) and cultivated lands (Peng and Gitelson, 2011).Furthermore, MTCI has been used to discriminate between C3 and C4 grasses (Foody and Dash, 2007) and to monitor vegetation phenology at the sub-regional (Boyd et al., 2011) and continental scales (Rodriguez-Galiano et al., 2015;Crabbe et al., 2016).Regarding canopy N detection, most applications were aimed at agricultural crops using MTCI values computed from in situ hyperspectral reflectance data (Tian et al., 2011;Clevers and Gitelson, 2013;Li et al., 2014).A few were directed towards sensing N concentration in natural environments using airborne data, e.g. in temperate forests (Wang et al., 2016), or spaceborne data, for example in grasslands (Ramoelo et al., 2012;Ullah et al., 2012) or sub-tropical forests (Cho et al., 2013).
Remote detection of foliage N status has been extensively studied at the leaf scale (Hansen and Schjoerring, 2003;Ferwerda et al., 2005;Pacheco-Labrador et al., 2014;Li et al., 2014) and a few studies have investigated the processes underlying the relationships between vegetation indices and canopy N. Detection of foliage N status with vegetation indices is attributed to the strong link between foliar nitrogen and chlorophyll content (Schlemmer et al., 2013) and is often based on the NIR and red-edge region of the spectrum, hence similar to the ones used for chlorophyll detection (Filella and Penuelas, 1994;Dash and Curran, 2004;Clevers and Gitelson, 2013).It is not clear however how this translates at the canopy level as the mechanisms possibly modifying the remote detection of foliage N status at the canopy scale are still not clearly understood.High correlation between canopy N and both NIR reflectance and albedo has been reported in boreal forests (Ollinger et al., 2008).However, the mechanism behind these findings is still controversial.Knyazikhin et al. (2013) argued that the observed correlation solely resulted from canopy structural differences between broad and needleleaf forests and was thus spurious.Other authors, although agreeing that canopy structure was a confounding factor to account for, stated that the NIRcanopy N relationship was not necessarily spurious as plant traits have been known to covary along the leaf economic spectrum (Wright et al., 2004;Ollinger et al., 2013;Townsend et al., 2013).
In this context, there are several knowledge gaps that we would like to address in this paper.First, although 1 km spatial resolution spaceborne MTCI time series are available from the ESA, MTCI has mainly been employed to detect canopy N in agricultural applications with in situ devices and rarely in a broader range of natural ecosystems and scales using spaceborne data.Due to its almost global coverage, MTCI time series could be applied to estimate canopy N at a larger scale.Moreover, to our knowledge limited research has been conducted to sense canopy N in Mediterranean ecosystems (Serrano et al., 2002) and even more so in Mediterranean forests.In addition, although in a temperate forest the reflectance spectrum of individual plant functional types (PFT) has been shown to be different (Wang et al., 2016), the relationship between MTCI and canopy N has seldom been studied and compared between PFTs.Moreover, investigating the influence of PFTs on this relationship might give further insight into the influence of structural effects in canopy N detection.Finally, the difference between sensing canopy N concentration (N[%], %N) and canopy N content (N[area], g m -2 ) has rarely been investigated.The relationship between MTCI and both of these variables has been studied separately (Clevers and Gitelson, 2013;Wang et al., 2016) The objective of our study is thus to investigate the relationship between the spaceborne MTCI remote sensing product and canopy N in Mediterranean forests at the regional scale.More specifically, the relationships between MTCI and both canopy N concentration and canopy N content are investigated and compared.We then also examine these relationships per PFT and at the species level.
Detection of canopy N is often limited to local scale studies due to the spatial restrictions associated with N data acquisition in the field and treatment of high spatial resolution remote sensing imagery with limited spatial coverage (Lepine et al., 2016).
Our case-study exploits the broadly and readily available MTCI time series at 1 km spatial resolution from the ESA ENVISAT mission and combines it with canopy N data, both concentration and content, from 1075 forest plots measured between 1988 and 2001 by the Catalonian National Forest Inventory (Gracia et al., 2004).First, we develop a methodology to overcome the time discrepancy between our two sets of data.Next, both data sets are resampled to the same, lower, spatial resolutions in order to overcome the initial scale discrepancy between MTCI spatial resolution (1 km) and the size of the forest plots (6 m).
Subsequently, we analyse the relationship between MTCI and both canopy N concentration and canopy N content variables, both at the resampled and initial spatial resolutions.These relationships are then stratified according to the plots PFT.The results are presented and discussed.Finally, we address the implications for future research and draw a conclusion.

Study area
Our study area corresponds to the region of Catalonia (Fig. 1 (Gracia et al., 2004).
The forest plots (Fig. 1) had a minimum diameter of 6 m, which varied depending on the tree density in order to include between 15 and 25 trees with a diameter at breast height of at least 5 cm.The plots were investigated for canopy N concentration (N[%], %N) defined as g of N per 100 g of dry leaf matter.Each leaf sample was constituted by the leaves of at least three different trees of the dominant tree species in the canopy.The species dominance was determined by the tallest individuals.
The leaves were collected from the upper central part of the crown using extensible loppers.All foliar cohorts present in the canopy were included in the leaf sample.
The samples were dried and then ground using a Braun Mikrodismembrator-U (B.Braun Biotech International, Melsungen, Germany).They were analysed for foliar N concentration using the combustion technique coupled to gas chromatography using a Thermo Electron Gas Chromatograph (model NA 2100, CE Instruments-Thermo Electron, Milan, Italy) (Gracia et al., 2004).
Along with the canopy N[%] data, we used foliar biomass data (g m -2 ) acquired during the same forest inventory (n = 2286).
The foliar biomass data were obtained for each plot from allometric equations relating the diameter of the branches to the leaves dry weight (Gracia et al. (2004), Fig. A1).These foliar biomass data were used to calculate canopy N content (N[area], g m -2 ) for each plot following Eq.( 1): where   [𝑎𝑟𝑒𝑎] is the canopy N content (N g per square meter of ground area, g m -2 ),   [%] is the canopy N concentration (%N) and fbiom is the foliar biomass (N g per square meter of ground area, g m -2 ).
Catalonian forests include both deciduous and evergreen broadleaf as well as evergreen needleleaf tree species.These three PFTs are referred to as Deciduous Broadleaf Forest (DBF), Evergreen Broadleaf Forest (EBF) and Evergreen Needleleaf Forest (ENF), respectively.The main tree species are Pinus halepensis Mill., Pinus sylvestris L., Quercus ilex L., Pinus uncinata Ramond ex DC., Pinus nigra J.F. Arnold, Quercus suber L., Quercus cerrioides Willk.& Costa., Quercus petraea Liebl.and Fagus sylvatica L. These species accounted for 92% of the sampled forest plots.The 15 tree species included in this analysis are listed in Table 1.Plots with a rare dominant tree species, i.e. species that were detected in only one single plot, were excluded from the analysis.This applied to plots with these dominant species: Abies alba Mill., Fraxinus augustifolia Vahl, Fraxinus excelsior L., Pinus radiata D. Don, Populus nigra L., Populus tremula L., Quercus robur L.
A proportion of 96% of the plots included in this analysis had a single dominant tree species.The plots with codominant species (n = 30) were sampled several (up to two) times, each time for each of the codominant tree species found on the plots.
After separate measurements, the obtained foliar N concentration values were averaged in order to obtain a single canopy N[%] value for each plot with several codominant species.Among these 30 plots with codominant species, 16 plots had codominant species from different PFT.Their PFT is thus labelled as mixed while the plots with several codominant species from the same PFT are labelled according to their PFTs.

MTCI product
The MERIS Terrestrial Chlorophyll Index (MTCI) was first developed to estimate chlorophyll content in canopies.Its calculation, presented in Eq. ( 2), is based on three reflectance bands, located around the red edge point (REP) (Dash and Curran, 2004): where  8 ,  9 and  10 represent the 8 th , 9 th and 10 th bands of MERIS, respectively.Following MERIS standard bands settings, the centres of the bands were located at 681.25 nm, 708.75 nm and 753.75 nm on the electromagnetic spectrum.
While the ESA ENVISAT satellite mission producing MERIS data came to an end in 2012, MERIS products and MTCI in particular are still relevant because the new ESA Sentinel-2 satellite mission has improved band settings compared to those of MERIS and increased the spatial resolution to 20 m (Drusch et al., 2012).The Sentinel-2 mission will also release a chlorophyll product that will continue the time series already available for MTCI.In this study, we put emphasis on ENVISAT-MERIS as our field data are closer to the MERIS acquisition period.
MTCI has been described to be sensitive to high chlorophyll content while presenting low sensitivity to soil brightness.MTCI 10 year time series are available almost globally and are already corrected for atmospheric influences and cloud cover (Curran and Dash, 2005).
MTCI level 3 imagery was obtained from the NERC Earth Observation Data Centre (NEODC, 2015) for the region of Catalonia between 2002 and 2012.The original data were provided by the European Space Agency and then processed by Airbus Defence and Space.The original data, following ENVISAT specifications, have a revisit time of three days and the processed imagery is available as an either weekly or monthly averaged product.The spatial resolution of the processed data is approximately 1 km.One MTCI imagery product covering the entire study area was obtained for every month between June 2002 and March 2012, except for October 2003, when no valid product was available.

Methodology to link canopy N data to MTCI values
There is a discrepancy between the timing of the ground truth sampling and the satellite image acquisition period.While the plot sampling campaigns were carried out between 1988 and 2001, the ENVISAT satellite mission was launched in 2002 and ended in 2012.To overcome the discrepancy, MTCI images were averaged by month over the 10 years of the satellite mission period.This process yielded twelve MTCI averaged images, one for each month.The averaged MTCI images were then linked to the forest plots based on the forest plot coordinates and sampling month, as the exact sampling date was known for each plot.The period between the 1 st of June and the 31 st of October was determined to be the growing season after a pre-analysis, where we studied yearly temporal variation of MTCI in several locations and forest types in Catalonia.This extended period was chosen to encompass the different vegetation phenology types corresponding to the contrasted climate conditions in this  , 2010).This map comprises 22 land cover classes and has a spatial resolution of 300 m.Using this map, we excluded forest plots that had undergone a land cover change since the sampling period and did not have a natural vegetation cover any more at the time of remote sensing image acquisition.To do so, the landcover map was first resampled to a spatial resolution of 1 km to be in accordance with MTCI spatial resolution.Then, the plots located on land area classified as either rainfed cropland, mosaic between croplands and natural vegetation, sparse vegetation or artificial surfaces were excluded from the analysis.

Relationship between MTCI and canopy N data at lower spatial resolution
In a first step, the relationships between MTCI and canopy N data values were investigated after resampling both datasets to the same, lower, spatial resolution.This was done because of the initial difference in support size between MTCI spatial resolution and the forest plots size (i.e. 1 km and 6 m, respectively).This enabled us to investigate the relationships between MTCI and canopy N data independently of differences in initial support size.
The Globcover 2009 land cover map was used to exclude from the computation the MTCI pixels located on land surface without natural vegetation cover.As for the forest plots, MTCI pixels whose land cover class corresponded to rainfed cropland, mosaic between croplands and natural vegetation, sparse vegetation or artificial surfaces were excluded from the upscaling analysis.The monthly averaged MTCI images obtained previously were resampled successively to 5 km, 10 km, 15 km, and 20 km.Forest plots data were then averaged over the newly obtained pixel.The relationship between the resampled MTCI values and canopy N data was analysed using linear regression.

Relationship between MTCI and canopy N data at initial higher spatial resolution
In a second step, the relationships between MTCI and canopy N data, both canopy N[%] and canopy N[area], were examined at the original spatial resolution of 1 km.This allowed us to investigate the influence of PFT and species on the relationships as this information was lost in the resampling process.The relationships between MTCI and canopy N at 1 km-spatial resolution were analysed with linear regression for the whole dataset, for each PFT separately as well as for individual species.

Statistical analysis
After applying the selection criteria explained above, i. 3 Results

Descriptive statistics
Descriptive statistical analysis of canopy N[%], canopy N[area] and foliar biomass were performed for each tree species included in the dataset (Table 1).The four most abundant species (Pinus halepensis, Pinus sylvestris, Quercus ilex and Pinus uncinata) dominated 667 plots i.e. almost 80% of the plots.The cumulated abundance percentages of ENF, EBF and DBF species were equal to 66 %, 22 % and 9 %, respectively.From this data, it is clear that the forests plots were mainly dominated by ENF species.On average, Pinus uncinata plots had the highest biomass values while Quercus suber plots showed the lowest mean value for this variable.Descriptive statistics were also analysed by PFT.The mean canopy N[%] was lowest for ENF species, 0.97 %N, and highest for DBF trees, 2.17 %N (Fig. 2a).Canopy N[%] value ranges were equal to 1.91 %N, 2.06 %N , 1.68 %N and 1.42 %N for DBF, EBF, ENF and mixed plots, respectively.The canopy N[area] statistics were analysed by PFT as well (Fig. 2b) and the averaged canopy N[area] values ranged from 1.82 g m -2 to 4.61 g m -2 .A Pearson correlation matrix (Fig. 3) was computed between the variables for the whole dataset and while the correlation between each pair of variables was significant, the relationship between canopy N[area] and foliar biomass was strongest (r = 0.88).This matrix also shows distribution histograms of the three variables.As canopy N[%] and canopy N[area] distributions are positively skewed, a logarithmic transformation was applied to these variables to fulfil linear model assumptions.Correlation matrices for each DBF, EBF and ENF plots are presented in the Appendix (Fig. A2 -4).

Relationship between MTCI and canopy N data at lower spatial resolution
The relationships between MTCI and both canopy N[%] and canopy N[area] were studied after resampling both datasets to the same, lower, spatial resolution.This was done to investigate the relationship between MTCI and canopy N data independently of differences in support size.The results showed that the relationship between MTCI and canopy N[%] was always stronger than the relationship for MTCI and canopy N[area] for each resampling factor.Moreover, the relationships between MTCI and either canopy N[%] or canopy N[area] were all highly significant (p<0.000).The r 2 values of the relationship between MTCI and canopy N[%] were equal to 0.33, 0.37, 0.35 and 0.43 for 5 km, 10 km, 15 km and 20 km resampled spatial resolution, respectively.The r 2 values of the relationship between MTCI and canopy N[area] were equal to 0.17, 0.20, 0.20 and 0.17

3.3
Relationship between MTCI and canopy N data at higher spatial resolution

Relationship between MTCI and canopy N concentration
The relationships between MTCI and canopy N data were studied after resampling both datasets to the same, lower, spatial resolution.The results showed that the linear regression between MTCI and canopy N[%] for the whole dataset (n = 846) was highly significant (p<0.000) and had a r 2 value of 0.32 (Table 2, Fig. 5a).The relationship between MTCI and Canopy N[%] was also investigated for each PFT individually (Fig. 5b-e).For DBF plots, the relationship between MTCI and canopy N[%] had an r 2 value of 0.25 (n = 80) and was significant.However, although significant, the r 2 of the relationship between MTCI and canopy N[%] for EBF and ENF plots were lower and equal to 0.03 (n = 186) and 0.10 (n = 564), respectively.
The relationship between MTCI and canopy N[%] was also significant for one individual species, Fagus sylvatica.The proportion of explained variance for this species was equal to 0.71 (n = 15).This result, although obtained on a restricted number of plots, showed that the significant relationships between MTCI and canopy N[%] not only existed when all DBF plots were included but also held for one individual DBF species.

Relationship between MTCI and canopy N content
Significant relationships between MTCI and canopy N[area] were found for the whole dataset as for EBF and ENF plots (Table 3).The scatterplots between MTCI and canopy N[area] are presented in Figure 6.The proportion of explained variance was higher for ENF plots compared to the other PFTs and compared to the overall relationship across all plots.The relationship between MTCI and canopy N[area] was also investigated for 10 individual species and one of them showed significant relationships: Quercus ilex.

Discussion
Our aim was to explore the relationship between the MTCI vegetation index and both canopy N[%] and canopy N[area] in Mediterranean forests at the regional scale in Catalonia, north eastern of Spain.This was done by using the ESA spaceborne MTCI remote sensing product and canopy N data from a forest inventory.The relationship was first investigated using MTCI and canopy N data resampled to the same, lower, spatial resolution.The relationship was then investigated across all plots and by PFT at MTCI initial spatial resolution of 1 km.

Relationship between MTCI and canopy N data at lower spatial resolution
This pre-analysis was undertaken to study the MTCI-canopy N relationships independently of the discrepancy between MTCI original spatial resolution (1 km) and the size of the forest plots (diameter of 6 m).By resampling both datasets to a lower spatial resolution, the obtained values were less impacted by small-scale variations because they were obtained by averaging several values over a larger area.The results showed that the relationship between MTCI and canopy N data was significant and consistent across all spatial resolutions investigated: 5 km, 10 km, 15 km and 20 km.This shows that, when the influence of the discrepancy between the original datasets was taken into account, MTCI and canopy N data were linked and that the MTCI-canopy N relationship was not strongly affected by the resampled spatial resolution.

4.2
Relationship between MTCI and canopy N data at higher spatial resolution

Canopy N concentration detection
The overall relationship between MTCI and canopy N[%] at 1 km spatial resolution for all the forest plots (n = 846) was significant and the r 2 value was equal to 0.32 (Fig. 5).This result showed that canopy N[%] could be related to MTCI in Mediterranean forests.The performance of the MTCI vegetation index to detect canopy N[%] in Mediterranean vegetation was similar to the results obtained from previous studies using spaceborne MTCI.For example, using MTCI computed from the spaceborne RapidEye sensor at 5 m spatial resolution, it was possible to detect canopy N[%] in grassland savannah and subtropical forest with similar coefficients of determination, r 2 = 0.35 and r 2 = 0.52, respectively (Ramoelo et al., 2012;Cho et al., 2013).More generally, while there is a consensus regarding MTCI ability for in situ leaf or canopy N[%] detection in a variety of crops (Tian et al., 2011;Li et al., 2014), there is no general agreement about MTCI ability for canopy N[%] detection across vegetation and sensor types.For example, MTCI computed from airborne data at 3 m spatial resolution could not be related to were significant for all the PFT taken separately (p-value<0.05).However, a higher proportion of variance was explained for DBF and mixed plots (r 2 = 0.25 and r 2 = 0.47 for DBF and mixed plots, respectively) compared to the other plant functional types (r 2 = 0.10 and r 2 = 0.03 for ENF and EBF trees, respectively) and the relationship between MTCI and canopy N[%] was especially weaker for EBF plots.Therefore, the relationship observed for all the forest plots may be mainly driven by DBF and mixed plots.This result is different from what was observed by Ollinger et al. (2008) in boreal forests, where canopy N[%] was related to NIR reflectance for both broadleaf and needleleaf plots taken separately.Moreover, the results obtained for ENF tree species are surprising as previous studies investigating the relationship between foliar N[%] and in situ measured spectra reported higher r 2 values, r 2 = 0.59 and r 2 = 0.81 in spruce and pine forest, respectively (Stein et al., 2014;Schlerf et al., 2010).
The differences in scale and methodology might explain the divergent results compared to previous findings.Indeed, in our study, the analysis is carried out at a much coarser spatial resolution using spaceborne data compared to the fine spatial scale obtained with in situ devices.Moreover, most of these studies were carried out in temperate forests and studies investigating canopy N[%] detection in Mediterranean regions are scarce.When investigating the relationship between canopy N[%] and MTCI at the species level, we also found that it was significant for Fagus sylvatica plots (r 2 = 0.71).In the literature, the relationship between MTCI and canopy N[%] is often not stratified by PFT or species (Sullivan et al., 2013;Wang et al., 2016).In this study, we showed that investigating this relationship for each PFT taken separately yielded additional insight.Indeed, to our knowledge the difference in explained variance between DBF and other PFTs in MTCI and canopy N[%] relationship has not been observed before.Moreover, we could show that the stronger relationship observed for DBF plots was consistent for Fagus sylvatica plots taken separately.

Canopy N content detection
The relationship between MTCI and canopy N[area], which was obtained by combining canopy N concentration values with biomass data, was significant across all plots (n = 841) (Table 3, Fig. 6).Although the r 2 value was lower for the relationship between MTCI and canopy N[area] (r 2 = 0.17) than for the relationship between MTCI and canopy N[%] (r 2 = 0.32), it is interesting to note that canopy N[area] can be related to spaceborne MTCI as remotely sensed detection of canopy N[area] is rarely investigated in forest environments (Mirik et al., 2005).In comparison, previous studies conducted in grasslands reported higher prediction accuracy e.g. by using spaceborne MTCI at 300 m spatial resolution or a simple ratio-type vegetation index computed from airborne imagery at 1 m spatial resolution, canopy N[area] was detected with r 2 values equal to 0.29 and 0.66, respectively (Mirik et al., 2005;Ullah et al., 2012).
The relationship between MTCI and canopy N[area] was only significant for ENF and EBF plots (Fig. 6b-e), with a higher proportion of explained variance for ENF plots (r 2 = 0.20).However, when this relationship was investigated at the species scale, significant results were found for Quercus ilex (EBF) plots.This is accordance with a previous study examining the detection of canopy N[area] in Quercus ilex trees by MTCI computed from in situ spectra (r 2 = 0.43) (Pacheco-Labrador et al., 2014).

Comparing results obtained for canopy N concentration and canopy N content detection
This analysis highlighted the difference between canopy N expressed as a percentage of leaf dry matter (canopy N[%]) and on an area basis (canopy N[area]) regarding detection by spaceborne MTCI for the different PFTs.Canopy N[%] of DBF and mixed plots showed higher correlation with MTCI compared to EBF and ENF plots while the relationship between canopy N[area] of ENF plots with MTCI was stronger than for any other PFTs.These differences between canopy N[%] and canopy N[area] detection by remote sensing can be related to previous findings showing that canopy N[area] but not canopy N[%] could be detected by MTCI in grassland (Ullah et al., 2012) and by a simple ratio index in heterogeneous rangelands (Mirik et al., 2005) at various spatial scales, 300 m and 1 m, respectively.In the literature, canopy N[%] is more often used to detect N state of foliage in forest while canopy N[area] is regularly employed in grasslands but also in crops (Clevers and Gitelson, 2013;Schlemmer et al., 2013).Our results showed that, for ENF plots, when biomass was accounted for, as in canopy N Other authors, although agreeing that canopy structural properties needed to be accounted for, suggested that a direct biochemical link between canopy N and reflectance data was not necessary to detect canopy N with reflectance data (Ollinger et al., 2013;Townsend et al., 2013).In this context, our analysis showed that the PFTs of the plots had an influence on the MTCI canopy N relationship in a specific type of ecosystem, namely Mediterranean forests.Further analysis is needed to disentangle the influence of structural properties on canopy N detection.

Perspectives for larger scale applications
The methodology applied in this paper is different from the usual methodology implemented to detect canopy N concentration in forests.Remote sensing of N in forest canopies by hyperspectral sensors is often coupled with intensive forest sampling measurements.This method has been effective at detecting canopy N concentration locally in a vast range of environments (Serrano et al., 2002;Smith et al., 2002;Townsend et al., 2003;Ollinger et al., 2008;Wang et al., 2016).Applying this technique at larger scales has already been explored.For example, Martin et al., (2008) compiled 137 field plots data from previous studies in various forest types and investigated the possibility to find a common detection algorithm.However, due to the different treatments required as well as the high spatial resolution (from 3 m to 30 m for Hyspex airborne and Hyperion spaceborne sensors, respectively, Wang et al., 2016;Smith et al., 2003), applying imaging spectroscopy at a broader scale might reveal laborious.A recent study in northern temperate forests explored the effect of spatial resolution on canopy N[%] estimation.The results showed that, although the prediction accuracy was reduced compared to what was achieved using PLS regression at higher spatial resolution, it was still possible to estimate canopy N[%] with r 2 between 0.34 and 0.81 using various vegetation indices computed from MODIS reflectance data at 500 m spatial resolution (Lepine et al., 2016)

Conclusion
In this study, we investigated the relationship between spaceborne MTCI from ENVISAT and both canopy N[%] and canopy N[area] at regional scale in Mediterranean forests.We found significant results across all plots both when the original data were resampled to 5 km, 10 km, 15 km and 20 km and for the original spatial resolution of 1 km.The relationship between MTCI and canopy N data was also significant for some individual PFTs and species.The r 2 values were 0.32 and 0.17 Biogeosciences Discuss., https://doi.org/10.5194/bg-2017-228Manuscript under review for journal Biogeosciences Discussion started: 29 June 2017 c Author(s) 2017.CC BY 3.0 License.
Biogeosciences Discuss., https://doi.org/10.5194/bg-2017-228Manuscript under review for journal Biogeosciences Discussion started: 29 June 2017 c Author(s) 2017.CC BY 3.0 License.region.The forest plots sampled outside of the growing season were excluded from the analysis.The Globcover 2009 land cover map was used to exclude forest plots located on unsuitable land surface.The Globcover map was created by ESA using MERIS reflectance data from 2009 (Bontemps et al., 2011).It was downloaded from the ESA data user elements website (ESA e. plots measured between June 1 st and October 31 st , exclusion of plots with infrequent species and selection based on Globcover 2009, 846 forest plots were available for analysis, including 841 plots with foliar biomass and canopy N content information.Descriptive statistics of canopy N[%], foliar biomass and canopy N[area] were then produced for each of the tree species and PFT included in the analysis.The linear regressions between MTCI and canopy N were then performed for both resampled and non-resampled datasets.Preliminary analysis showed that using a Biogeosciences Discuss., https://doi.org/10.5194/bg-2017-228Manuscript under review for journal Biogeosciences Discussion started: 29 June 2017 c Author(s) 2017.CC BY 3.0 License.natural logarithm transformation (log) of the canopy N variables was necessary to fulfil linear regression model assumptions, namely normality and homogeneity of variance of the residuals.The minimum number of data points needed to carry out the linear regression analysis was fixed at 10.The spatial analyses were done with the PCRaster software (Karssenberg et al., 2010).The statistical analyses were performed in the R environment (R Development Core Team, 2014) and the ggplot2 package was used for the graphics (Wickham, 2009).
canopy N[%] from a mixed temperate forest(Wang et al., 2016).In this context our finding brings new insight into MTCI N[%]     sensing capabilities at a much coarser spatial resolution (1 km) compared to what has been done before.Investigating the influence of the PFTS on the overall relationship highlighted the difference between DBF, EBF and ENF types of vegetation regarding canopy N[%] detection by spaceborne MTCI.The relationships between MTCI and canopy N[%] Biogeosciences Discuss., https://doi.org/10.5194/bg-2017-228Manuscript under review for journal Biogeosciences Discussion started: 29 June 2017 c Author(s) 2017.CC BY 3.0 License.
[area], the relationship between MTCI and canopy N[area] was stronger compared to canopy N[%].Biogeosciences Discuss., https://doi.org/10.5194/bg-2017-228Manuscript under review for journal Biogeosciences Discussion started: 29 June 2017 c Author(s) 2017.CC BY 3.0 License.4.4 Influence of the plant functional type on the MTCI canopy N relationship The relationships between MTCI and both canopy N[%] and canopy N[area] were influenced by the PFT of the plots.The relationship between MTCI and canopy N[%] was stronger for DBF and mixed plots compared to EBF and ENF plots while the opposite was true for the MTCI-canopy N[area] relationship.In the ongoing discussion about the mechanisms underlying the remote detection of canopy N, some authors argued that the differences in canopy structure between different PFTs were causing the observed relationship between canopy N and remote sensing data, rendering it spurious (Knyazikhin et al., 2013).
for the overall relationships between MTCI and either canopy N[%] or canopy N[area], respectively.We highlighted the differences between PFTs and both canopy N[%] and canopy N[area]: the relationship between MTCI and canopy N[%] was stronger for DBF and mixed plots while canopy N[area] was more linked to MTCI for ENF plots.Such differences in relationships between MTCI and either canopy N[%] or canopy N[area] were already observed in grasslands ecosystem.Our results showed that MTCI could be related to canopy N for some individual PFTs, indicating an influence of the PFTs on the MTCI-canopy N relationship.The methodology developed in this study could be investigated at larger scales in different types of ecosystem.While this could already be undertaken using the ENVISAT MTCI 10 years time series as it is almost globally available, ESA new Sentinel-2 satellite launched on 23 June 2015 yields reflectance data at improved spatial and temporal resolution than ENVISAT-MERIS.Canopy N estimates collected through larger scales applications could be exploited in vegetation modelling studies including both the C and N cycles.6 Biogeosciences Discuss., https://doi.org/10.5194/bg-2017-228Manuscript under review for journal Biogeosciences Discussion started: 29 June 2017 c Author(s) 2017.CC BY 3.0 License.

Figure 1 .
Figure 1.Map showing the forest plots location in the region of Catalonia, north eastern Spain.DBF = Deciduous Broadleaf Forest, EBF = Evergreen Broadleaf Forest, ENF = Evergreen Needleleaf Forest, mixed = mixed forest.

Figure A 1 .
Figure A 1. Allometric relationship between the branch diameter (mm) and the foliar dry mass (g) for Pinus sylvestris (n = 208).This is an example of the relationships used to obtain the foliar biomass data of the forest plots.Modified from (Gracia et al., 2004).

Figure A 2 .
Figure A 2. The upper right part of this figure shows the Pearson correlation matrix between canopy N[%] (%N), canopy N[area] (g m-2) and foliar biomass (g m-2) variables for deciduous broadleaf forest plots (DBF), n = 80.The diagonal presents the histogram of the variable on the x-axis, while the y-axis represents the number of counts.The lower left part of this figure represents the scatterplots between the variables.

Figure A 3 .
Figure A 3. The upper right part of this figure shows the Pearson correlation matrix between canopy N[%] (%N), canopy N[area] (g m -2 ) and foliar biomass (g m -2 ) variables for evergreen broadleaf forest (EBF) plots, n = 186.The diagonal presents the histogram of the variable on the x-axis, while the y-axis represents the number of counts.The lower left part of this figure represents the scatterplots between the variables.

Figure A 4 .
Figure A 4. The upper right part of this figure shows the Pearson correlation matrix between canopy N[%] (%N), canopy N[area] (g m -2 ) and foliar biomass (g m -2 ) variables for evergreen needleleaf forest (ENF) plots, n = 563.The diagonal presents the histogram of the variable on the x-axis, while the y-axis represents the number of counts.The lower left part of this figure represents the scatterplots between the variables.
. In this context, the methodology applied in this article could be a valuable alternative to explore canopy N detection at larger scale.Using published data from an extensive field plot inventory, we were able to relate both canopy N[%] and canopy N[area] to MTCI at different spatial resolutions.As MTCI time series are readily and almost globally available, it could eventually be possible to apply our approach at a broader scale in different types of biomes.The results obtained for DBF species and Fagus sylvatica in particular suggest that this method may be efficient at estimating canopy N in temperate forests.Obtaining reliable canopy N estimates over larger scale and for diverse and globally distributed type of vegetation would be necessary to be able to use such estimates to improve and calibrate global vegetation models.Biogeosciences Discuss., https://doi.org/10.5194/bg-2017-228Manuscript under review for journal Biogeosciences Discussion started: 29 June 2017 c Author(s) 2017.CC BY 3.0 License.

Table 1 . Descriptive analysis of canopy nitrogen (N) concentration (N[%], g 100g -1 ), foliar biomass (g m -2 ) and canopy N content (N[area], g m -2 ) by tree species. PFT = Plant Functional Type, DBF = Deciduous Broadleaf Forest, EBF = Evergreen Broadleaf Forest, ENF = Evergreen Needleleaf Forest, mixed = mixed forest, min = minimum value, max = maximum value, mean = averaged value, sd = standard deviation, a codominant plots refer to the plots where two tree species were dominant in the canopy, b foliar biomass data was lacking for five of the plots. Foliar biomass and canopy N content statistics are thus measured on a restricted number of plots.
Biogeosciences Discuss., https://doi.org/10.5194/bg-2017-228Manuscript under review for journal Biogeosciences Discussion started: 29 June 2017 c Author(s) 2017.CC BY 3.0 License.