Remote Sensing Estimating Canopy Nitrogen Content in a Heterogeneous Grassland with Varying Fire and Grazing Treatments: Konza Prairie, Kansas, Usa

Quantitative, spatially explicit estimates of canopy nitrogen are essential for understanding the structure and function of natural and managed ecosystems. Methods for extracting nitrogen estimates via hyperspectral remote sensing have been an active area of research. Much of this research has been conducted either in the laboratory, or in relatively uniform canopies such as crops. Efforts to assess the feasibility of the use of hyperspectral analysis in heterogeneous canopies with diverse plant species and canopy structures have been less extensive. In this study, we use in situ and aircraft hyperspectral data to assess several empirical methods for extracting canopy nitrogen from a tallgrass prairie with varying fire and grazing treatments. The remote sensing data were collected four times between May and September in 2011, and were then coupled with the field-measured leaf nitrogen levels for empirical modeling of canopy nitrogen content based on first derivatives, continuum-removed reflectance and ratio-based indices in the 562–600 nm range. Results indicated that the best-performing model type varied between in situ and aircraft data in different months. However, models from the pooled samples over the growing season with acceptable accuracy suggested that these methods are robust with respect to canopy heterogeneity across spatial and temporal scales.


Introduction
Research efforts into grassland ecosystems are substantial [1][2][3][4], given that grasslands are the potential vegetation covering approximately 36% of the earth's surface [5], and one of the largest vegetative provinces in North America [6].Along with fire and climate, herbivore grazing plays a central role in altering the structure and function of grassland ecosystems [7][8][9].Grassland heterogeneity is a major outcome of fire and grazing activities, which in turn affects future forage pattern of grazers.This interplay between grassland heterogeneity and grazing pattern is of great interest [10,11] due to its significant influences on grassland processes through nutrient redistribution and cycling.Quantitative estimates of canopy biochemical properties are essential for understanding the forage quality distribution and heterogeneity in grassland ecosystems.Among the more important of these biochemical properties is foliar nitrogen content.Nitrogen (N) is an indispensable element in the composition of amino and nucleic acids in all living organisms and an often limiting nutrient in natural ecosystems including grasslands [12].N cycling is a key biogeochemical flux through the earth system, where vegetation is involved in uptake, storage and exchange with other components of the system.Understanding the concentration and distribution of N within vegetative canopies is therefore important for addressing a wide variety of applied and systemic questions in biospheric science.
There is widespread interest in developing methods to estimate the distribution of foliar quality of vegetation canopies using remote sensing [13,14].Traditional methods for measuring nitrogen in vegetative canopies based on laboratory chemical analysis of field-collected foliar samples are complex and expensive.This is especially true when a large number of samples are needed to adequately characterize N distribution over a relatively large area.To address such questions using traditional methods, samples typically consisting of a few grams of canopy material must be collected from point sites distributed throughout the area of interest.Each individual sample represents only a few square centimeters of the canopy, thus some sort of spatial interpolation/extrapolation scheme must be used to estimate values between the point samples if nitrogen is to be mapped over an extensive area.Interpolation schemes can themselves be a significant source of error.
In general, remote sensing based methods for retrieving canopy nitrogen have used high spectral resolution (hyperspectral) reflectance spectra collected from leaf or canopy surfaces [13,[15][16][17][18]. Remote sensing of canopy nitrogen is complicated by the interaction of light in the VNIR/SWIR (400-2500 nm) spectral region with nitrogen in foliar tissues.Nitrogen in foliar tissues directly affects the reflectance spectrum resulting from absorption features due to molecular nitrogen in these reflectance spectra [19].However, the relevant spectral regions for direct analysis of canopy N are in the SWIR spectral region, and are thus not accessible to spectrometers or imaging system which lack sensitivity in this region.Other commonly used methods for retrieving estimates of N from vegetation canopies depend on indirect estimates that exploit the fact that foliar nitrogen is concentrated in plant chloroplasts [20,21].Chlorophyll absorption bands are very prominent in the visible region of the spectrum, and are readily accessible to current remote sensing detectors.The efficacy of using the indirect methods in remote sensing studies for estimating nitrogen concentration via the use of chlorophyll absorption features is now well known and widely used [22][23][24][25][26].
A main challenge of retrieving canopy-level N using remote sensing is that the structure and composition of vegetation canopies is complex, affecting the reflectance of the canopy and therefore the signal received by the sensor.Even relatively homogeneous canopies such as crops are structurally complex mixtures of living and senescent plant tissues, stems and other components.These diverse materials, which in sparse canopies may also include the soil background, are integrated within the field of view of the remote sensor, producing a mixture of spectral elements in the received signal.In heterogeneous canopies (e.g., grasslands) with a diverse mixture of species and growth forms, the signal received at-sensor can be an even more complex mixture spectral responses from diverse reflective materials [13,27].
Available methods for normalizing spectral signals to minimize the effects of canopy heterogeneity have mainly concentrated on removing background effects in order to -isolate‖ the biochemical signal of interest.These methods include: (1) use of the first derivative of the spectral reflectance curve and (2) continuum removal [28][29][30].The first derivative of reflectance characterizes the changes of the original reflectance by which the background effects are countervailed [31].Continuum removal provides a common baseline that allows the spectral features to be standardized and compared [21,32].Additionally, hyperspectral vegetation indices (VIs) integrating spectral values at different wavelengths can normalize spectral features [33][34][35] and are thus potentially applicable to the canopy heterogeneity problem.The standardized difference between two spectral values is the most commonly used form of VI, due to its inherent scaling of the spectral signals and reduced sensitivity to background variations [27,36].
The results we present here are part of a larger project aimed at understanding the relationship between canopy nutrient quality and grazer behavior in tallgrass prairie.The interaction between distribution of foliar nutrients (e.g., N) in a grassland canopy and the forage pattern of ungulate grazers is of special interest since these grazing activities alter the structure and composition of vegetation and therefore influence vegetation quality [37,38].Grazed sites often yield higher nutritional value with greater plant species diversity than do ungrazed sites [39][40][41] due to the capability of plants to overcompensate with regrowth for low levels of herbivory [42].Vegetation quality is one of the major factors that influence grazing strategies of herbivores, and the sites with higher nutrient concentration are likely to be more frequently used by grazers [1,3].Canopy N is a key index of grassland forage quality, and detailed information on its spatial and temporal distribution derived from remote sensing is essential to the goals of this research.
The problem of remote sensing nitrogen in tallgrass prairie canopies is conceptually similar to other canopies, although in practice it is complicated by the heterogeneity problem discussed above.Among grassland canopies, tallgrass prairie is especially heterogeneous, consisting of a diverse, intermingled mix of herbaceous and woody species.Several interacting biophysical factors, including fire, topography, and grazing determine this complex spatial distribution of canopy properties [38,43].Tallgrass prairie canopies can have more than twenty species present in a single square meter, each with different growth forms and leaf structure, thus producing a closely intermingled mixture of species.In addition, these intermingled species often have different photosynthetic pathways.Most of the graminoids (grasses and sedges) are C4 plants, they photosynthesize via a pathway in which the first intermediate product is a 4-carbon molecule.The majority of forbs (flowering species) are C3 photosynthesizers, whose first intermediate product is a 3-carbon molecule.In general, C4 plants tolerate temperature extremes better than C3s [44].All these factors influence the spectral response of the canopy, making the aggregate reflectance signal from even a small area of the canopy a complicated mixture of individual components, and potentially making the retrieval of biophysical information from tallgrass prairie much more difficult.Our objective is to develop and evaluate empirical models using hyperspectral remote sensing data that can be used for estimating and mapping canopy nitrogen in tallgrass prairie canopies.

Study Site
This research was conducted at the Konza Prairie Biological Station (KPBS), near Manhattan, KS, USA (39°05′N, 96°35′W, Figure 1).KPBS is located in the Flint Hills, the largest extant area of tallgrass prairie in North America.The site covers an area of 34.87 km 2 , and is characterized by a continental climate with warm, wet springs, hot summers and dry, cold winters.Approximately 75% of the annual precipitation (~826 mm) is concentrated during the period from April to September, supporting vegetation growth in the prairie.The vegetation in KPBS is characterized by C4 perennial grasses interspersed with C3 forbs and a few C3 grasses, resulting in great species diversity.Grasses compose 80%-90% of total vegetation cover; dominant grass species include Andropogon gerardii, Sorghastrum nutans, Bouteloua gracilis, Panicum virgatum, and Schizachyrium scoparium.Forbs constitute a minor component of the canopy but account for much of the species diversity.Common forbs include Aster ericoides, Psoralea tenuiflora, Solidago missouriensis, Soldiago rigida, Liaris aspera, Vernonia baldwinii and Ambrosia psilostachya [38,45].A long term watershed-level experiment at KPBS is designed to investigate the effects of fire and grazing on the structure and function of the tallgrass prairie vegetative community.The site is divided into more than 50 watersheds, each with varying combinations of fire return frequencies (1,2,4,10, and 20 years) and one of three grazing treatments (grazed by American bison (Bison bison), grazed by domesticated cattle (Bos taurus), and ungrazed).Since this research was conducted as part of a larger study of grazing behavior in Bison bison, we selected three watersheds (with a total area of ~3.40 km 2 ) in the bison grazing area for our analysis (see Figure 1 for locations of selected watersheds).In situ hyperspectral data and foliar samples were collected in three watersheds (N1B, N4D and N20B) with fire frequencies of one year, four years and 20 years, respectively.Variations in fire frequency influence the proportion of green vegetation, senescent materials, stems, and litter in the canopy, which in turn influences spatial patterns of grazing by bison [46][47][48].All of these fire-related factors contribute to a complex, spatially variable canopy.Thus, our selected combination of fire and grazing treatments yields a study area with widely varying species composition and canopy density, resulting in a very heterogeneous canopy.

Data Collection
Data were collected on four dates spanning the 2011 growing season.Each data collection effort consisted of a field component (Table 1), coordinated with an aircraft overflight.Field data were collected from fifteen 10 m × 10 m plots randomly distributed in each of the three designated watersheds, giving a total of 45 sampling sites.Random selection of sites ensured that the samples we collected were representative of species composition over the entire watershed.Technical problems limited the usable data from the September field survey to only 27 sites.In each plot, subsamples consisting of six randomly selected points were distributed around the center of the plot (Figure 2).Two types of nondestructive data sampling were performed at each site, (1) spectral reflectance; and (2) PAR (above and below canopy).In addition, a destructive sample of plant tissue was collected for laboratory analysis of nitrogen content.Hyperspectral reflectance data were taken for each subplot using an Analytical Spectral Devices (ASD) FieldSpec Pro portable spectrometer (Analytical Spectral Devices, Boulder, CO, USA).Spectra were collected from a height of 1.5 m above the target using a 25° field of view optic, yielding a ~0.4 m 2 circular field of view, with a diameter of ~0.7 m on the target.The acquired reflected spectra range from 350 to 2500 nm, with spectral resolution ranging from ~3 nm at the 700 nm wavelength to 10-12 nm for wavelengths of 1050-2500 nm.Individual spectra were converted to reflectance by normalization to a Spectralon reflectance standard.For each plot, a single spectrum was calculated by averaging the reflectance values of each wavelength from the spectra collected from each subplot.On two of the four sampling dates (29 June-1 July and 1-2 August 2011), PAR samples were collected from the same IFOV as the spectral data using an Accupar LP-80 line ceptometer.Geographic coordinates for the center of each plot were collected using a portable GPS receiver.Immediately following the non-destructive sampling, five foliar samples from each of the six subplots, representative of the plant species (grasses and forbs) in each subplot, were clipped from the same IFOV as the spectral and PAR data.In the laboratory, these samples were mixed and ground after drying at 60 °C for 48 h.Dried ground samples from three randomly selected subplots of the six, sufficient to characterize the N status in the whole plot, were then used for foliar nitrogen concentration determination by laboratory elemental analysis using a Costech ECS 4010 (Costech Analytical Technologies Inc., Valencia, CA, USA).Foliar nitrogen content for each plot was determined by averaging the values from the three samples.
Hyperspectral imagery were collected concurrently with each field sampling session using an AISA Eagle camera mounted on a Piper Warrior aircraft operated by the Center for Advanced Land Management Information Technology (CALMIT) of the University of Nebraska-Lincoln.On each date, 4-5 lines were flown, covering the whole area of KPBS and including the three watersheds used in this analysis.The aircraft was flown at an altitude yielding a spatial resolution of 2 m × 2 m.Spectral resolution varied between ~3 nm for the May-August flights, and ~10 nm for the September flight.The camera was configured with a spectral range from 435 to 950 nm (note that the ASIA Eagle is not sensitive beyond 970 nm).

Canopy Nitrogen Calculation
Since our interest was in estimating canopy nitrogen content, converting leaf nitrogen (Nleaf) samples to canopy scale was a key step in data analysis.Nitrogen content can be scaled from the leaf to the canopy level using the canopy-integrated method [49], which relates canopy nitrogen (Ncan) to Nleaf with respect to a leaf area index (LAI): LAI can be estimated from above and below canopy PAR measurements using: where θ is the solar zenith angle.As earlier noted, direct ceptometric PAR measurements were not available for all of the data collection periods.In addition, PAR-method was limited in senescent canopies, which usually overestimated the green LAI.We therefore used an alternative method of calculating LAI, based on the Normalized Difference Vegetation Index (NDVI): 3.706 0.08 where NDVI was calculated from the ASD reflectance values by integrating the wavelengths of 780-800 nm for broadband red reflectance, and 670-678 nm for the broadband NIR.This empirical relationship was validated using data of the two dates for which both ceptometry and spectral reflectance were available, and was found to be reasonably accurate (R* = 0.85, based on jackknife analysis with sample size = 27) [50].

Spectral Reflectance Indices Calculation-Field-Collected Spectra
Preparation of the field measured spectral reflectance data for analysis consisted of numerical differentiation of the spectral curves, continuum removal, and calculation of relevant vegetation indices.Both the first derivative calculation and continuum removal are sensitive to random noise in the spectral signal, so the raw reflectance spectra were smoothed prior to further analysis using a three-point moving average filter.Because the actual spectral resolution of the spectrometer is ~3 nm, we resampled the field spectra by averaging reflectance at every three nanometer increments into one value, similar to the way in which the aircraft sensor sampled spectra.The first derivative of reflectance, a transformation which has been shown to be useful for extracting biophysical information from spectral curves [28,31], was calculated in the spectral region of 562-600 nm.This spectral region has a subtle absorption feature (Figure 3a) sensitive to leaf chlorophyll concentration [51,52], and leaf chlorophyll concentration has been demonstrated to be linked with foliar nitrogen more or less [20,21,53,54].As earlier noted, we did not use the nitrogen absorption feature centered at 1510 nm because the spectral range of the aircraft-mounted sensor (435-950 nm) does not include the SWIR spectral region.The ASD spectrometer does include this region, but we elected not to use it in order to keep the methods applied to both the field and aircraft data sets as parallel as possible.Derivatives were calculated using: where Di is the first derivative at the wavelength i; SRi is the smoothed reflectance value at the wavelength i, and λi is the value of the wavelength i. Continuum removal standardizes the canopy spectral signals using a linear hull (continuum line) to fit over the top of a given absorption feature in the spectrum [29,30].The continuum removed spectral value at a wavelength is calculated as the ratio of the value in the reflectance spectrum to the corresponding value in the continuum removal line at the given wavelength (Figure 3b).We used the same absorption feature for continuum removal (562-600 nm) that we used for the first derivative calculation, for the same reasons.The continuum removed value at each wavelength was calculated as: where CRi is the continuum removed spectral value at the wavelength i; SRi is the smoothed reflectance at the wavelength i; and SLi is the value at the wavelength i in the continuum line fitted over the absorption spectral feature.The standardized difference vegetation index (SDVI) was developed for various derivative and continuum removed spectral features using: where Si is the spectral feature at the wavelength i. SDVI derives from the theoretical underpinning of ratio-based vegetation indices such as NDVI [55] that scales the difference in two spectral features which behavior diversely.

Spectral Reflectance and Indices Calculation-Aircraft-Collected Spectra
Aerial imagery was processed using the ENvironment for Visualizing Images (ENVI) software (Exelis Visual Information Solutions, Inc., Boulder, CO, USA).Prior to analysis, the FLAASH model implemented within ENVI was used for atmospheric correction.For each of the four FLAASH runs, the atmospheric model was set to midlatitude summer and the aerosol model to rural.Visibility for each day was estimated using meteorological optical range (MOR) data from the Topeka, KS, National Weather Service station (located approximately 80 km from the study site).Visibilities were 16 km for the May, July and August acquisitions, and 14 km for the September acquisition.Since our analysis was restricted only to grassland canopies, supervised Maximum Likelihood classification was used to classify land cover on the four images captured during May-September into four classes, including gallery forests (i.e., wooded areas along drainages), roads, sparse grasses and dense grasses.The gallery forest areas and the roads were masked from the imagery and not included in subsequent analysis.Since the gallery forests were expanding during the growing seasons, the forest mask varied over time.A union of masks for the four data collection sessions was built and applied to make sure that the mask looked the same and all forested pixels were removed.
Plots for in situ data collection were located on aerial imagery using the plot center coordinates collected in the field.Spectral data from airborne imagery for each plot were obtained by calculating the mean of spectral reflectance data within a window of 5 pixel × 5 pixel (10 m × 10 m) centered around the plot center on the imagery.These averaged values were then subjected to the same analysis (smoothing, continuum removal, numerical differentiation, index calculation) as the field data.The spectral region used was 562-600 nm, as in the field data analysis.

Field Data-Descriptive Statistics
Descriptive statistics of in situ Nleaf, LAI and N can are summarized in Table 2. Ncan was determined by the product of measured Nleaf and LAI, where LAI was derived from NDVI based upon a validated empirical exponential relationship between LAI and NDVI (Formula (3) in Section 3.2.1).The use of empirical LAI-NDVI relationship should be prudent due to NDVI saturation at high LAI values (LAI > 2.5), but possibly appropriate in this study given the low ranges of LAI overall (all LAI values < 2.5).Relatively high ranges of Nleaf appeared in May, and of LAI appeared in July.Average N can was higher in May and July than in other months.Various fire intensity influenced nitrogen spatial distribution in different watersheds.Watershed N1B with high fire frequency of one year had higher levels of N can than the other two watersheds N4D and N20B throughout the growing seasons from May-September.Analysis of correlations (Pearson's r) between N can and the various spectral indices was used to evaluate the strength of relationships between canopy N and spectral reflectance across the 562-600 nm absorption feature for each watershed individually, and for the three watersheds together.For each pooled dataset of three watersheds, thirty samples were randomly selected from the total samples of 45 in the May-August data (18 from the total 27 in September data) for correlation analysis.For the first derivative data (Figure 4), significant wavelengths related to nitrogen content included the regions from 563-566 nm and 593-596 nm, where the absolute r-values were greater than 0.75 across all watersheds and all dates.A common region with relatively low absolute r-values was at wavelengths of 569-575 nm centered at 572 nm, corresponding to the region with maximum absolute values of the first derivative data, where the original reflectance values decreased fastest.Absolute r-values for the pooled datasets (May-September) were greater than 0.86 (Figure 4e), which indicates a generally strong correlation between foliar N and the first derivative data in the given absorption spectral feature across various spatial and temporal scales.
For the continuum removed spectra (Figure 5), absolute r-values at wavelengths of 569-593 nm in May-August had a range from 0.83 to 0.89 (Figure 5a-c).Spectral region of 578-593 nm in September with the absolute r of 0.72-0.76(Figure 5d), was strongly related to N can across the three watersheds.Correlations between the nitrogen and the spectral dataset incorporating data for the three watersheds throughout the entire growing season were relatively strong in the region of 569-584 nm with the absolute values of r greater than 0.90 (Figure 5e).
Based on the results of the correlation analyses, absolute values of r in the region of 562-569 nm start with low ranges and vary dramatically.This can be accounted for the continuum removal line that presses close to the reflectance spectral curve at the left end of the absorption feature, resulting in highly -compressed‖ continuum removed values in this region (Figure 3b).It was better to exclude this region in analysis as the highly -compressed‖ continuum removed values reveal little useful spectral variation information.We therefore calculated the SDVI using the continuum removed spectral features at the wavelengths of 572 nm and 581 nm-CR572 and CR581.The wavelength at 572 nm is in the starting portion of the -normal‖ continuum removed spectral region; 581 nm is around the middle point between 562 and 600 nm, which is the nadir in the continuum removed absorption trough.Further correlation analysis supports the use of SVDI formed at these wavelengths as an indicator of canopy N, with absolute r-values between N can and the SDVI varying from 0.78 to 0.96 (Table 3).

Canopy Nitrogen Modeling From in situ Data
Thirty samples were randomly selected from the May-August data (18 from the September data) for model development.The remaining 15 samples (nine from September) were reserved for validation.Descriptive statistics of the in situ datasets of canopy nitrogen contents for modeling and validation (Table 4) demonstrate that the ranges of the modeling datasets were in general slightly larger than that of the corresponding validation datasets, indicating that the modeling datasets cover the full range of N contents well.Similar values of the mean and standard deviation (SD) between the modeling and validation datasets suggest similar central tendencies and dispersions between the two types of datasets.This verifies a reasonable selection of modeling and validation values.Monthly models of canopy nitrogen retrieval using hyperspectral data for May, July, August and September, and a general model for the whole growing season were derived using regression analysis.All of possible wavelength combinations were examined using the first derivative and continuum removed spectra, subject to restrictions imposed by the modeling method.Chief among these choices was to limit the number of independent variables.Although a model may fit the data better when incorporating more independent variables, an excessive number may result in an overfit model that will not perform well in validation.For our data, the need to limit model inputs is even more pronounced since the independent variables are spectral features concentrated in a narrow region, which might be strongly correlated with each other and thus might not be all statistically significant in the model.We therefore used combinations of no more than three independent variables for each model.In addition to the multivariate models based on individual spectral bands, other models were derived using only the SDVI.All of the regression models were validated by comparing the measured nitrogen from the samples reserved for validation with the values estimated using the remote sensing methods.
The best-performing models from the regression analysis are summarized in Table 5.Note that while the form of the regression models remains the same, the empirical coefficients of each model changed in each month of data collection.This is hardly surprising, given the empirical nature of the modeling and the continual changes in the canopy over the growing seasons.The coefficient of determination (R 2 ) was generally lower for the validation data than for the model-development data.Again, this is not unexpected.The similarity of the error metrics from the model development step and the validation step suggests that the models have little overfitting and should therefore be robust when applied to new data.In general, there was little difference between the two best multivariable band-based models and the univariate SDVI model.Model performances were also similar throughout most of the growing season, although there was a notable decrease in model performance in September.This result could be due to differences in canopy phenology.At the time of the September data collection the canopy was beginning to senesce, so the fraction of green material in the canopy was likely less than in the previous data collection periods.It could also be an artifact of the smaller sample size, which resulted in a validation set with samples size = 9.
The results from pooling all the data (May-September) were similar to, and in some cases better than, the results from the individual months.This is a particularly important outcome, since it suggests that one model may be sufficient to predict Ncan across the entire growing season, eliminating the need to use multiple models depending on the growth state of the canopy.The larger sample size available for both model development and validation (the result of combining all of the measurements) might also explain some of the better validation results in the pooled sample.

Canopy Nitrogen Retrieval from Airborne Imagery
Results from the in situ spectrometry (Sections 4.1 and 4.2) suggest that SDVI calculated using continuum removed reflectance in the spectral region of 562-600 nm yielded the best model for retrieving Ncan.While the validation metrics for the SDVI model differed little from the multivariate models, the intrinsic parsimony of the SDVI model (because it is based on a single independent variable) strongly support its further use in Ncan retrieval models.However, we chose to analyze the data from the aerial imagery as we did the in situ data.Here, we tested all forms of retrieval model in order to confirm if SDVI still yielded the best prediction models when applied to a data set with different spatial and spectral resolution.We therefore repeated the analysis described in Section 4.2, using the atmospherically corrected reflectance data collected during the four aircraft overflights.
Results of modeling and validation for nitrogen retrieval from airborne imagery (Table 6) showed that reasonably accurate Ncan estimation models can be determined from datasets of either the first derivative and continuum removed reflectance or the SDVI values.Feasible models for September suggested that the wider spectral bands (10 nm vs. 3 nm, see Section 3.1) were still adequate for nitrogen estimates.Despite variations in selection of wavelength combinations for optimal nitrogen modeling in different months, general nitrogen retrieval models could be developed for the growth seasons during May-August.
The models developed from the first derivative and the continuum removed spectral data showed higher R 2 values for regression modeling and lower ranges of RMSE for model validation than the SDVI model, indicating that multivariate models performed better for the canopy N estimates based on the aircraft data.As with the in situ data, there was little difference between the accuracy of each month's model and the overall accuracy of the pooled (all-months) models.Again, this is a significant result, as it shows that Ncan can be extracted across the growing season using the same retrieval model, simplifying the processing of data collected at various times during the growing seasons.

Comparison of Results and Nitrogen Mapping using Aerial Imagery
Even though the in situ data were resampled with similar band passes to the aerial imagery, differences in other physical mechanisms and measurement conditions between the two sensors possibly make it still difficult for the two sets of measurements to be interchangeable.Wavelength combinations and empirical model coefficients vary between optimal in situ models and the corresponding aerial imagery models.This is not unexpected given the intrinsic features of empirical methods.However, similar results of model performance and validation suggest that these methods are feasible and stable when applied to respective datasets of field spectra and aerial imagery.In other words, the field model is applicable to the field data, and the aircraft model to the aircraft data, but the two cannot be interchanged.
Maps of Ncan for the entire study site were made using the best-performing equation among the pooled retrieval models.Because the pooled model developed from the first derivative of spectral reflectance at 583 and 592 nm had the lowest RMSE, we used it to create maps of Ncan for May-August (Figure 6a-c).The Ncan map for September was created using the first derivative at 582 nm due to the different spectral resolution of aerial imagery in September (Figure 6d).These Ncan maps support the statistical comparison across various spatial and temporal scales (Table 2 in Section 4.1), and also provide more spatial detail and phenological characteristics.

Assessment of N Relation to Fire, Vegetation Density and Topography
The patterns in the Ncan maps were expected to be affected by fire frequency, vegetation density and topography over the study site.To evaluate these influences on Ncan distribution, Ncan estimations on each map were stratified into the low, medium and high levels by 3-quantiles, and then compared to the land cover type (sparse grasses vs. dense grasses, classified by supervised Maximum Likelihood in Section 3.2.3)and a digital elevation model (DEM, with spatial resolution of 2 m × 2 m) [56].The effect of fire frequency can be seen in the Ncan stratifications between watersheds (Figure 7).Ncan values in watershed N1B with fire frequency of one year are generally higher than in the other two watersheds.This is especially noticeable in the May and September imagery, where the effects of canopy and treatment heterogeneity are most prominent.The July image shows relatively even distribution of Ncan, consistent with a period in which the canopy is more uniformly green and actively photosynthesizing.The spatial relationships between Ncan, canopy density and topography observed in the three watersheds of our study (Figures 7 and 8) are fairly straightforward, similar to those observed throughout KPBS [57].Values of Ncan coincide directly with canopy density.Canopy density, in turn, is closely related to topography through its effect on soil thickness and moisture.The densest canopies are found in flatter, low-lying areas, where the soil depth is greatest and soil moisture tends to be highest.Thinnest canopies occur on hill sides, where slope is greatest and the thin, rocky soil layer is unable to retain the moisture needed to support a denser canopy.Intermediately thick canopies occur mainly on hilltops, drainage divides, and other topographic highs, where slopes are low and soils depths are greater than on slope sides but less than in the bottom areas.This topographic effect on Ncan distribution is most noticeable in May, but becomes less so as the canopies become more developed during the growing seasons.

Conclusions
Our results show that canopy nitrogen can be estimated using empirical models based on hyperspectral remote sensing data collected from a tallgrass prairie canopy subjected to varying fire treatments and characterized by a heterogeneous mix of species and canopy structures.Spectra collected from both ground-based and aircraft-mounted spectroradiometers were used to develop empirical retrieval models.A simple model based on the SDVI was best for estimating Ncan from the ground-based data, whereas a multivariate regression model performed best using the aircraft data.However, the differences between the models were small, and in general models from both data types were comparably accurate.
Comparing our results with similar studies recorded in the literature suggests that the performance metrics for our best empirical models were generally comparable to those of similar analyses conducted over other canopy types.In past efforts, remote estimation of canopy N content yielded the best results when applied to agricultural ecosystems of various types [35,[58][59][60].This is perhaps to be expected.
Although crop canopies are not simple monocultures they are more homogeneous than most natural vegetation canopies and therefore present a more uniform surface from which to extract biophysical information.Forest canopies, a natural surface but one typically characterized by a less diverse range of species actually visible to the sensor, also tended to provide slightly better results in remote analysis of nitrogen when compared to ours [61][62][63].Among those studies done in grasslands or other more heterogeneous environments, our results were quite similar to those reported by Boegh et al. [64]-a study in which the problem of canopy heterogeneity was considered directly-and to those reported by Asner et al. [13], Kooistra et al. [65], and Mitchell et al. [66].
While our results agree with and support previous work on retrieving canopy biochemical concentrations via remote sensing, there are also findings of further relevance.Most importantly, our results show that empirical methods for N retrieval can be robust when applied to a canopy with both spatial and temporal heterogeneity.In this case, the spatial heterogeneity was introduced as part of the research protocol by including three experimental units (i.e., watersheds) whose differing burn treatments yield very different canopy conditions.In addition, all three of these watersheds studied here were subjected to grazing by large ungulates, further contributing to a diverse and complex canopy.For each month of data collection, the same empirical model was able to estimate Ncan in each watershed with comparable accuracy.This is a significant result since it shows that Ncan (and possibly other elements of stoichiometric interests) can be mapped across diverse canopies without the need for stratification and the use of multiple models.In addition, our results show that empirical models can perform well across the entire growing season, again obviating the need for multiple empirical models calibrated to account for phenology.

Figure 1 .
Figure 1.Study sites at Konza Prairie Biological Station (KPBS).KPBS include more than 50 watersheds with different fire frequencies.Samples for this study were collected in three watersheds, N1B, N4D and N20B, with fire frequencies of one year, four years and 20 years respectively.The total area of the three watersheds is ~3.40 km 2 .

Figure 2 .
Figure 2. Plot design.Fifteen 10 m × 10 m plots with six subplots spirally arranged around the plot center were set up for data collection in each of the three watersheds.Hyperspectral measurements for each subplot were taken from a 0.4 m 2 circular field of view with a diameter of 0.7 m.Five foliar samples were clipped in each field of view for chemical analysis of leaf nitrogen concentration.

Figure 3 .
Figure 3. (a) The absorption feature related to leaf chlorophyll concentration in the spectral region of 562-600 nm.(b) Continuum removal applied to the absorption feature.

Figure 4 .
Figure 4. Correlations between canopy nitrogen and the first derivative spectral data at the wavelengths of 562-600 nm for three watersheds N1B, N4D and N20B respectively and as a whole in (a) May, (b) July, (c) August, (d) September, and for (e) four months respectively and the whole growing seasons in May-September across various watersheds.

Figure 5 .
Figure 5. Correlations between canopy nitrogen and the continuum removed spectra at the wavelengths of 562-600 nm for three watersheds N1B, N4D and N20B respectively and as a whole in (a) May, (b) July, (c) August, (d) September, and for (e) four months respectively and the whole growing seasons in May-September across various watersheds.

Figure 6 .
Figure 6.Maps of Ncan distribution in selected watersheds at KPBS in (a) May, (b) July, (c) August, and (d) September.White areas are where gallery forest pixels were masked prior to analysis.

Figure 7 .
Figure 7.Comparison between Ncan stratification and grass density in selected watersheds at KPBS in (a) May, (b) July, (c) August, and (d) September.Overlaps of high Ncan level and dense grasses are in the dark shade; low Ncan level and sparse grasses are in the light.

Figure 8 .
Figure 8.Comparison between Ncan stratification and (a) DEM in selected watersheds at KPBS in (b) May, (c) July, (d) August, and (e) September.Hillslopes are highlighted in the light shade, and flat lowlands are in the dark shade.Coincidence between Ncan distribution patterns and topographic features is most noticeable in May.

Table 1 .
Datasets of in situ measurements.

Table 2 .
Descriptive statistics of in situ Nleaf, LAI and Ncan.

Table 3 .
Correlations between canopy nitrogen N can and the SDVI for the field-collected data.

Table 4 .
Descriptive statistics of the in situ N can datasets for modeling and validation.

Table 5 .
Regression models for canopy nitrogen recovery using field-measured derivative (D), continuum removed reflectance (CR) spectra and SDVI.

Table 6 .
Regression models for canopy nitrogen retrieval using derivative (D), continuum removed reflectance (CR) spectra and SDVI derived from aircraft image.