Scale-dependency of Arctic ecosystem properties revealed by UAV

In the face of climate change, it is important to estimate changes in key ecosystem properties such as plant biomass and gross primary productivity (GPP). Ground truth estimates and especially experiments are performed at small spatial scales (0.01–1 m2) and scaled up using coarse scale satellite remote sensing products. This will lead to a scaling bias for non-linearly related properties in heterogeneous environments when the relationships are not developed at the same spatial scale as the remote sensing products. We show that unmanned aerial vehicles (UAVs) can reliably measure normalized difference vegetation index (NDVI) at centimeter resolution even in highly heterogeneous Arctic tundra terrain. This reveals that this scaling bias increases most at very fine resolution, but UAVs can overcome this by generating remote sensing products at the same scales as ecological changes occur. Using ground truth data generated at 0.0625 m2 and 1 m2 with Landsat 30 m scale satellite imagery the resulting underestimation is large (8.9%–17.0% for biomass and 5.0%–9.7% for GPP600) and of a magnitude comparable to the expected effects of decades of climate change. Methods to correct this upscaling bias exist but rely on sub-pixel information. Our data shows that this scale-dependency will vary strongly between areas and across seasons, making it hard to derive generalized functions compensating for it. This is particularly relevant to Arctic greening with a predominantly heterogeneous land cover, strong seasonality and much experimental research at sub-meter scale, but also applies to other heterogeneous landscapes. These results demonstrate the value of UAVs for satellite validation. UAVs can bridge between plot scale used in ecological field investigations and coarse scale in satellite monitoring relevant for Earth System Models. Since future climate changes are expected to alter landscape heterogeneity, seasonally updated UAV imagery will be an essential tool to correctly predict landscape-scale changes in ecosystem properties.


Introduction
Climate change is altering ecosystems worldwide and is leading to a greening trend in the Arctic [1,2]. This includes changes in plant growth, productivity and gas exchange, creating feedbacks to the global carbon cycle [3,4]. Ecosystem changes are typically monitored with very high resolution at small spatial scales by ecological field surveys and experiments using plots (i.e. 0.01-2 m [5];) or at coarse resolution using large spatial scales of satellite imagery appropriate for regional, biome-wide or global estimates and Earth System Models (i.e. 30 m-12 km; [e.g. 1,6,7]. At the plot scale, ecologists have observed numerous changes in Arctic ecosystems in response to warming, including increases in plant height, shrub extent, biomass, growing season length, and decreases in Arctic plant species diversity and bare ground surfaces [4,[6][7][8]. At the coarse resolution satellite scale, a positive Arctic greening trend inferred from the normalized difference vegetation index (NDVI [9];) was observed since the 1980s [1,[10][11][12].
However, the mismatch in spatial scales makes it difficult to link and interpret observed changes derived at either scale. Many ecosystem properties are key biophysical variables in Earth System Models and model benchmarking [13,14], but cannot be directly measured using remote sensing techniques or rely on substantial ground truthing, such as biomass and gross primary productivity (GPP). Biomass is determined through harvesting of in situ living plants, in the Arctic usually at <0.5 m scale [5,15]. It represents a very simple, but fundamental ecosystem variable that relates to carbon storage, respiration and energy processes [16]. GPP is a measure for plant carbon dioxide (CO 2 ) fixation through photosynthesis and measured using flux chambers (0.1-1 m 2 ) over specific vegetation communities or using eddy covariance towers covering heterogeneous footprint areas of ∼100-10 000 m 2 [17]. Spatial estimates of both variables are desired and scaling from field plots to remotely sensed vegetation indices is often performed using NDVI. However, if the resolution of the ground truth product and algorithm does not match the remote sensing product it will cause a scaling bias. This is due both the non-linearity of the relation between variables and NDVI, and heterogeneity within pixels [18][19][20][21]. This scaling biases can be corrected using a range of approaches, but this is especially hard in heterogeneous landscapes and requires sub-pixel information [21][22][23][24][25][26].
UAVs provide a novel resolution advantage compared to satellites, but does higher resolution improve the quantification of ecosystem variables? We investigate whether UAVs are a reliable instrument to replicate fine scale heterogeneity and to follow seasonal changes in ecosystem properties in the Arctic. Second, we analyze whether fine scale heterogeneity in NDVI has a meaningful impact on the estimation of GPP and biomass when scaled to coarser resolution comparable to common satellite platforms. Third, we investigate differences in the observed scaling bias associated with seasonality and location.

Study areas
Four study areas with an extent of 15.4-21.4 ha (∼350 m × ∼500 m) covering extensive tundra terrain were studied near Abisko, northern Sweden (figure 2). The areas are located paired at (low) and above (high) the tree line formed by mountain birch forest (Betula pubescens ssp. tortuosa) and reflect at large the variation in climatic and environmental conditions of vegetated tundra in Fennoscandia. Two areas are located in the Abisko valley (Nisson forest; NF; low/Nisson tundra; NT; high) with a continental climate (1961-90: −0.8 • C; 304 mm), and two are located 30 km northwest in Vassijaure (VJ; forest; low) and Katterjåkk (KJ; tundra; high) representing an oceanic climate with significantly more snow (1961-90: −2 • C; 844 mm). Vegetation at the lower sites consists of alternating tundra, forest and wetland patches and partly vegetated patterned ground. The higher sites mostly feature dry heath tundra with a significant amount of sparsely vegetated barren terrain, patterned ground, snowfields and meadows (see also table S1 (available online at stacks.iop.org/ERL/15/094030/mmedia)).

Research design and plot data
Monitoring was performed in ∼13 d intervals over the growing season in 2017. During each interval, the NDVI was measured using multispectral imagery from a UAV and at 33 ground control plots (GCPs) per study area. Three GCPs were only added midseason. Measurements during each interval stretched over 2-6 d and are subsequently grouped as date clusters.
NDVI was measured at the plots using a handheld Spectrosense 2+ light meter with a footprint of 1 m 2 . Values of −0.2-0.05 generally reflect snow, water and barren environments, while values <0.05 reflect increasing amounts of photosynthetic active tissue in plants [1,47]. NDVI is calculated as a ratio of reflectance in the red and near-infrared (NIR) band: Aboveground biomass was harvested at peak season from representative 0.25 m × 0.25 m (0.0625 m 2 ) plots next to each GCP. GPP at 600 µmol C m −2 s −1 photosynthetic photon flux density (GPP 600 ) was derived from legacy chamber measurements at a scale of 1 m 2 (footprint) unrelated to the GCPs described above [42,48]. It was collected mostly in the Abisko valley and at a more oceanic site between both study areas in June and August in 2004-2005 with corresponding NDVI plot values from a handheld NDVI meter (see [48] for more detail). We fitted exponential regression models to link handheld NDVI to biomass and GPP 600 , respectively. The resulting equations were used to derive maps of biomass and GPP 600 by applying the function on a per pixel basis to NDVI maps generated from the UAV survey. For biomass, this was done at peak season and for GPP 600 over the entire season. See the supplement for more detail on the GCP data collection, the regression and statistical analysis.

UAV and satellite data
UAV data was collected with an eBee fixed-wing mapping drone (senseFly, Switzerland), in combination with a Sequoia multi-spectral sensor (Parrot, France).
The sensor system provides red and NIR imagery enabling NDVI measurements. The UAV was flown at an elevation of 106 m providing a ground resolution of ∼12 cm after photogrammetric image processing using Pix4Dmapper (Pix4D,Switzerland). RGB imagery from a Sony WX and Canon G9X camera was used as complementary data. We retrieved satellite NDVI imagery from Sentinel 2 (10 m spatial resolution), Landsat 8 OLI (Landsat, 30 m) and MODIS (250 m) for 2017. For detailed comparison, we focused on NDVI data for 2017-07-25, as it was the closest date to peak season with cloud free imagery for all study areas and satellite platforms. UAV NDVI data at each GCP was extracted using a 1 m 2 footprint area comparable to handheld NDVI values. See table 1 for a comparison of used data and the supplement for more detail.

Analysis of scale
We investigate the scaling bias when relating NDVI measurements to GPP and biomass from plot scale derived formulas to maps of NDVI by reducing the spatial resolution of the UAV imagery step-wise. We start at the original resolution of ∼12 cm, followed by 25 cm and thereafter in 10 cm intervals: 0.3-10 m, and 50 cm intervals (10.50-75 m). The resolution is reduced using an averaging algorithm for resampling using GdalUtils within R. Averaging preserves the mean and median of the pixels values and best simulates lower resolution remote sensing imagery [54], Table 1.  and it has been used in similar investigations e.g. [18,22,24,26], making our conclusions comparable to previous work. At each resolution and study area, the plot scale derived relationship equations between handheld NDVI and biomass and GPP 600 respectively are applied to the NDVI map on a per pixel basis in order to develop maps of biomass and GPP 600 . The mean NDVI, biomass and GPP 600 are then extracted using the outline of each study area. The calculation of NDVI, biomass and GPP 600 respects the proportional cover of bordering pixels at reduced resolutions and provides consistent results up to a resolution of 75 m limited by the area extent. The workflow is described in detail in the supplement with code examples.

Seasonal NDVI trends across products
The measured mean NDVI values for the individual UAV flights across the season closely match both cloud free satellite imagery and GCP NDVI values (figure 2). Winter and spring values in satellite imagery are dictated by a continuous snow cover (NDVI <0). In the early season (June), when we start to have UAV and hand-held GCP data, NDVI values differ substantially among areas and flights. Early season UAV estimates are lower than satellite NDVI data, but match once snow has melted ( figure 2(b)). At the continental sites (NF and NT), NDVI values were around 0.5 right after snow melt with small patches remaining in depressions. In early June, the oceanic sites (VJ and KJ) were still snow-covered apart from exposed ridges (NDVI: <0.2). Despite the four areas being highly contrasting in land cover (figure 1), the peak NDVI values from the UAV sensors varied only between 0.7 (KJ) and 0.78 (VJ). NF and NT reached peak NDVI already in late July. VJ and KJ reached peak only towards mid and end August respectively. NDVI values drop synchronously to about 0.6 during September. October values from satellite imagery varied substantially among areas, but are likely affected by intense shadow effects.

Plot scale NDVI compared to UAV and satellite derived NDVI
NDVI derived from hand-held sensors and UAVs shows good linear agreement across the growing season ( figure 3(a)). The agreement is weakest in the shoulder seasons and most robust at vegetation peak in August. For NDVI values <0.50, both measurements closely match the expected 1:1 line. At NDVI <0.50, positive and negative disagreement increases. Positive outliers are mostly from barren and water affected land covers. Underestimation is related to fast snow melt between UAV and hand-held measurements in spring. At peak season, the UAV data agrees best with hand-held GCPs, followed by Sentinel 2, Landsat and MODIS ( figure 3(b)). For UAV, Sentinel 2 and Landsat, the agreement is best at high NDVI values and spreads for low values. For MODIS, we see a strong artifact caused by the difference in scale.

Resolving fine scale heterogeneity
The UAV imagery resolves fine scale (<1 m) land cover patches, vegetation heterogeneity and their associated NDVI values not resolved by satellite imagery ( figure 4). Sources of sub-meter heterogeneity are multi-fold: patterned ground in form of nonsorted circles and solifluction lobes of varying size with alternating patches of barren ground and vegetation or flooded center mostly at the continental sites (figure 4 NF, NT); by glacial erratics and partly exposed boulders in all areas, but especially at KJ (figure 4 KJ); mixed, discontinuous vegetation cover (all areas); sparse birch forest at the treeline creates a matrix of high NDVI birch tree canopies over understorey background values (figure 4 NF, VJ); by Compared with satellite imagery, we find that UAV imagery resolves tundra heterogeneity and produces clean pixels of the land cover at the scale of our GCP NDVI, biomass and GPP measurements (0.0625-1 m 2 ; figure 4). At Sentinel 2 scale (10 m), we observe spectral mixing over different land covers, such as barren ground and vegetation patches, erratics, or small ponds. At Landsat scale (30 m), spectral mixing dissolves snowbeds, water paths and blockfields, and edges of decameter scale patches including for lakes and wetlands. At MODIS scale (250 m), all detail is lost and pixels average over several land covers. Furthermore, the spectral value of a given GCP location is dependent on the pixel grid of the product.

Relating NDVI to ecosystem properties
The relationship between NDVI and biomass (gm −2 ) from hand-held sensors at our GCPs is approximated using an exponential model (r 2 = 0.464; p = <0.001; N = 131): The relationship between NDVI and GPP at 600 µmol C m −2 s −1 is approximated as (r 2 = 0.741; p = <0.001; N = 105): Higher r 2 values could be achieved by establishing relationships for individual land cover and vegetation communities. However, a specific vegetation classification is beyond the focus of this study. See the supplement for graphical illustration and confidence intervals of these functions. We also present in figure S2 a comparison of the UAV NDVI to biomass and the hand-held NDVI to biomass relationships confirming that both transfer functions are near identical.

Importance of small-scale heterogeneity
We apply these plot scale derived functions to raster maps of NDVI and calculate the mean value for the NDVI, biomass and GPP 600 , and sequentially repeated this for NDVI maps with a reduced spatial resolution. While the extracted mean NDVI is stable, plant biomass and GPP 600 estimates decrease substantially with decreasing spatial resolution creating an upscaling bias (figures 5 and 6).
The scaling bias of biomass and GPP 600 estimates follows a curve (figures 6; S3; table S2). This curve is initially very steep from 0.1 m to 1 m followed by a strong bending phase until ∼10 m and a tail phase at coarser scales ∼10-75 m. Both, the magnitude of the decline with spatial resolution and the difference in this magnitude among areas are large. The scaling bias is more pronounced at high altitude tundra sites (NT and KJ; figure 6) compared to low altitude sites (NF; VJ), while the oceanic sites (VJ; KJ) show a stronger scaling bias compared to the continental sites (NF; NT). The low altitude sites located at the treeline show similar scaling bias. The high altitude sites differed mainly for coarse resolutions with NT showing a flatter response curve. When the spatial resolution of NDVI matches the ground truth sampling (0.25 m; 0.0625 m 2 ), the biomass is estimated to 549-732 g m −2 , compared to 553-736 g/ m 2 or +0.5%-0.8% at UAV resolution (0.12 m) representing a downscaling bias, but 434-646 g m −2 or -11.5%-21.0% at lowest resolution (75 m) representing an upscaling bias. Likewise, when the spatial resolution of NDVI matched the ground sampling (1 m 2 ) of GPP 600 , it is estimated to 7.65-9.62 µmol C m −2 s −1 , compared to 7.81-9.74 µmol C m −2 s −1 or +1.4%-2.2% Figure 5. Effect of reducing the spatial resolution of UAV imagery exemplified for highly heterogeneous tundra with alternating vegetation cover and barren ground dictated by patterned ground (compare figure 4 NT). Top to bottom rows: RGB imagery, NDVI, estimates of biomass and GPP 600 from NDVI at the specific scale. Note that mean NDVI remains stable with decreasing resolution, while estimates of biomass and GPP 600 drop at lower resolution.

Seasonal patterns of gross primary production
When we apply the plot scale derived GPP 600 function to raster maps of NDVI across the entire growing season, we can show temporal shifts in the upscaling bias (figure 7). Seasonal patterns vary considerably with environmental conditions. The effect is strongest at peak season, while shoulder seasons show flatter response curves. At the continental sites (NT and NF), we captured the first days after snow melt with few remaining snow patches, indicating that at fine resolution (<5 m) the effect and magnitude of curve changes are related to phenology. This is supported by shallower response curves at the oceanic sites (KJ and VJ) early in the season with significant snow cover. At coarser resolutions, the curves are flattest at NT, which is most homogeneous at that scale. The scaling bias differs between early and late season for the oceanic, snow affected sites (VJ and KJ).

Comparing NDVI values: UAV, satellite and ground control data
We show that UAVs with multispectral sensors can monitor seasonal NDVI patterns with high spatial resolution. When aggregated to whole study areas, we find good agreement of NDVI values from UAV, ground control points (GCPs) and different satellites. NDVI curves show an expected trajectory over the growing season with similar peak estimates for all areas ( figure 2(a)), but shorter growing season for oceanic sites. While we found a strong relationship between the UAV NDVI and the satellite NDVI (r 2 = 0.98) when aggregated at the study area ( figure 2(b)), the measured agreement to the ground control decreases with increasing resolution of the NDVI product. GCP NDVI values <0.5 agree persistently better with UAV and satellite NDVI than values <0.5 and UAV NDVI matches GCP NDVI better than satellite NDVI (figures 3(a), (b)). This is clearly an effect of pixel resolution dissolving fine scale land cover heterogeneity ( figure 4). However, UAV data matches the scale of land cover and remaining more nuanced discrepancies to GCPs relate to a number of environmental and technical factors. These include short term environmental changes (rain, snow melt, phenology when GCP, UAV or satellite measurements occur on different dates), illumination differences between dates, and shadows caused by vegetation and topography, and sensor bandwidths (table 1). This is more important for early and late season data presumably due to rapid changing environmental conditions, lower solar angles and more varied weather. We also notice light scattering from vegetated to non-vegetated areas causing overestimation of low NDVI (<0.5) values ( figure 3(a)). For UAVs, low light conditions and wind cause motion blur and shorter flights with less imagery. Some of these effects have been investigated in detail elsewhere and remain an important field of study [30,31,55]. Satellite imagery is less sensitive in many respects, but limited by resolution, flexibility, cost, clear sky conditions, and is affected by shadows and atmospheric interference.

UAV imagery resolves tundra heterogeneity
The UAV imagery reveals numerous environmental sources of compositional (land cover and vegetation types) and configurational heterogeneity (patterns and patches of varying sizes) strongly expressed in tundra terrain (figure 4). Patterned ground caused by periglacial landforms is a widespread Arctic phenomenon [56,57], and has been linked to surface vegetation [58], subsurface SOC variability [59,60] and extensive occurrence of ponds [61]. We also resolve micro-wetlands (<10 m 2 ), small snow patches and boulders, which create spectral contrast. This is not resolved in the Sentinel 2, Landsat and MODIS imagery due to spectral mixing (figure 4) and reflects the fragmented nature of land cover in tundra environments and omission of important ecosystem properties at lower resolutions, such as edges, linear features and patch sizes [36]. However, many of these properties are intrinsically linked to key ecosystem variables and C dynamics.

Scaling bias of biomass and gross primary production (GPP)
The described scaling bias is caused by the nonlinearity in the relationship between remote sensing products and ecosystem properties, and because of the heterogeneity within pixels [18,23]. We find that both these mechanisms are highly relevant in tundra because the GPP 600 ∼ NDVI and biomass ∼ NDVI relationships are strongly non-linear, and vegetation is unevenly distributed. We detected strong scale-dependency with significant under-and overestimation. The upscaling error is up to 17.0% for biomass and 9.7% for GPP 600 at Landsat resolution, and up to 21.0% for biomass and 12.5% for GPP 600 at 75 m resolution. These discrepancies are clearly large enough to matter in a long-term Arctic change context. For biomass, this matches the magnitude of estimated increases of 60-120 g m −2 caused by 5 • C warmer summer temperatures [62] or a measured 19.8% increase between 1982 and 2010 at Arctic scale [63]. For GPP, this is in the same range of the increases found in open top chambers that warm vegetation between 1 • C and 2 • C [64,65].
Previous investigations have shown that such scaling bias exists at meter scale [20]. Using an UAV we found that the scaling bias extents even to centimeter scale, and differs dramatically between areas and seasons. By comparing scaling bias curves, we can extract factors influencing this scale-dependency (figures 6 & 7). At UAV scale (0.12 m), seasonal variability indicates that phenology plays a key role. Reducing the resolution to 1 m dissolves patterned ground, and shrub or tree branches and strong spectral contrasts caused by rocks and ponds (figure 5). At medium scale (1-10 m), scale-dependency is affected by the spatial structure of vegetation communities and micro-topography. At coarser scale 10-75 m, differences point at topography affected land cover patches (tundra, snowbeds, forest patches, lakes) and homogeneous areas turn scale invariant with a flat response curve (NT), while a high heterogeneity from snow patches and vegetation leads to the steep and seasonally contrasting response curves (e.g. KJ at mean GPP 600 ∼ 3-4; figure 7).

Monitoring changes, bias correction and satellite validation
In the context of climate change and Arctic greening, plant communities show a diverse trait response [7]. Since different plant communities typically intermingle at a centimeter to decimeter scale in Arctic tundra, satellite based remote sensing pixels will be a mix of different vegetation types, that will at least partly respond in contrasting ways to climate change. The most straight forward way to scale up results from plot scale (e.g. warming chambers, experiments) to landscapes will be using UAV resolution input data. Further, when we instead use coarse resolution satellite products to estimate ecosystem changes such as Arctic greening, information at UAV resolution is needed for to correct for upscaling biases and interpret the results. Satellite validation of ecosystem variables is important for long-term monitoring [66,67] and should include UAVs in heterogeneous environments to validate statistical relationships derived at coarser scale.
The observed scaling bias can, partly be corrected. A power law relationship between ecosystem properties and scale explains part of the bias, possibly related to fractal properties of land cover and vegetation [68,69]. Bias correction algorithms include using sub-pixel information like land cover and topography [22,24,25], Taylor Series and extensions [21,23,26], and wavelet-fractal technique [23]. Using land cover with specific transfer functions can provide an alternative path to scale ecosystem properties, but also leads to generalization. Upsccaling biases can also be reduced by using vegetation indices less sensitive to saturation than NDVI, such as EVI2 [70]. To be efficient, all of these approaches require information on heterogeneity at sub-pixel scale. However, neither will fully remove scaling biases. Especially in areas where the heterogeneity differs strongly between sites and seasons, and will change in response to climate change, biases will be large and hard to estimate even after corrections have been applied. We thus recommend UAVs for upscaling of ecological processes in such heterogeneous landscapes to get accurate estimates of the focal processes at landscape scale.

Conclusions
Mapping ecosystem properties is key to understand the impact of a warming climate. We show that UAVs provide a great opportunity to monitor seasonal changes in NDVI and to resolve the heterogeneity of tundra terrain at submeter scale. Clearly, ground controls, UAV imagery and different satellite platforms measure the same environmental signal, when aggregated at study area scale (15-21 ha). However, the upscaling bias when linking biomass and GPP 600 from plots to remotely sensed data is at a magnitude comparable to the expected effects of decades of climate changes and we show that it is heavily influenced by the local heterogeneity and seasons. The scaling bias is also relevant to many other ecosystem properties, such as soil respiration, net ecosystem exchange or LAI [20,41,71], since it derives from the non-linear relationship between remote sensing product and ecosystem processes, and spatial heterogeneity. Furthermore, interpretation of biomass and GPP changes from coarse resolution satellite products is limited without high-resolution information on sub-pixel heterogeneity. This is particularly important in the Arctic, due to its highly fragmented land cover. Finally, a warmer climate is expected to influence this sub-pixel heterogeneity and thus the upscaling bias. UAV derived remote sensing products can resolve these issues by providing remote sensing data at the same scale as the ground truth data, facilitate interpretations of the observed coarse pixel changes and be the basis for developing more accurate scaling functions needed to track and validate satellite derived long-term changes in ecosystem properties.