Integrated Analysis of Productivity and Biodiversity in a Southern Alberta Prairie

: Grasslands play important roles in ecosystem production and support a large farming and grazing industry. An accurate and efﬁcient way is needed to estimate grassland health and production for monitoring and adjusting management to get sustainable products and other ecosystem services. Previous studies of grasslands have shown varying relationships between productivity and biodiversity, with most showing either a positive or a hump-shaped relationship where productivity peaks at intermediate diversity. In this study, we used airborne imaging spectrometry combined with ground sampling and eddy covariance measurements to estimate the spatial pattern of production and biodiversity for two sites of contrasting productivity in a southern Alberta prairie ecosystem. Resulting patterns revealed that more diverse sites generally had greater productivity, supporting the hypothesis of a positive relationship between production and biodiversity for this site. We showed that the addition of evenness to richness (using the Shannon Index of dominant species instead of the number of dominant species alone) improved the correlation with optical diversity, an optically derived metric of biodiversity based on the coefﬁcient of variation in spectral reﬂectance across space. Similarly, the Shannon Index was better correlated with productivity (estimated via NDVI (Normalized Difference Vegetation Index)) than the number of dominant species alone. Optical diversity provided a potent proxy for other more traditional biodiversity metrics (richness and Shannon index). Coupling ﬁeld measurements and imaging spectrometry provides a method for assessing grassland productivity and at a larger scale than can be sampled from the ground, and allows the integrated analysis of the productivity–biodiversity relationship over large areas.


Introduction
Grasslands occur on all continents except Antarctica and occupy nearly 25% of the land surface of the Earth [1].In North America, the Great Plains is the broad expanse of flat land located east of the Rocky Mountains and west of the Mississippi River in the United States and Canada.This land is mainly covered in steppe and grassland, and is often referred to as "prairie."In the United States, prairie area is 1.62 ˆ10 6 km 2 , or 21 percent of the total area of the country [2].In Canada, prairie area is 5 ˆ10 5 km 2 , 5 percent of the total area of the country [3].In Alberta, approximately 7% of the area is covered by dry mixed grass prairie [4], most of which is found in southern Alberta and is contiguous with the Great Plains.
In North America, prairie is an important biome for agriculture and grazing (rangeland).Prairies typically are less productive than croplands but play important roles in carbon sequestration and food supplement for animals.Healthy prairies also sustain soil quality, maintain biodiversity, and provide clean water [2].Among all the ecosystem services provided by prairies, biomass yield is typically the factor that most interests humans due to the direct relationship with the prairie's major products, food and fiber [5].Prairie management typically focuses on maintaining or improving the output of consumable products along with other ecosystem services [4].
Prairie productivity is often estimated through biomass harvesting, which is expensive and time consuming [6][7][8].Traditional methods of estimating prairie biomass can be subjective because they depend on the experience and skill of the field staff [7].The estimation is also affected by the sample size and sampling method [9].Accurate and fast methods of prairie production estimation that quickly acquire information over large areas are needed for range managers to evaluate the condition of prairies and adjust their management regimes [7].
The classic Light Use Efficiency (LUE) model, which originated with Monteith's work [10,11], is often used to address the spatial and temporal dynamics of production from remote sensing [12].The LUE model considers that primary production largely depends on the absorbed photosynthetically active radiation (APAR) captured by plants and the efficiency with which the absorbed energy is converted to fixed carbon.APAR is affected by the solar irradiance, leaf pigment concentration and canopy structure (e.g., leaf area index and leaf angle distribution).The efficiency term is affected by the light energy distribution within the leaf, and is often influenced by pigment composition and activity [13].The LUE model parameterized by remote sensing measurements generally works well for grasslands, where the APAR term dominates and the efficiency term can often be treated as a constant, particularly for short periods [14,15].While some challenges remain in scaling, this provides a simple route to determining grassland productivity for comparison with diversity.
There is a long history of using remote sensing to estimate vegetation properties from satellite or airborne platforms [16][17][18].A number of studies have used remote sensing to estimate the percent cover, production and biophysical properties of grasslands [6,7,19].Many of these studies have shown that the Normalized Difference Vegetation Index (NDVI) is highly correlated with green biomass, leaf area index, and radiation absorption by green vegetation (APAR green ) in grasslands [8,19], consistent with the simple LUE model approach, and allowing ready production estimates over large grassland areas.However, a major limitation in the application of satellite remote sensing in prairie management is the mismatch between the information that range managers need and the spatial scale of most satellite remote sensing.Due to the coarse pixel size (often several tens of meters to kilometers) [7] and broad bands of most satellite images, satellite remote sensing can best provide land cover or production estimates over large landscapes, or over regional or continental scales [2], but satellite data often have trouble providing detailed species distribution and prairie health information at the scale of small management units.
The high spatial and spectral resolution of airborne imaging spectrometry can reflect fine-scale features by detecting physical and biochemical properties of plants in specific narrow bands at small (sub-meter) image pixel sizes.Consequently, airborne imaging spectrometry provides an effective tool to map prairie patterns related to photosynthetic function [14] and species distribution patterns in grassland ecosystems [20].One way to derive detailed spatial and temporal patterns of productivity from such imagery is to calibrate the image data against ground measurements using eddy covariance and field spectral measurements [14].Airborne remote sensing, with its high spatial resolution, provides an ideal tool to address these issues of scale when addressing productivity and diversity.
There is now a large body of literature indicating that biodiversity affects ecosystem production [21] and stability [22,23].Particularly in the face of disturbance, the relationship between diversity and productivity, typically measured as biomass, has been controversial for grassland ecosystems.Some studies report a positive relationship between biomass and biodiversity, whereas others indicate a "hump-shaped" relationship, with biodiversity peaking at intermediate biomass levels [24].This relationship can be affected by resource levels (e.g., fertilizer or irrigation levels) and the degree and nature of disturbance.High productivity and low diversity sites are often highly managed via irrigation or fertilizer application [24].For rangelands, grazing levels can confound this relationship, and highest species richness often occurs at light to moderate grazing levels, which is consistent with the "intermediate disturbance hypothesis" [25,26].Heavy to very heavy grazing can cause a decline in species richness, often due to the invasion of exotic species that lead to a decrease in species diversity while maintaining high levels of productivity [4].
Airborne remote sensing provides an efficient and inexpensive way to assess biodiversity.Traditionally, three basic methods have been used to estimate biodiversity with remote sensing: (1) mapping habitat for key species; (2) mapping species distribution [27][28][29][30] or community composition [31]; and (3) assessing species richness [32][33][34], α-diversity [33] or β-diversity [35] through spatial variation in vegetation optical properties (optical diversity in space), sometimes referred to as "spectral heterogeneity" [36].Using this latter approach, which we call "optical diversity," a variety of methods have been used to capture optical variation as a way to estimate biodiversity.For example, Féret and Asner [37] defined "spectral species" by applying a clustering model to high spatial resolution airborne imagery to map α-diversity and β-diversity in tropical forests.
The spectral species approach assumes that there are unique, definable spectral types ("species") that can be distinguished in image processing, allowing an estimate of biodiversity.A similar concept can be applied in a more abstract level of image information content, without resorting to identifying particular spectral types.We propose that the information content of the imagery itself, which can be expressed as optical diversity, relates to the spatial variation in the spectral data [32,34] and provides an indicator of relative diversity.The advantage of this approach is that it provides an objective method that can conceivably be applied to any image or ecosystem, without having to define categorical species or spectral types, which are likely to change seasonally or spatially.
In this study, we used flux data and field optical data to help calibrate airborne imagery and map ecosystem productivity in a grazed prairie ecosystem in southern Alberta, Canada.Airborne NDVI measurements were calibrated against CO 2 flux measurements and above-ground biomass to estimate landscape productivity, with particular attention to potential errors associated with spatial scale due to various flux footprint assumptions.We also sampled relative levels of biodiversity using optical diversity and field measurements, and by combining airborne data with a vegetation distribution map.Combining three metrics of biodiversity (airborne data, vegetation map and field sampling) to evaluate broad relationships between diversity and productivity across the landscape provided a unique test of the diversity-productivity hypothesis [21] over a large area (10 km 2 ) of this managed prairie ecosystem.

Study Site
The Mattheis Research Ranch is located 150 km east of Calgary, Alberta, Canada (Centre Latitude: 50.90 ˝N, Centre Longitude: 111.88 ˝W) (Figure 1).It covers approximately 5000 hectares, 4/5 of which is native prairie.The vegetation distribution map and code for major vegetation communities (based on dominant vegetation type) in the Mattheis Research Ranch are provided in Supplementary Materials (Figure S1 and Table S1).The landscape included two calibration sites, designated "E3" and "E5" (Figure 1 [38]), where eddy covariance, vegetation biomass, species composition, and field optical measurements (APAR and spectral reflectance) were made.The landscape in this region is covered mainly by dry mixed grass prairie.The study area consisted primarily of grazed prairie, with some nearby wetlands and cultivated areas (largely out of the immediate study area).The grazing regime is a rotational system with relatively short grazing periods of 7-12 days, followed by long recovery periods (6-8 weeks) between grazing periods at moderate stocking.In 2012, grazing commenced in early May-late November, and cattle were outside the direct study area during the overflight.The grazing effects on species composition and structure are not visible from the historical aerial photos, indicating a history of a light grazing regime for this site [39].

Optical Phenology and Flux Measurements
Two flux towers were established within the ranch at two calibration sites labeled "E3" (50.8672 ˝N, 111.9045 ˝W) and "E5" (50.9056 ˝N, 111.8823 ˝W), located approximately 4.5 km apart.A phenology station was set within 10 m of the flux tower at each site.The optical phenology stations measured reflectance of broadband solar radiation and PAR every fifteen minutes.Each optical station had a data logger (H21-001, Onset Computer Corporation, Bourne, MA, USA) and two-band radiometer.One band had two PAR sensors (S-LIA, Onset Computer Corporation, Bourne, MA, USA) and the other band consisted of two PYR (pyranometer) sensors (S-LIB, Onset Computer Corporation, Bourne, MA, USA).Both the upward and downward sensors have cosine foreoptics, providing a relatively large optical footprint (approximately 100 m 2 ) around the flux tower.These data were used to calculate a proxy NDVI [40,41]: NDVI proxy " pρ PYR ´ρPAR q{pρ PYR `ρPAR q (1) where ρ PYR is the reflectance of the PYR band (solar radiation), calculated as the ratio of the reflected solar radiation to the incoming solar radiation across a spectral range from 300 to 1100 nm and ρ PAR is the reflectance of the PAR band that is the ratio between the reflected PAR and the incoming PAR within 400-700 nm.The proxy NDVI values were subsequently used to derive APAR for the LUE model, as described below (Figure 2).Eddy covariance data were collected at both sites using an open-path infrared gas analyzer (IRGA; LI-7500, LI-COR, Lincoln, NE, USA) and a three-dimensional sonic anemometer (CSAT3; Campbell Scientific, Logan, UT, USA), at 2.9 m (E3 site) and 3.0 m (E5 site) above ground.The raw flux data recorded at 10 Hz was aggregated at thirty-minute time intervals as a measure of net ecosystem productivity (NEP = ´Net Ecosystem Exchange, NEE) using the EddyPro (LI-COR, v. 5.1) software package.In our study, we used the sign convention of positive NEP fluxes indicating a net uptake of CO 2 by the biosphere and negative NEP fluxes signifying a net release of CO 2 to the atmosphere, with the reverse convention used for NEE.
Flux and optical data from the E3 site was used to create the APAR-NEE model, instead of data from E5, which had an extensive flux data gap in the growing season.The flux data at E3 site was collected continuously from 16 May to 19 October 2012.At the E3 site, there were three major data gaps in the growing season: May (10 days), June (16 days) and August (8 days).
The peak flux in the summer was used to separate the season into green-up and senescence phases in order to avoid possible effects of hysteresis [15].Considering the date of the airplane flight (17 August 2012), a model of proxy NDVI-NEP in the second half of the growing season was used.

Biomass Harvest
To calibrate productivity estimates, biomass was measured monthly at both calibration sites (E3 and E5) during the growing season in 2012 (11 measurements for each site from May to August, 44 measurements in total).A 30 cm diameter metal sampling ring was placed around the vegetation to be harvested.All above-ground vegetation was cut to a height of <1 cm above the ground.Forb and graminoid components were separately collected and taken back to the lab for further processing.The biomass was sorted, oven-dried, and weighed (g/m 2 ) as separate categories (green and brown).A "green fraction" factor was calculated based on the biomass ratio between green and total (green and brown) dry mass, and this was subsequently used in APAR calculation.

APAR Determination
Fraction of Absorbed Photosynthetically Active Radiation (f APAR) measurements were taken at the same time as the biomass harvest in the field using a light bar (Accupar, Decagon Inc., Pullman, WA, USA), and were calibrated against proxy NDVI measurements (above) to parameterize the f APAR term of the light-use efficiency model.These f APAR measurements were multiplied by the green fraction obtained from biomass harvests and calibrated against the proxy NDVI values to derive a continuous f APAR green .This NDVI-derived f APAR green was subsequently multiplied by the Photosynthetic Photon Flux Density (PPFD) from the phenology stations to estimate radiation absorbed by green canopy material (APAR green ).

Ground NDVI
Spectral reflectance was measured with a dual channel spectrometer (Unispec DC, PP Systems, Amesbury, MA, USA) at biomass and proxy NDVI sampling sites (Figure 2).The detectors collected irradiance and radiance from 350 to 1130 nm with an approximate spectral resolution (FWHM) of 10 nm.The upward-looking channel included a fiber optic and a cosine head to record the solar irradiance.The downward-looking channel included a fiber optic and a field-of-view restrictor that limited the field of view to approximately 15 degrees.In this application, the spatial resolution of each sample on the ground was approximately 0.5 m 2 .
Both upwelling radiance and downwelling irradiance were measured over the vegetation target and a white reference calibration panel (Spectralon, Labsphere, North Sutton, NH, USA), and used to correct for the atmosphere variation and calculate surface reflectance [42].The reflectance (ρ) at wavelength (λ) was calculated as In this equation, L target,λ indicates the radiance measured at each wavelength (λ, in nm) by a downward-pointed detector sampling the surface ("target"), while E target,λ indicates the irradiance measured simultaneously by an upward-looking detector sampling the downwelling radiation.L panel,λ indicates the radiance measured by a downward-pointed detector sampling the calibration panel, and E panel,λ indicates the irradiance measured simultaneously by an upward-pointed detector sampling the downwelling radiation.
Ground NDVI was calculated from spectrometer measurements using Equation (4).
NDVI 680,800 " ρ 800 ´ρ680 ρ 800 `ρ680 where ρ 800 and ρ 680 represent the reflectance at 800 and 680 nm, respectively.Ground spectral reflectance measurements were taken at the same time and over the same 1-ha landscape region as the biomass harvest and proxy NDVI.A linear model between time-series ground NDVI and green biomass was created and later applied to the airborne data for subsequent determination of biomass and APAR green from airborne data (Figure 2).

Airborne Data
An imaging spectrometer (Headwall A Series, Headwall Photonics Inc., Fitchburg, MA, USA) was used to collect airborne data on 17 August 2012.The instrument was mounted on a fixed-wing aircraft (Piper Navajo, Piper Aircraft, Vero Beach Florida) from a height of 1220 m and a speed of 213-222 km/h.The pixel size on the ground was approximately 1.1 m.The imaging spectrometer provided 400~1000 nm hyperspectral images with 3 nm spectral resolution (Full width at half maximum, FWHM).
Several steps were taken to process the raw digital numbers to reflectance spectra.(1) A flat field correction was done to remove patterns inherent in the detector array; (2) A spectral adjustment was taken using the 760 nm oxygen Fraunhofer line as a reference band to correct for wavelength shifts; (3) Three 9 ˆ9 meter calibration targets (white, charcoal, and black) made from polyester fabric (Odyssey, J. Ennis, Edmonton, AB, Canada) were used in the surface reflectance correction.These targets were placed on the ground during the airplane overflight to calculate coefficients for calculating surface reflectance using the empirical line correction [43].Finally, reflectance spectra were degraded to a spectral resolution of 10 nm for better data quality using Equation (5).
FWHM o and ρ λ,FW HMo were the spectral resolution and reflectance of the original spectra, FWHM d and ρ λ,FW HMd were the spectral resolution and reflectance of the target spectra, K λ´λ was a Gaussian kernel function while ş 8 ´8 K λ´λ d λ " 1 [44].Position and rotational attributes (pitch, roll, and yaw) of the airplane during the flight were recorded by an inflight GPS and inertial measurement unit (IMU).This information was used to apply the geo-referencing correction.The images were resampled to 1 meter 2 spatial resolution when applying the geo-referencing correction.Airborne narrow-band NDVI was calculated using Equation (4).
We used a north-south flight line (10 km long ˆ1.3 km wide) that covered both flux tower sites to estimate the productivity in this region (Figure 1).We applied a "space-for-time substitution" strategy by applying an NDVI-green biomass calibration and LUE model (calibrated by data collected from the E3 and E5 calibration sites during the second half of the growing season) to airborne NDVI data (spatially distributed over the landscape) to map the green biomass and NEE distribution in this prairie ecosystem.The airborne NDVI was calibrated using the narrow-band NDVI-green biomass relationship obtained with a ground spectrometer at 11 ground calibration sites located in a 1 ha area at both E3 and E5 sites.The airborne NDVI was calibrated against proxy NDVI to generate the APAR-NEE relationship, which assumed a fixed efficiency (ε) for this prairie ecosystem in the LUE model.Because this NDVI was linearly correlated to both green biomass from harvests (R 2 = 0.82) and Net Ecosystem Exchange (NEE) from flux measurements (R 2 = 0.77) it provided a proxy metric of productivity for this grassland ecosystem (Table 1).

Sensitivity Analysis of Footprint
To help understand the area sampled for CO 2 exchange and its influence on the LUE model, flux footprints were calculated for each valid 30 min period using a parameterization of a Lagrangian stochastic model as described in [45].We ran the footprint model using a fixed boundary layer depth of 1000 m and roughness length and displacement height values derived from canopy height estimates.The cross-wind integrated footprints were summed up into a time-integrated footprint matrix according to wind direction.Aggregation was done in two ways: on a monthly basis at 5 h midday intervals (11:00-16:00) and intervals comprising the entire daytime period (6:00-21:00).
To help design an effective data integration scheme, a sensitivity analysis of the footprint was conducted with the airborne image.We resampled the footprint to the same size (1 ˆ1 m) of the airborne image pixel and calculated NDVI weighted according to the footprint using the following equation.
N was the total pixel number, w i was the decimal percent concentration of each grid cell to the cumulative footprint function values that summed up to 1.
To evaluate the sensitivity of the study results to footprint assumptions when comparing optical to flux data, four weighted NDVI results were calculated according to different footprint configurations: (1) midday average footprint calculated using the sampling area calculated from the footprint model; (2) daytime average footprint calculated using the sampling area calculated from the footprint model; (3) 200-m radius circles centered at flux towers; and (4) a 300 ˆ300 m squares centered on each flux tower; and (5) the 100 ˆ100 m (1-ha) ground sampling region.The latter three footprints allowed us to consider the possible error associated with arbitrary assumptions of footprint shape on upscaling in the absence of an explicit flux footprint model.The radius (200 m) was calculated based on the assumption of a 100:1 fetch-to-height ratio in flux footprint analysis for simple calculation [46,47], and the 300 ˆ300 m represented a square approximation of the same region, whereas the modeled footprints represented a more realistic sampling footprint depictions based on actual flux data (Figure 3).

Vegetation Map Analysis
The vegetation map analysis used a dominant vegetation cover map (1:15,000) at Mattheis Research ranch [39].This map was derived from a combination of aerial photos and field sampling begun in 2010, two years prior to the flight, but published one year following the flight [39].An inherent assumption of our analysis was that dominant vegetation cover was stable over this study period.For validation, 134 5 ˆ5 m plots were selected in the ranch, 66 of which were inventoried by species and percent cover of each species, while only species richness was counted for the other 68 plots.

Biodiversity Estimation
Biodiversity was evaluated three ways: (1) via optical diversity expressed in the airborne imagery; (2) using the vegetation map depicting vegetation types by dominant species (the vegetation types are shown in Figure S1 and Table S1 in Supplementary Materials); and (3) by independent field measurements of species composition for a subset of two of the calibration sites (E3 and E5) (Figure 2).For the first two methods, we used 50 randomly placed circular sampling areas, each having a 200 m radius, along the flight line (Figure 1).Mean reflectance, standard deviation and coefficient of variation of all wavelengths and pixels for each circle were calculated.To provide a metric of optical diversity, we used the average coefficient of variation (CV), calculated as the average CV for all wavelengths from 400 nm to 800 nm (Equation ( 7)).
CV circle " ř 800 λ"400 ˆstdpρ λ q meanpρ λ q ṅumber o f bands In Equation ( 7), ρ λ is the reflectance at wavelength λ, std(ρ λ ) and mean (ρ λ ) are the standard deviation and mean value of reflectance at wavelength λ across all the pixels in one circle, respectively.This metric of optical diversity was applied to each of the 50 randomly distributed sampling circles to compare to productivity estimated from the LUE model, as well as to the entire flight line, using a "moving-window" method, to develop an optical diversity map of the entire region.For the moving-window method, we applied a 10-m lag and calculated the CV for both a circular sampling window (200 m radius) and a square sampling window (1 ha), to illustrate the effect of sampling scale on the resulting optical diversity image.
We calculated the number of dominant species (richness) and Shannon Index of dominant species (Shannon Index) [48] within each of 50 sampling circles as two metrics of diversity from the vegetation map.Shannon index, which combines the effect of richness and evenness, a commonly used measure of species diversity in ecology, was calculated as: where p i is the proportion of the number ith species.Because this map listed dominant vegetation species (and not all species per se), this provided a metric of relative plant species richness by broad community association for each pixel, instead of a full accounting of every species present, which was impractical over such a large region.To confirm that this method correlated with a more complete count of species richness, we also compared the results obtained from the vegetation maps to actual species counts obtained at the two calibration sites in 2015 (assuming a relatively stable species composition between years).Actual species counts were obtained with a subsample of 13 30-cm diameter plots, chosen from a grid covering 1 ha per site.

Model Results
In our study, narrow-band NDVI yielded a high correlation with green biomass and NEE (Table 1), confirming that NDVI provided a useful metric of productivity for this prairie ecosystem.NDVI correlated better with green biomass than with total biomass (data not shown) due to the presence of the prior year's dead or non-green canopy materials, which intercepted solar radiation but contributed little to NDVI.
The relationship between proxy NDVI and net CO 2 flux showed an optimal fit using a 5-h aggregation period around solar noon (data not shown), so this aggregation period was used to develop the NEE model applied here (Table 1).The NDVI-flux relationship from the first half season data (greening phase, before peak biomass) was different from the second half season (July-September senescence phase, after peak biomass) (data not shown).Considering the airborne data was collected in mid-August, which was located around the middle of the July-September senescence phase, data from this second half were used to calibrate the linear model, which yielded an R 2 value of 0.7701 (Table 1).The linear shape of this equation supports the assumption of a constant LUE applied over this time frame for this prairie ecosystem [15,19] and allowed us to use NDVI as an index of productivity.

Sensitivity Analysis of Footprint
The footprint analysis (Figure 3) revealed subtle differences between footprint assumptions that had a relatively small effect on the productivity model.In most cases, the primary source areas for flux measurements were concentrated on wind directions from north-northwest and south for E3 and northwest and southeast for E5.The estimated footprint lengths peaked most frequently between 15 m and 60 m and did generally not exceed 125 m in the along-wind direction.
The footprint sensitivity analysis showed that NDVI weighted was less sensitive to the month of footprint (data not shown) than to the differences between the two sites.On the other hand, assumptions about footprint size and shape (commonly used in multi-scale analyses) had a slightly larger effect on the calculated NDVI than the temporal variation in flux footprint (Table 2).For example, the 1 hectare measurements matching the calibration site areas overestimated the result compared to the modeled footprint analysis by 4% at the E3 site.Two alternative methods of estimating the footprint (200 radius circles and 300 ˆ300 m squares) overestimated the weighted NDVI calculated using the actual percentages of the isopleths at the E5 site (by 2% and 1.7% respectively), but the difference was small at the E3 site (0.6% and ´0.9%) (Table 2).Regardless of the method, E3 had a higher NDVI than E5, showing this site to be more productive.

Diversity Estimation at Calibration Sites
To evaluate optical diversity, we calculated the coefficient of variation (CV) using spectral data from both calibration sites, E3 and E5 (Figure 4).Site E3 revealed a higher overall CV than E5, indicating a higher optical diversity for this site.Similarly, we calculated the CV spectra for each 200-m radius sampling circle indicated in Figure 1, and representative values of high and low optical diversity (CV) (Figure 5a) are shown from this analysis (Figure 4).For subsequent analyses, CV values for each spectrum were averaged across all wavelengths to obtain a single metric of optical diversity for each sampling area, which was then compared to other metrics of diversity and productivity (Figures 6  and 7).1).The positions of low and high richness sites are shown in Figure 5a.The average CV (averaged across wavelengths) provided a metric of optical diversity for subsequent comparison with other metrics of diversity and productivity (Figures 6 and 7).1).Two different sampling scales, a 1-ha square (b) and a 200 m circle (c), were selected for calculating CV as a metric of optical diversity using a sampling lag of 10 m.The vegetation diversity analysis revealed differences between the two sites.Tables 2 and 3 showed the NDVI, vegetation community richness, Shannon index, and field sampling species richness at the two calibration sites (E3 and E5).Species counts from a subset of plots at the two primary calibration sites (E3 and E5) supported this interpretation of biodiversity based on vegetation types; the relative relationship for species richness between the two sites obtained with field sampling matched that of the vegetation map, with higher species richness occurring at the E3 site (Table 3).Similarly, both the mean NDVI and standard deviation of NDVI were higher at the E3 site than at the E5 site, indicating higher productivity and variability in biomass.
Table 3. Biodiversity metrics within a 200-m radius from the flux towers.Richness refers to the number of dominant vegetation types sampled within this 200-m radius based on the vegetation map, and the Shannon Index considered both number and evenness of these vegetation types.Species richness is based on an actual plant species count of a subset of locations within a 1-ha area surrounding the flux towers, as shown in Figure 3.

Extrapolating to a Larger Region
Using data from a larger region (10 km 2 ) of the flight line allowed us to extend the analysis and explore the diversity-productivity relationship over a much larger area.Over the whole flight line, places with higher optical diversity (higher coefficient of variation) tended to have higher mean NDVI, indicating higher green biomass and productivity (Figure 5).The clarity of the diversity images varied with sampling scale, with smaller sampling areas (1-ha square, panel b) providing a sharper image than coarser sampling areas (200-m radius circle, panel c), but both showed similar spatial patterns of optical diversity.
From the analysis of the entire flight line, diversity metrics calculated using the 200-m sampling circles (Figure 1) applied to the vegetation cover map showed that places with higher richness or Shannon index had higher optical diversity, expressed as coefficient of variation (Figure 6), and higher NDVI (Figure 7).The Shannon index showed a stronger correlation with both coefficient of variation and NDVI than with richness in our study (Figures 6 and 7), suggesting an important influence of evenness on these relationships.The fit of the regression for the NDVI-optical diversity relationship can be improved (R 2 increased from 0.46 to 0.58) by applying a logarithmic transformation to the optical diversity data (Figure 7a).In all cases, the relationships between diversity metrics (coefficient of variation, relative richness, or Shannon index) and productivity (NDVI) were highly significant (p < 0.001), supporting a positive relationship between productivity and biodiversity for this prairie landscape.

Green Biomass and NDVI
Like other studies of grasslands [15,19] we found significant correlations between NDVI and biomass.In our study, green biomass showed a clear linear, not exponential, relationship with NDVI, which at first glance seems to contradict results from other grassland studies that have often reported non-linear relationships [19,49].However, this difference was most likely due to the relatively low green biomass values found in our study.The highest green biomass measured in our study was 135.75 g/m 2 , which was a small fraction (2%) of the reported biomass from one study spanning several vegetation types [19] and 1/3 of the maximum green biomass of another study from Alberta prairie that included a relatively wet year [49].The lower green biomass measurement in our study only captured the lower part of the exponential relationship between NDVI and green biomass, which may have appeared linear by missing the rapidly increasing part of what has been described as an exponential function [19].Presumably, if our study spanned many sites and years having higher productivity, a similar non-linear relationship would have emerged between green biomass and NDVI.

Sensitivity Analysis of Footprint
The actual flux footprint varied temporally with wind speed and direction, in agreement with previous studies [50].The exact position of the flux tower or other ground sampling is not an issue when the surface is homogeneous [50].In our study, the footprint area was larger than 1 ha and there was little difference in aggregated NDVI between the circular and square footprint estimates (200 m radius and 300 ˆ300 m 2 ) that were insensitive to atmospheric conditions and the aggregated half hourly footprint model predictions.The footprint analysis allowed us to evaluate potential upscaling errors, and confirmed that, in this case, other approximations of the footprint (e.g., circular or square areas surrounding the flux tower) provide reasonably accurate approximations of the actual flux footprint over a relatively flat and uniform landscape such as the grassland at our site.By contrast, the difference in NDVI values between sites were closer to 20% (Table 2) suggesting that our upscaling method can accurately depict relative differences in productivity across this landscape.These results indicate that, depending upon the degree of landscape heterogeneity, the highly dynamic footprint of eddy covariance measurements [50] can add uncertainty to the calibration, but that any resulting error would have been small for this relatively homogenous site.

Biodiversity-Ecosystem Function
Productivity represents one of the most important ecosystem metrics, and the species richness-productivity relationship has long been of interest in ecology [21,51].Although productivity cannot assess biodiversity directly, it integrates other major drivers of biodiversity that can include land use, climate, nitrogen decomposition, biotic change and increasing atmospheric CO 2 [52].Other prairie studies from the US have shown that high biodiversity plots tend to have higher biomass [53], and our results are consistent with this finding.On the other hand, a recent study of 30 grassland sites around the world (19 countries and 6 continents, including samples from this landscape) suggests that the global diversity-productivity relationship may be more "hump-shaped", particularly when intensively managed ecosystems are included [24].This decline in diversity at high productivity was not seen in our study, and may reflect the greater input of water or nutrients and other management practices associated with many highly managed ecosystems.Compared to many other grassland sites, our study site was a low-productivity grassland.The airborne data offered a means of evaluating the productivity-diversity relationships over a much larger area of this natural prairie ecosystem than would be feasible from the ground, and provided an objective means to address the biodiversity-productivity relationship over large landscapes.
The NDVI-diversity relationship can change through time and be affected by environmental variables and phenology.Cool-season grasses usually dominate the productivity in the early season.Both cool-season and warm-season grasses affect the productivity in the peak season (August) but can be easily influenced by precipitation in summer [54].In our study, NDVI was strongly associated with green biomass, which related to diversity, but it can be influenced by other factors, e.g., canopy growth stage, water stress and flowering patterns.A full evaluation of the seasonal effects on the NDVI-diversity relationship was beyond the scope of this study.However, in a parallel study of a manipulated prairie ecosystem, Wang et al. [55] explored the dynamics of the NDVI-diversity relationship over a full growing season and found that the optimal NDVI-productivity relationship occurred in early August, similar to the time of the overflight in this study."Optical diversity", which indicates the variation in optical signals detected by remote sensing, has been proposed to relate to more conventional metrics of biodiversity [56].Instead of mapping species per se, optical diversity presumably detects different functional and structural properties, which vary with species or functional groups ("optical types") [56,57].This spectral heterogeneity has been related to variation of species or optical types of forests [36,37].In our study, the coefficient of variation was applied as a metric of optical diversity based on information content present in the reflectance spectra themselves, and expressed as the coefficient of variation of reflectance calculated over many pixels.This optical diversity metric showed significant correlations with conventional species diversity indices (richness and Shannon index) (Figure 6), suggesting that optical diversity metrics from airborne remote sensing can provide useful diversity metrics over large regions.This method, based on a statistical assessment of spectral variability, does not directly identify the cause of the strong links with productivity and other diversity metrics.However, it does provide an objective method that can be readily applied to remotely sensed imagery.
The surrogacy hypothesis is sometimes invoked to relate species richness in one taxon, or diversity at one level, to diversity at another level [58].For example, high genetic richness is related to high species richness and environmental heterogeneity is related to species richness.Similarly, in this study, we used the number of vegetation types by dominant species as a surrogate for species richness by assuming species richness will be higher in places with more vegetation types.Our current working hypothesis is that optical diversity is influenced by both leaf traits and canopy structure, as further influenced by the seasonal expression of these leaf and canopy features [57].Because leaf traits and canopy structure vary between species, optical diversity can provide a surrogate (or proxy metric) for traditional metrics based on species richness and evenness.Both optical diversity and species diversity can be influenced by environmental heterogeneity, such as soil texture and microtopographic variability.In our study, subtle gradients in soil and microtopography could have affected the pattern of diversity across the landscape.Fully understanding the mechanisms underlying the optical diversity-biodiversity relationship should be an objective of future work.
The positive NDVI-diversity relationship (Figure 7) is consistent with another recent study using ground NDVI measurement assessing the biodiversity-productivity relationship with a manipulated prairie ecosystem (Cedar Creek Ecosystem Science Reserve, MN, USA) [55].Together, these findings support the positive diversity-productivity hypothesis for these two prairie ecosystems, and are consistent with previous studies of prairies using more traditional field sampling methods [51,53].
More work is needed on the scale-dependence of this method, as sampling scale clearly influences the resulting optical diversity patterns (Figure 5).The coefficient of variation is influenced by the spatial scale ("grain" or sampling size) that relates to the actual vegetation distribution.In our study, we used different sampling scales both in our footprint analysis (Figure 3) and in our depiction of optical diversity (Figure 5).The results at two different scales showed a similar pattern of CV along the flight line (Figure 5).However, it is not yet clear how to decide the "best" scale that balances a large sampling region with the need for fine spatial detail.Multi-scale analysis is needed to investigate the correlation between optical diversity and biodiversity, and remote sensing provides one tool for such an investigation.
Both species richness and evenness influence ecosystem services and the optical diversity metrics measured here.Recent studies have shown that both species richness and evenness have positive effects on the diversity-productivity relationship in ecosystems [59][60][61].In accordance with another recent study [35], our measure of optical diversity showed a better correlation with the Shannon index than richness per se, suggesting that vegetation evenness (or heterogeneity) influences both optical diversity and the productivity patterns detected by NDVI.Adding evenness can add additional information on community structure, which can apparently affect the variance of the optical signal beyond the effects of species richness alone.Further work is needed to explore the exact impact of richness and evenness on the optical signals detected by spectral reflectance.Similarly, more study is needed to understand the relative importance of factors influencing optical diversity that may include canopy structure, leaf traits, and phenology [57].
This study demonstrates a method for integrated analysis of productivity and diversity using a combination of airborne, field sampling and flux measurements.Integrating flux measurements and remote sensing with the LUE model provides a method for assessing ecosystem health and productivity in continuous temporal and spatial dimensions.This can also be a useful approach for evaluating both biodiversity and carbon uptake together, and thus for assessing overall ecosystem health via the provision of goods and services.Airborne campaigns can acquire surface reflectance measurements over large areas without disrupting the flux footprint and can assess ecosystem status over large areas at high spatial resolution.Additionally, airborne imagery can provide help in selecting the ideal position of the flux station to best represent the target ecosystem.However, obtaining high frequency time series of airborne data is still a challenge, largely due to the high cost of airborne acquisition (which can easily run $20-$30K per field campaign).However, these high costs are largely due to the fixed costs of maintaining an aircraft, pilot and flight team over a period of several days.Such costs could be greatly lowered to end users through subsidizing flight costs or through "volume pricing" if groups of users could share the cost of data acquisition, much in the way that satellite imagery can be provided for "free" or at low cost by volume pricing or via government or corporate support.The Google model of providing "free" global imagery to all via Google Earth, and NASA's provision of satellite data through the Earth Observing System (EOS) are both examples of cost-effective business models for remote sensing.

Conclusions
Remote sensing provides an efficient approach to estimating prairie production and biodiversity over large regions.Differences of biomass and ecosystem production across a 10 km prairie transect were shown clearly with airborne images.Regardless of the diversity method used, high biodiversity areas tended to have higher production, in this grassland ecosystem.These relationships were sensitive to both richness and evenness, and the addition of evenness improved the relationship with remotely sensed optical diversity, assessed as the coefficient of variation of reflectance.We propose that optical diversity provides a potent proxy for other more traditional biodiversity metrics (richness and Shannon index).Further work is needed to further understand the proximal drivers and scale-dependence (spectrally, spatially and temporally) of the biodiversity-optical diversity relationships.
This study demonstrates the benefit of coupling traditional field sampling, eddy covariance footprint analysis and airborne remote sensing to estimate rangeland productivity and biodiversity.The CV provides a simple, objective metric of optical diversity that is significantly correlated with other traditional diversity metrics.However, like other commonly used indices in remote sensing of biodiversity studies, we currently lack a full mechanistic understanding of the optical diversity-biodiversity relationship.Future work could also include more extensive airborne campaigns coupled with continuous satellite observation, along with more detailed field studies to detect changing prairie ecosystem function and composition over larger areas and longer time series.At the same time, experimental approaches employing remote sensing methods at multiple spatial, spectral and temporal scales, and across additional ecosystems, are also needed to test the general applicability of the findings reported here.Such experiments are key to developing a defensible operational approach for wide application in prairies for the purpose of assessing rangeland health, production, biodiversity and carbon sequestration.

Figure 1 .
Figure 1.Location of Mattheis Research Ranch (50.9038 ˝N, 111.8799 ˝W) and airborne true-color image with 50 randomly selected 200-meter-radius circles in the flight line used for evaluating the productivity-diversity relationship.E5 and E3 labels showed the two flux towers within this area used for calibrating the productivity map, and inset photographs illustrate the flux tower and phenology station for each site.The Canada map was provided by the Statistics Canada [38].

Figure 2 .
Figure 2. Experimental design and data used in this study.The linear relationship between NDVI, biomass and NEE (Net Ecosystem Exchange) allowed us to use NDVI (Normalized Difference Vegetation Index) as a proxy for productivity (*) that was then compared to several metrics of biodiversity indicated by asterisks (*).

Figure 3 .
Figure 3. Footprints for two flux tower calibration sites (E3 and E5) at Mattheis Ranch in August, 2012 (see "Methods" for details on footprint calculations).Left: Midday (upper) and daytime (lower) flux footprints (black lines) at E3. Right: Midday (upper) and daytime (lower) flux footprints (black lines) at E5. Also shown for comparison are a square (100 ˆ100 m) field sampling region used for field biomass and reflectance calibration, and two alternate footprint assumptions: a circular area (200 m radius), and a 300 ˆ300 m square.Sampling methods were overlaid on a false color image from the airborne imaging spectrometer.

Figure 4 .
Figure 4. Sample coefficient of variation (CV) spectra along the flight line.E3 and E5 indicate the CV spectra of the 200-m radius circle around the two calibration sites.High and low richness spectra indicate the site with highest and lowest CV values of all the 50 200-m radius circles along the flight line (Figure1).The positions of low and high richness sites are shown in Figure5a.The average CV (averaged across wavelengths) provided a metric of optical diversity for subsequent comparison with other metrics of diversity and productivity (Figures6 and 7).

Figure 5 .
Figure 5. Airborne image of NDVI (a) and optical diversity (b and c) along the flight line.The letters "H" and "L" in the flight line (a) indicate the position of high and low diversity sites illustrated in Figure 4. Calibration sites E3 and E5 are also shown.Since NDVI was linearly related to biomass and NEE, the NDVI maps also indicate relative productivity according to those metrics (Table1).Two different sampling scales, a 1-ha square (b) and a 200 m circle (c), were selected for calculating CV as a metric of optical diversity using a sampling lag of 10 m.

Figure 6 .
Figure 6.Optical diversity (coefficient of variation) versus relative richness (a) and Shannon index (b) calculated with the vegetation cover map.Richness represents relative species richness based on dominant species indicated on the map rather than a full field count of all species.Each data point represents a 200-m radius circular sampling area (Figure 1).

Figure 7 .
Figure 7. Optical diversity (coefficient of variation) versus NDVI (a), Mean NDVI versus relative species richness (b) and Shannon Index (c).Both linear and logarithmic fits are shown for panel (a).Each data point represents a 200-m radius circular sampling area (Figure 1)

Table 1 .
Equations derived from field calibration and subsequently used to map green biomass and NEE using airborne NDVI.p values < 0.01 for both relationships.

Table 2 .
Estimated NDVI weighted according to different footprint assumptions.Footprint areas are illustrated in Figure3.