High-Resolution Spatio-Temporal Estimation of Net Ecosystem Exchange in Ice-Wedge Polygon Tundra Using In Situ Sensors and Remote Sensing Data

: Land-atmosphere carbon exchange is known to be extremely heterogeneous in arctic ice-wedge polygonal tundra regions. In this study, a Kalman ﬁlter-based method was developed to estimate the spatio-temporal dynamics of daytime average net ecosystem exchange (NEEday) at 0.5-m resolution over a 550 m by 700 m study site. We integrated multi-scale, multi-type datasets, including normalized difference vegetation indices (NDVIs) obtained from a novel automated mobile sensor system (or tram system) and a greenness index map obtained from airborne imagery. We took advantage of the signiﬁcant correlations between NDVI and NEEday identiﬁed based on ﬂux chamber measurements. The weighted average of the estimated NEEday within the ﬂux-tower footprint agreed with the ﬂux tower data in term of its seasonal dynamics. We then evaluated the spatial variability of the growing season average NEEday, as a function of polygon geomorphic classes; i.e., the combination of polygon types—which are known to present different degradation stages associated with permafrost thaw—and microtopographic features (i.e., troughs, centers and rims). Our study suggests the importance of considering microtopographic features and their spatial coverage in computing spatially aggregated carbon exchange. spatiotemporal of a multiscale datasets. difference index high-temporal along the ﬁlter, across an ice-wedge-dominated the and the daytime net ecosystem exchange the site-speciﬁc relationship NEEday for is aggregated the tower footprint,


Introduction
Climate change has an amplified effect in the Arctic, which has been reported in both observational and modeling studies, e.g., [1,2]. For example, air temperature in the Arctic is expected to increase twice as much as the global average [3]. At the current pace, a significant portion of permafrost-which has sequestered organic carbon over millenniacould be lost in this century [4]. Concern has been mounting over the possibility that increased microbial activity in previously frozen soils could metabolize soil carbon and release greenhouse gases, leading to a positive feedback to climate [5,6]. At the same time, increasing temperature and associated changes in plant species could increase the photosynthetic activity of plants, enhancing the terrestrial CO 2 sink [7].
Ice-wedge polygons, which cover much of the high-Arctic regions in the USA, Russia and others, e.g., [8,9], have unique geomorphic characteristics created by frost cracks in the ground and the subsequent growth of wedge-shaped ice within these frost cracks [10,11]. The growth of ice wedges moves soil and creates microtopography, e.g., [12]. Typically, icewedge polygons have long narrow depressions (i.e., troughs) along frost cracks and above ice wedges. Low-centered polygons (LCPs) have distinct rims or ridges along the troughs, which create low-elevation areas in the centers (Figure 1a). High-centered polygons (HCPs) have elevated centers surrounded by the troughs and are considered to represent the degraded stage of polygon evolution (Figure 1c), after the ridge of the low-center polygons collapses due to melting of the ice wedges. There are also flat-centered polygons (FCPs) which are at an intermediate stage of degradation between LCPs and HCPs. Recently, Liljedahl et al. [13] reported a wide-spread degradation of ice-wedge polygons across the pan-Arctic region associated with permafrost thaw, with the shifts occurring from LCPs to HCPs. This shift could result in increased water drainage and reduced water inundation [13], which could have a significant impact on vegetation-community dynamics, and on carbon exchange [14].
Land 2021, 10, x FOR PEER REVIEW 2 of 19 [10,11]. The growth of ice wedges moves soil and creates microtopography, e.g., [12]. Typically, ice-wedge polygons have long narrow depressions (i.e., troughs) along frost cracks and above ice wedges. Low-centered polygons (LCPs) have distinct rims or ridges along the troughs, which create low-elevation areas in the centers (Figure 1a). High-centered polygons (HCPs) have elevated centers surrounded by the troughs and are considered to represent the degraded stage of polygon evolution (Figure 1c), after the ridge of the lowcenter polygons collapses due to melting of the ice wedges. There are also flat-centered polygons (FCPs) which are at an intermediate stage of degradation between LCPs and HCPs. Recently, Liljedahl et al. [13] reported a wide-spread degradation of ice-wedge polygons across the pan-Arctic region associated with permafrost thaw, with the shifts occurring from LCPs to HCPs. This shift could result in increased water drainage and reduced water inundation [13], which could have a significant impact on vegetation-community dynamics, and on carbon exchange [14].  Extensive studies have sought to characterize the spatio-temporal heterogeneity of carbon fluxes in ice-wedge polygon tundra regions, using (1) chamber-based measurements, e.g., [9,12,[15][16][17][18][19][20], (2) flux towers, e.g., [21] and (3) airborne-based measurements, e.g., [22]. These studies have identified the environmental factors that dictate the spatio-temporal variability of carbon exchange, including soil moisture, vegetation and soil types. In particular, the chamber-based studies have found that these factors are variable on the scale of several meters associated with microtopography, e.g., [9,12]. Wainwright et al. [12], for example, observed that soil CO 2 effluxes are significantly higher in the centers of HCPs than in the ones of LCPs, which suggests that future polygon degradation could lead to an increase in soil CO 2 efflux.
Recent advances in remote sensing-particularly using airborne platforms-have made it possible to capture the spatial variability of the land surface at submeter resolution. This is particularly important for characterizing heterogeneous ice-wedge polygonal ground where ecosystem properties vary laterally over the scales of meters [9,12,23,24]. In addition, unmanned aerial vehicles (UAVs) are increasingly being used to characterize sub-meter topography and vegetation, e.g., [25,26]. Since several studies have established relationships between certain vegetation indices-such as normalized difference vegetation indices (NDVIs)-and carbon exchanges [21,27], high-resolution multispectral optical remote-sensing techniques have been used to facilitate the mapping of carbon exchanges over the ice-wedge polygon tundra. However, airborne campaigns require human control and can be costly, which limits data acquisition to typically one or at most a few times during the growing season.
In parallel, various ground-based sensing technologies have recently been deployed to measure the temporal variability of terrestrial system behavior continuously and autonomously over time. These sensing systems couple various sensors and platforms such as in situ soil and optical sensors, flux towers, geophysics and weather stations, e.g., [24,28]. Such time-lapse sensing and remote data acquisition have been enabled by recent advances in telecommunication technologies and low-cost microcontrollers/computers for automated data acquisition remotely. In particular, flux towers have made a significant contribution to our understanding of temporal (from diurnal to seasonal and annual) changes in land-atmosphere carbon exchanges [29,30]. In addition, robotic technologies have further enabled the development of automated mobile tram systems that can autonomously move and collect terrestrial datasets at high spatial and temporal resolutions. Healey et al. [31], for example, developed such a system with multiple sensors (e.g., normalized difference vegetation index, albedo, radiation), which could move autonomously along a 50-m transect using an above-ground wireline. Such ground-based monitoring systems have a significant advantage in the Arctic region where frequent cloud coverage limits the use of optical satellite images.
In spite of our increasing ability to acquire such diverse datasets using different platforms, significant challenges remain when attempting to capture spatio-temporal ecosystem dynamics continuously over time and space, and to compute the overall or spatially aggregated carbon exchanges. One particular challenge is the integration of different datasets over different spatial scales (i.e., footprints, resolutions, spatial coverages) and sampling frequencies, as well as the integration of direct and indirect information. Chamber measurements are the most direct measurements of carbon fluxes with a footprint of less than one square meter, which can resolve the microtopographic effects. Their spatio-temporal coverage is, however, often limited; typically to several to several tens of locations in the domain and several times in the growing season. The flux towers measure carbon fluxes in a continuous time manner, but their measurements are a spatial aggregate over the tower footprint without resolving microtopography. Although remotely sensed NDVI and other indices can provide information over space and time, they are only proxies for carbon exchanges, and have significant uncertainties. Based on the spatio-temporal coverage of these measurements, it is desirable but challenging to integrate temporally extensive but spatially limited datasets (such as from ground-based sensors) with temporally sparse but spatially extensive datasets (such as from airborne remote sensing platforms). Although there have been several multiscale data-integration approaches available for spatial integration [32][33][34], the integration both in time and space is still at an early stage of development [35,36]. There have recently been several attempts to estimate the spatiotemporal distribution of carbon exchange at the flux-tower or site scale based on remote sensing data and chamber measurements [37] and mechanistic models [38][39][40][41]. However, the new robotics-based ground-based sensor system sensors mentioned above will open opportunities for improving the spatio-temporal estimation.
This study aims to develop a new approach for integrating multiscale multi-type datasets (i.e., airborne images, point measurements and mobile sensors) to quantify the spatiotemporal variability of land-atmosphere carbon exchanges associated with ice-wedge polygon geomorphology. We develop and apply a Kalman filter-based approach to integrate multiscale spatio-temporal datasets. Although the Kalman filter has been extensively applied in various scientific disciplines, including remote sensing and in situ sensor data applications, e.g., [42,43], it has not been used in the context of multiscale and multi-type data integration. The Kalman filter provides a flexible framework that can accommodate multiple datasets with high computational efficiency, as well as integrate the site-specific knowledge and data correlations. In particular, we make use of the recent observations indicating that Arctic tundra properties are highly correlated with one another, due to tightly coupled hydrological, geochemical and ecological processes, e.g., [12,23,44]. Our datasets include airborne LiDAR and multispectral images (which provide submeter information about the variability of topography and vegetation across the site) and a normalized difference vegetation index (NDVI) obtained from a mobile tram system, capable of providing high-temporal information along a 68 m transect during the growing season. Using the Kalman filter, we first estimate NDVI across an ice-wedge-dominated Arctic tundra site over the growing season, and then estimate the daytime net ecosystem exchange (NEEday) through the site-specific relationship between NDVI and chamber-based measurements (collected between 11:00 am-20:00 pm). The estimated NEEday for the area is aggregated over the flux tower footprint, and then compared with the flux-tower NEEday data (averaged between 11:00 am-20:00 pm).
Using the estimated spatio-temporal maps of NEEday, we quantify the difference in NEEday depending on the different polygon types that have been considered as proxies for the different stages of ice-wedge degradation [13,45] as well as microtopographic features (i.e., polygon troughs, rims and centers). We compute the spatially averaged NEEday in each polygon type and feature to investigate the spatial aggregation effects. We demonstrate our approach at the research station developed by the U.S. Department of Energy's Next Generation Ecosystem Experiment (NGEE)-Arctic project located near Utqiaġvik, Alaska, USA. In the last fifty years, the land-cover changes are known to be minimal in this region compared to other Arctic sites [46], although Liljedahl et al. [13] reported some ice-wedge degradation in this region based on ground observations. This site includes different polygon types [12], which provides us with a unique opportunity to test our approach and to evaluate the impact of polygon geomorphology on spatially aggregated carbon exchanges.

NGEE-Arctic Site
The NGEE-Arctic site is located within the Barrow Environmental Observatory near Utqiaġvik (formerly called Barrow), AK, USA. Thaw lakes, drained thaw-lake basins (DTLBs), and ice-wedge polygonal tundra are the main landscape in this region. The NGEE site is located in the interstitial region of DTLBs-as it has been classified by Hinkel et al. [47]-which has been unaffected by the thaw-lake cycle for more than 5000 years [22,47]. Mean annual air temperature in this region is −11.3 • C, and mean annual precipitation is 106 mm [23]. The growing season typically starts from snowmelt in late May or early June, and ends with freezing and snowfall in early-to-mid September.
The site contains different types of ice-wedge polygons ( Figure 2). The polygon types were categorized by Wainwright et al. [12] into low-centered polygons (LCPs), flatcentered polygons (FCPs) and high-centered polygons (HCPs). Each polygon type includes distinct microtopographic features-troughs, rims and centers ( Figure 1). The vegetation characteristics depend on polygon types and features [9,48]. In well-developed ice-wedge polygons, vascular plants (such as Carex aquatilis) typically fill the wet troughs and centers of LCPs. In contrast, mosses and lichens dominate relatively dry areas, such as the rims of LCPs and the centers of HCPs and FCPs. [22,47]. Mean annual air temperature in this region is −11.3 °C, and mean annual precipi-tation is 106 mm [23]. The growing season typically starts from snowmelt in late May or early June, and ends with freezing and snowfall in early-to-mid September.
The site contains different types of ice-wedge polygons ( Figure 2). The polygon types were categorized by Wainwright et al. [12] into low-centered polygons (LCPs), flat-centered polygons (FCPs) and high-centered polygons (HCPs). Each polygon type includes distinct microtopographic features-troughs, rims and centers ( Figure 1). The vegetation characteristics depend on polygon types and features [9,48]. In well-developed ice-wedge polygons, vascular plants (such as Carex aquatilis) typically fill the wet troughs and centers of LCPs. In contrast, mosses and lichens dominate relatively dry areas, such as the rims of LCPs and the centers of HCPs and FCPs.

Remote Sensing Data
We used the same airborne LiDAR data as Wainwright et al. [12], which were collected using a manned aircraft on 4 October 2005. The LiDAR provided a digital elevation map (DEM) over the entire site at 0.5 m by 0.5 m resolution ( Figure 2). We measured the ground surface elevation in September 2011 at 1286 points using a Real-Time Kinematic Global Positioning System (RTK GPS), and found that the root-mean-square error (RMSE) of the Lidar DEM was 6.08 cm compared to the RTK GPS data.
High-resolution Red Green Blue (RGB) images were acquired on 7 August 2013 using a manned aircraft ( Figure 2). The image (at 32 cm by 32 cm resolution) was converted to a map of greenness index (gI) using: where R, G and B refer to the red, green, and blue chromatic channels, respectively. The gI has been used in Wainwright et al. [12] and Dafflon et al. [26] as an effective index to characterize the vegetation variability based on the RGB image. Before the conversion, we used low-pass filtering with a cubic kernel to remove the centimeter-scale details that are beyond the scope of our work to resolve.

Tram and Flux Datasets
The tram system consisted of a motorized cart that carried several sensors including energy/radiation sensors ( Figure 3). The cart was autonomously pulled along a 68 m long traverse, stopping at 0.5 m intervals to collect data at 137 stations. The footprint of each measurement is considered to be on the scale of approximately 0.5 m. Its route spanned different polygon types and features, including HCPs and LCPs. Measurement height

Remote Sensing Data
We used the same airborne LiDAR data as Wainwright et al. [12], which were collected using a manned aircraft on 4 October 2005. The LiDAR provided a digital elevation map (DEM) over the entire site at 0.5 m by 0.5 m resolution ( Figure 2). We measured the ground surface elevation in September 2011 at 1286 points using a Real-Time Kinematic Global Positioning System (RTK GPS), and found that the root-mean-square error (RMSE) of the Lidar DEM was 6.08 cm compared to the RTK GPS data.
High-resolution Red Green Blue (RGB) images were acquired on 7 August 2013 using a manned aircraft ( Figure 2). The image (at 32 cm by 32 cm resolution) was converted to a map of greenness index (gI) using: where R, G and B refer to the red, green, and blue chromatic channels, respectively. The gI has been used in Wainwright et al. [12] and Dafflon et al. [26] as an effective index to characterize the vegetation variability based on the RGB image. Before the conversion, we used low-pass filtering with a cubic kernel to remove the centimeter-scale details that are beyond the scope of our work to resolve.

Tram and Flux Datasets
The tram system consisted of a motorized cart that carried several sensors including energy/radiation sensors ( Figure 3). The cart was autonomously pulled along a 68 m long traverse, stopping at 0.5 m intervals to collect data at 137 stations. The footprint of each measurement is considered to be on the scale of approximately 0.5 m. Its route spanned different polygon types and features, including HCPs and LCPs. Measurement height above the ground surface varied depending on surface properties between 1 m (snow cover) and 2 m depending on microtopography (further information can be found in [49]). We used the NDVI data measured by an upward-and downward-looking spectral Red/NIR sensor. The NDVI data were collected at irregular time intervals (up to 8 times a day, depending on the time of year/season) from 5 June 2014, and to 21 September 2014. We used the mid-day NDVI data collected daily around noon. Although occasional Land 2021, 10, 722 6 of 19 problems arose with instruments or communications, we were able to successfully collect data over 65 days.
above the ground surface varied depending on surface properties between 1 m (snow cover) and 2 m depending on microtopography (further information can be found in [49]). We used the NDVI data measured by an upward-and downward-looking spectral Red/NIR sensor. The NDVI data were collected at irregular time intervals (up to 8 times a day, depending on the time of year/season) from 5 June 2014, and to 21 September 2014. We used the mid-day NDVI data collected daily around noon. Although occasional problems arose with instruments or communications, we were able to successfully collect data over 65 days. In addition to the measurements obtained by the tram sensors, we manually collected point measurements, including soil moisture (top 20 cm) and chamber CO2 flux data, along the tram transect on eight different occasions during 2014. The soil moisture (dielectric permittivity) and fluxes were recorded at 21 of the 137 tram stations. Soil moisture was estimated using dielectric permittivity measurements from a time domain reflectometer. Small-scale net CO2 fluxes between tundra surface and atmosphere were measured with a transparent closed-dynamic chamber connected to a Los Gatos Research, Inc. (San Jose, CA, USA) portable Greenhouse Gas Analyzer in the same manner as in Wainwright et al. [12] and Vaughn et al. [17]. For each measurement, the chamber was placed on a PVC base (30 cm diameter, installed ~15 cm deep) approximately 24 h prior to measurement. Fluxes were calculated from the slope of the linear section of greenhouse gas concentrations versus time. The chamber flux measurements were done during the day mostly between 11:00 and 20:00 (73% samples are between 11:00 and 15:00).
Tower-based eddy covariance (EC) CO2 fluxes (AmeriFlux ID US-NGB [18]) were measured on a 3.75 m tall tripod in the vicinity of the tram setup, employing a Gill R3-50 sonic anemometer (Gill, Lymington, UK) together with a LI-COR LI-7500A CO2/H2O gas analyzer (LI-COR, Lincoln, NE, USA). We followed the standard calibration and data QA/QC procedures as widely applied across the Ameriflux network (ameriflux.lbl.gov, accessed on 7 July 2021). The tower-based NEE was averaged NEE between 11:00 and 20:00 each day. We assume that NEE in our study is representing the daytime average NEE (NEEday).
To compute the contribution of each location to the flux tower, we used a footprint model developed by Kormann and Meixner [50] in determining the spatial extent of the EC system's field of view (i.e., footprint). The footprint depends on the measurement height, surface roughness, atmospheric stability, wind direction and friction velocity. Our In addition to the measurements obtained by the tram sensors, we manually collected point measurements, including soil moisture (top 20 cm) and chamber CO 2 flux data, along the tram transect on eight different occasions during 2014. The soil moisture (dielectric permittivity) and fluxes were recorded at 21 of the 137 tram stations. Soil moisture was estimated using dielectric permittivity measurements from a time domain reflectometer. Small-scale net CO 2 fluxes between tundra surface and atmosphere were measured with a transparent closed-dynamic chamber connected to a Los Gatos Research, Inc. (San Jose, CA, USA) portable Greenhouse Gas Analyzer in the same manner as in Wainwright et al. [12] and Vaughn et al. [17]. For each measurement, the chamber was placed on a PVC base (30 cm diameter, installed~15 cm deep) approximately 24 h prior to measurement. Fluxes were calculated from the slope of the linear section of greenhouse gas concentrations versus time. The chamber flux measurements were done during the day mostly between 11:00 and 20:00 (73% samples are between 11:00 and 15:00).
Tower-based eddy covariance (EC) CO 2 fluxes (AmeriFlux ID US-NGB [18]) were measured on a 3.75 m tall tripod in the vicinity of the tram setup, employing a Gill R3-50 sonic anemometer (Gill, Lymington, UK) together with a LI-COR LI-7500A CO 2 /H 2 O gas analyzer (LI-COR, Lincoln, NE, USA). We followed the standard calibration and data QA/QC procedures as widely applied across the Ameriflux network (ameriflux.lbl.gov, accessed on 7 July 2021). The tower-based NEE was averaged NEE between 11:00 and 20:00 each day. We assume that NEE in our study is representing the daytime average NEE (NEEday).
To compute the contribution of each location to the flux tower, we used a footprint model developed by Kormann and Meixner [50] in determining the spatial extent of the EC system's field of view (i.e., footprint). The footprint depends on the measurement height, surface roughness, atmospheric stability, wind direction and friction velocity. Our footprint model is a crosswind-integrated model, which computes the cumulative normalized contribution and how much of the measured flux comes from a particular distance (along with a particular wind direction), and also estimates the relative contribution of NEE from each pixel.

Integration Strategy
Our goal is to estimate NDVI and NEE over space and time, using (1) the temporally and spatially sparse chamber flux data, (2) the temporally extensive (but spatially limited) tram NDVI data and (3) the spatially extensive (but temporally limited) airborne greenness index (gI) data. Our estimation framework consists of a graphical model connecting multiple datasets at different spatial and temporal scales (Figure 4). Each node represents a variable (i.e., dataset or a target variable), while each connection represents a correlation between two variables [51]. We represent this connection based on the site-specific and general knowledge that a variety of properties are correlated through tightly coupled above/belowground processes described above, including subsurface properties, vegetation, geomorphology and carbon fluxes [12,26]. This estimation framework relies on a two-step approach: (1) the spatio-temporal dynamics of NDVI are first estimated at daily time steps based on the temporally extensive tram NDVI data and spatially extensive gI from the airborne imagery, using the Kalman filter method (described below), and (2) the NEE is then estimated from NDVI based on the correlations that are found between the tram NDVI and chamber measurements.
footprint model is a crosswind-integrated model, which computes the cumulative normalized contribution and how much of the measured flux comes from a particular distance (along with a particular wind direction), and also estimates the relative contribution of NEE from each pixel.

Integration Strategy
Our goal is to estimate NDVI and NEE over space and time, using (1) the temporally and spatially sparse chamber flux data, (2) the temporally extensive (but spatially limited) tram NDVI data and (3) the spatially extensive (but temporally limited) airborne greenness index (gI) data. Our estimation framework consists of a graphical model connecting multiple datasets at different spatial and temporal scales (Figure 4). Each node represents a variable (i.e., dataset or a target variable), while each connection represents a correlation between two variables [51]. We represent this connection based on the site-specific and general knowledge that a variety of properties are correlated through tightly coupled above/belowground processes described above, including subsurface properties, vegetation, geomorphology and carbon fluxes [12,26]. This estimation framework relies on a two-step approach: (1) the spatio-temporal dynamics of NDVI are first estimated at daily time steps based on the temporally extensive tram NDVI data and spatially extensive gI from the airborne imagery, using the Kalman filter method (described below), and (2) the NEE is then estimated from NDVI based on the correlations that are found between the tram NDVI and chamber measurements. Previous studies have reported that NDVI is correlated with NEE in terms of seasonal changes without accounting for diurnal variability [21]. NEE has been estimated from NDVI and other satellite-derived metrics in the Arctic ecosystems [52,53]. In this study, we develop a site-specific function to calculate NEEday from NDVI based on the flux chamber measurements. In addition, we assume that NDVI is correlated with gI because (1) they both represent plant vigor or photosynthetic activities [54], and (2) Dafflon et al. [24] documented the correlations between NDVI and gI at our site. The challenge here is that we have only an estimate of gI at one point in time. With respect to seasonal temporal dynamics, we take advantage of the observations by Dafflon et al. [24] showing that spatial correlations were persistent in the ice-wedge tundra over time such that the lower elevations were consistently wetter and greener (represented by higher gI) over the growing season relative to the higher elevation regions. Additionally, it is known that polygon structures did not change significantly between 2013 and 2014. As such, we use the onetime gI map to represent the spatial variability. Since the correlations between the one- Previous studies have reported that NDVI is correlated with NEE in terms of seasonal changes without accounting for diurnal variability [21]. NEE has been estimated from NDVI and other satellite-derived metrics in the Arctic ecosystems [52,53]. In this study, we develop a site-specific function to calculate NEEday from NDVI based on the flux chamber measurements. In addition, we assume that NDVI is correlated with gI because (1) they both represent plant vigor or photosynthetic activities [54], and (2) Dafflon et al. [24] documented the correlations between NDVI and gI at our site. The challenge here is that we have only an estimate of gI at one point in time. With respect to seasonal temporal dynamics, we take advantage of the observations by Dafflon et al. [24] showing that spatial correlations were persistent in the ice-wedge tundra over time such that the lower elevations were consistently wetter and greener (represented by higher gI) over the growing season relative to the higher elevation regions. Additionally, it is known that polygon structures did not change significantly between 2013 and 2014. As such, we use the onetime gI map to represent the spatial variability. Since the correlations between the one-time gI and daily NDVI changes over time, we compute the temporally variable correlation parameters between gI and NDVI at the tram locations over the growing season.
The Kalman filter algorithm represents a spatio-temporal process by a linear statespace model [55]. We assume a state vector x t = {x i,t | i = 1, . . . , n} representing NDVI at all pixels (the number of pixels is n) at Day t. We denote a data vector by z t = {z j,t | j = 1, . . . , m}, which is a collection of data at Day t, including (1) the NDVI measurements at tram locations (1 ≤ j ≤ m T ; the number of the tram measurements is m T ), and (2) the gI data at pixels (m T + 1 ≤ j ≤ m T + m VI = m; the number of gI data pixels is m VI ). At each time step, the data vector is described as a function of the state vector by the state-observation equation: where H is a m-by-n data transfer matrix representing the connection between data and state vectors, K is the matrix defining the intercept term, and u is a unit vector (a vector containing all elements equal to one with the length of m). Physically, the observation equation (Equation (2)) describes all the data as a function of target variables (NDVI in our case), so that-within the estimation framework-this equation transfers the information from the data to the target variables. H and K include the correlation parameters describing the data value as a function of NDVI at each pixel. At the tram pixels, the data are the direct measurements of NDVI: where w T follows an independent identically distributed (i.i.d.) zero-mean normal distribution with a fixed variance, representing a measurement error. For the non-tram locations, we assume that the j-th gI data is linearly correlated with NDVI at the i-th pixel such that: where a V (t) and b V (t) are the time-varying slope and intercept parameters in the linear correlations, and w V (t) represents the uncertainty in the correlations. We assume that w V (t) follows a zero-mean normal distribution with the variance (σ V ) of the residual in the linear correlations. In Equation (2), the transfer matrix H is represented by: , if j-th gI data corresponds to i-th pixel (m T + 1 ≤ j ≤ m).
The matrix K is an m-by-m matrix where K j,j = 0, if j-th tram station corresponds to i-th pixel (1 ≤ j ≤ m T ) = b V (t), if j-th gI data corresponds to i-th pixel (m T + 1 ≤ j ≤ m).
The parameters a V (t), b V (t), and σ V are co-estimated at each time step during the Kalman filter. At each time step, the linear correlation is evaluated between the NDVI and gI, at the tram locations, and then the correlation parameters and the residual variance are calculated.
The temporal evolution of NDVI-the change in NDVI at all the pixels from date t to t + 1-is described as a transition of the state from t to t + 1 in a state-transition equation as: The matrices A and B determine the magnitude of changes, while the error vector v(t) represents the deviation from the linear model such as natural fluctuation. Although the time evolution of NDVI can be described by a mechanistic model, we use a data-driven, first-order autoregressive (AR1) model in this study, adapted from the spatio-temporal integration approach developed by Chen et al. (2013). In our model, the general form of Equation (7) is changed to: This means that the NDVI at the i-th pixel changes from x i,t to x i,t+1 with the magnitude v i (in Equation (7), the matrix A is an identity matrix, and the matrix B is zero). This implementation is the same as Rastetter et al. [56] who implemented the temporal evolution of leaf area index (LAI) using the Kalman filter. Rastetter et al. [56] assumed that the change of LAI in each time step was small and hence the systematic change component could be approximated to be zero, and that the observation equation (equivalent to Equation (2) above) could transfer the information related to the temporal evolution in LAI from the data within the Kalman filter. Similarly, we assume that the daily change in NDVI is relatively small so that v i follows an i.i.d. zero-mean normal distribution with a fixed variance. The variance of v i is determined from taking the variance of the NDVI difference (x i,t+1 -x i,t ) at all the tram locations. Although the systematic increase and decrease are not included in this formulation, such temporal dynamics are captured based on the tram information represented in Equation (3) with the Kalman filter.
Once the state and observation equations are defined, the standard Kalman filtering algorithm can be used to estimate the spatio-temporal dynamics of x t and its estimation variance (the algorithm is described in Appendix A). The Kalman filter consists of prediction and update at each time step. In the prediction step, the state transition equation (Equation (7)) is used to predict the state vector (i.e., NDVI at all the pixels) at the next time step x t+1 from the previous time step x t . In the update step (or correction step), the estimate of the state vector is improved or corrected by including the tram NDVI and airborne gI data through the observation equation (Equation (2)). In this update step, the information transfers from the data to the target variables (e.g., the estimates of NDVI). The covariance of the state vector is estimated at each step; each diagonal element of the covariance matrix is the estimation variance σ xi,t , representing the uncertainty in the estimates.
After NDVI values are estimated over space and time, NEEday at each pixel is estimated based on the NDVI-NEEday relationship. The NEEday at pixel i and time step t (y i,t ) is a function of NDVI such that: We assume an i.i.d. zero-mean normal distribution to represent the error and uncertainty ε t in this relationship. To account for the uncertainty in this relationship, we used the variance propagation formula, in which the variance of y i,t is |f'(x i,t )| 2 σ xi,t + σ y , where f' is the first derivative of the function f, and σ xi,t is the estimation variance of NDVI at pixel i and time step t derived from the Kalman filter, and σ y is the variance of the residual in the NEEday-NDVI relationship. Once the NEEday map is created, we computed the weighted average of NEEday (Equation (9)) within the EC flux tower footprint to compare with the NEEday from the EC flux tower data. The weights were calculated based on the contribution of each pixel to the tower data in the footprint model.
Our formulation provides a general framework for incorporating multi-type datasets, including time-varying data and stationary data together in the data vector z t . The Kalman filtering approach has been used extensively for remote sensing and environmental science applications as well as for integrating ecosystem or land surface models with eddycovariance flux tower data [39,56]. Our unique contribution is to seamlessly integrate spatio-temporal disparate datasets sampled over a range of spatial and temporal scales, which addresses the particular challenge of integrating spatially extensive but temporally sparse data with high frequency but spatially sparse data.

Results
The tram data reveals the spatio-temporal variability of NDVI along the 68-m transect (Figure 5a), a variability associated with both the polygon features and types which are associated with distinct plant species as described in Section 2.1. After the snowmelt, NDVI increases more significantly at the polygon troughs and centers of LCPs than the polygon rims and the centers of FCPs. The peak NDVI around the day-of-year (DOY) 221 is higher at the troughs and the centers of LCPs. In parallel, the co-located soil moisture and chambermeasured NEEday show their dependency on microtopographic features (Figure 5c,d). Soil moisture does not change significantly over time during the growing season (Figure 5b). However, the spatial variability of soil moisture (top 20 cm) is significant, with the troughs and the LCP centers fully saturated and the polygon rims dry. In Figure 5d, NEE exhibits small spatial variability in the beginning and end of the growing season (DOY 182 and 256), but the spatial variability becomes larger in the peak season (approximately Day 221) with a significantly higher NEE (i.e., carbon sink) at the troughs and the LCP centers. The highest NEE values (Figure 5d) during the growing season correspond to the highest NDVI timing (i.e., DOY 221 in Figure 5b). polygon rims and the centers of FCPs. The peak NDVI around the day-of-year (DOY) 221 is higher at the troughs and the centers of LCPs. In parallel, the co-located soil moisture and chamber-measured NEEday show their dependency on microtopographic features (Figure 5c,d). Soil moisture does not change significantly over time during the growing season (Figure 5b). However, the spatial variability of soil moisture (top 20 cm) is significant, with the troughs and the LCP centers fully saturated and the polygon rims dry. In Figure 5d, NEE exhibits small spatial variability in the beginning and end of the growing season (DOY 182 and 256), but the spatial variability becomes larger in the peak season (approximately Day 221) with a significantly higher NEE (i.e., carbon sink) at the troughs and the LCP centers. The highest NEE values (Figure 5d) during the growing season correspond to the highest NDVI timing (i.e., DOY 221 in Figure 5b).  The tram-measured NDVI and chamber-measured NEEday data are significantly correlated with each other over time (Figure 6a). The NDVI-NEEday correlations at different sampling times are temporally consistent along the same trend. Since some nonlinearity is observed in their relationship, we used the Spearman rank correlation that yielded a correlation coefficient of −0.73 (p-value < 0.01). We fitted a quadratic function to represent this trend curve (R 2 = 0.57 and p-value of <1 × 10 −13 ). On the other hand, the correlations between the airborne-based gI (one-time acquisition) and the tram-based NDVI (multiple acquisitions) change over time (Figure 6b). The Spearman correlation coefficients are −0.10 (p-value 0.26), 0.63 (p-value < 0.01), 0.65 (p-value < 0.01) and −0.01 (p-value 0.92) on Days 182, 206, 221 and 256, respectively. While the correlations are not significant in the shoulder seasons, they are higher in the middle of the growing season when the NDVI and NEE are most spatially heterogeneous.
The tram-measured NDVI and chamber-measured NEEday data are significantly correlated with each other over time (Figure 6a). The NDVI-NEEday correlations at different sampling times are temporally consistent along the same trend. Since some nonlinearity is observed in their relationship, we used the Spearman rank correlation that yielded a correlation coefficient of −0.73 (p-value < 0.01). We fitted a quadratic function to represent this trend curve (R 2 = 0.57 and p-value of <1 × 10 −13 ). On the other hand, the correlations between the airborne-based gI (one-time acquisition) and the tram-based NDVI (multiple acquisitions) change over time (Figure 6b). The Spearman correlation coefficients are −0.10 (p-value 0.26), 0.63 (p-value < 0.01), 0.65 (p-value < 0.01) and −0.01 (pvalue 0.92) on Days 182, 206, 221 and 256, respectively. While the correlations are not significant in the shoulder seasons, they are higher in the middle of the growing season when the NDVI and NEE are most spatially heterogeneous. Using the Kalman filtering approach and the data-derived correlations, we estimated NDVI and NEEday over the growing season between the complete snowmelt (DOY 185) and the first snow accumulation (DOY 250). Although both NDVI and NEE have a diurnal cycle, we estimated only the daytime values corresponding to the chamber measurements. There were missing days in the tram data, owing to instrument, power or communication failures. We interpolated the tram data for each date over the estimation period.
The temporal evolution of NDVI and NEEday over the site (Figure 7) shows a spatial variability consistent with the tram data ( Figure 6). In the early growing season ( Figure  7a,d), both NDVI and NEEday are low across the site, and their spatial variability is small. Both NDVI and NEEday increase significantly in July (between DOY 185 and 210), as does their spatial variability. The troughs and LCP centers have higher peak NDVI and NEEday than the HCP centers and rims in the middle of the growing season (DOY 210; Figure 7b,e). Later in the growing season (Figure 7c,f), both NDVI and NEEday decrease significantly at the HCP centers and rims, while they continue to be high at the LCP centers and troughs. Using the Kalman filtering approach and the data-derived correlations, we estimated NDVI and NEEday over the growing season between the complete snowmelt (DOY 185) and the first snow accumulation (DOY 250). Although both NDVI and NEE have a diurnal cycle, we estimated only the daytime values corresponding to the chamber measurements. There were missing days in the tram data, owing to instrument, power or communication failures. We interpolated the tram data for each date over the estimation period.
The temporal evolution of NDVI and NEEday over the site (Figure 7) shows a spatial variability consistent with the tram data ( Figure 6). In the early growing season (Figure 7a,d), both NDVI and NEEday are low across the site, and their spatial variability is small. Both NDVI and NEEday increase significantly in July (between DOY 185 and 210), as does their spatial variability. The troughs and LCP centers have higher peak NDVI and NEEday than the HCP centers and rims in the middle of the growing season (DOY 210; Figure 7b,e). Later in the growing season (Figure 7c,f), both NDVI and NEEday decrease significantly at the HCP centers and rims, while they continue to be high at the LCP centers and troughs.
We confirmed the performance of the estimation based on the tram NDVI data at the ten locations (selected randomly) excluded in the estimation (i.e., validation data). Validation results showed the root-mean square error (RMSE) of NDVI to be 0.050. In addition, we compared the distributed NEEday estimates (based on the Kalman filter and our chamber measurements) with the flux tower data by computing the weighted average of NEEday within the tower footprint ( Figure 8). The comparison shows that the estimated NEEday closely follows the flux tower data, capturing the seasonal changes and general trend. The correlation statistics between the EC flux tower data and estimated NEEday are R of 0.53 and p-value of <1 × 10 −4 . Since the EC data has large scatter, we smoothed data with the moving average of the 15-day window to represent the seasonal variability. Different from earlier studies (e.g., Rasteter et al., 2010), the tower data were not included in the estimation, and we did not calibrate any parameters according to the flux tower data. We would note that the flux-tower data often have large deviation and scattering (Fox et al., 2008). Land 2021, 10, x FOR PEER REVIEW 12 of 19 We confirmed the performance of the estimation based on the tram NDVI data at the ten locations (selected randomly) excluded in the estimation (i.e., validation data). Validation results showed the root-mean square error (RMSE) of NDVI to be 0.050. In addition, we compared the distributed NEEday estimates (based on the Kalman filter and our chamber measurements) with the flux tower data by computing the weighted average of NEEday within the tower footprint ( Figure 8). The comparison shows that the estimated NEEday closely follows the flux tower data, capturing the seasonal changes and general trend. The correlation statistics between the EC flux tower data and estimated NEEday are R of 0.53 and p-value of <1 × 10 −4 . Since the EC data has large scatter, we smoothed data with the moving average of the 15-day window to represent the seasonal variability. Different from earlier studies (e.g., Rasteter et al., 2010), the tower data were not included in the estimation, and we did not calibrate any parameters according to the flux tower data. We would note that the flux-tower data often have large deviation and scattering (Fox et al., 2008).  We confirmed the performance of the estimation based on the tram NDVI data at the ten locations (selected randomly) excluded in the estimation (i.e., validation data). Validation results showed the root-mean square error (RMSE) of NDVI to be 0.050. In addition, we compared the distributed NEEday estimates (based on the Kalman filter and our chamber measurements) with the flux tower data by computing the weighted average of NEEday within the tower footprint ( Figure 8). The comparison shows that the estimated NEEday closely follows the flux tower data, capturing the seasonal changes and general trend. The correlation statistics between the EC flux tower data and estimated NEEday are R of 0.53 and p-value of <1 × 10 −4 . Since the EC data has large scatter, we smoothed data with the moving average of the 15-day window to represent the seasonal variability. Different from earlier studies (e.g., Rasteter et al., 2010), the tower data were not included in the estimation, and we did not calibrate any parameters according to the flux tower data. We would note that the flux-tower data often have large deviation and scattering (Fox et al., 2008). We computed the temporal average of NDVI and NEE over the growing season, and then took the spatial average within the different polygon geomorphic classes (i.e., the combination of the polygon types and features), shown in Table 1. As shown in the table, between the highest (LCP troughs) and lowest values (HCP centers), the average NDVI varies by 37%, while the average NEEday varies by 183%. The polygon troughs have higher NDVI and NEEday than the other features across the different polygon types. Within the LCPs and FCPs, the trough region has the highest average NDVI and NEEday, while the rim area has the lowest. Within the HCPs, the center region has the lowest NDVI and NEEday. The LCP centers have higher average NEE than the HCP centers by 75%. In addition, we computed the spatial coverage of each polygon feature within each polygon type. The HCP region has a large trough area of 50%, while only 14% of the LCP region is covered by the troughs. The high-elevation region in microtopography associated with low NDVI and NEE (i.e., LCP and FCP rims and HCP centers) occupied 50% of the area in each polygon type. Table 1. Spatial average of NEEday and NDVI temporal average (Av.) over the growing season in each polygon type-feature combination. The standard deviation (SD) is also shown. The spatial coverage is the area of each polygon feature within each polygon type (The total area is normalized to 1).  Table 2 shows the spatially averaged values of each polygon type for the average NDVI and NEEday within. Even though there is significant variability depending on the polygon types and features shown in Table 1, none of the average and peak NDVI and NEEday varies significantly as a function of polygon types in Table 2 (within the standard deviation). The average NEEday is slightly higher in the HCPs than others, even though the HCP centers have lower NEEday ( Figure 5). This is because HCPs has the highest areal coverage of troughs (Table 1) where both NDVI and NEEday are significantly high.

Discussion
In this paper, we have developed an approach to integrate multi-type and multiscale datasets and to estimate the spatio-temporal dynamics of NDVI and NEEday over the site at a resolution that can resolve microtopography associated with ice-wedge polygons. The dynamic maps of NEEday can be compared with the flux-tower data through the weighted average of the NEEday maps over the flux-tower footprint. Our comparison shows a reasonable agreement between the estimates and flux tower data in terms of capturing the seasonal dynamics and trend without any fitting processes. The particular advantage of using the Kalman filter is its ability to account the uncertainty associated with the correlations such that we can determine the estimation variance at each time step. This allows us to use the observations that are only indirectly associated with NEE such as the one-time gI map.
The estimated peak NDVI over the site is 0.3-0.55 within the dry areas and 0.5-0.85 in the wet regions, which is similar to the values found at nearby sites by Engstrom et al. [44] and Beck and Goetz [57]. The peak NEEday ranges from −5.0 to 0.0 umol/m 2 /s across the site, the range of which is larger than that of the flux tower-measured NEE across different drained lake basins in this region [21]. This can possibly be attributed to the spatial averaging effect that occurs for the flux tower data, which is unable to resolve the spatial variability of fluxes associated with microtopography in ice-wedge polygons. Nor can such microtopographic variability be captured by the low-resolution satellite-based estimations of NDVI and NEEday.
This study also provides improved understanding of spatiotemporal variability in ecosystem functioning and carbon exchanges over ice-wedge polygon tundra, specifically through the correlation between NDVI and NEEday, and their relationship to polygon features. Consistent with previous studies, e.g., [9,12], we find that NDVI and NEEday are significantly affected by ice-wedge polygon geomorphology. This variability is primarily dictated by dry features (i.e., polygon rims, and HCP centers) and wet features (i.e., troughs and HCP/FCP centers), owing to the soil moisture limitations in ecosystem functioning, as shown by previous studies, e.g., [9,12]. However, we also find that troughs have significantly higher plant productivity and CO 2 uptake (i.e., high NDVI and NEEday) than other features across different polygon types, and also that the troughs have higher values than LCP centers, even though they both had high soil moisture. This is possibly because low centers tended to have deep ponded water and anoxic conditions that inhibit vegetation growth [13]. In terms of temporal variability, both NDVI and NEEday were low at the dry rims and the center of HCPs throughout the season, while they were more dynamic at the troughs and center of LCPs. Although both NDVI and NEE changed over time, their correlation was consistent along the same curve. In addition, the spatial heterogeneity of NDVI and NEEday was persistent during most of the growing season such that the high NDVI and NEEday regions tended to be the same throughout the season. This is why a single image of gI was able to describe the spatial heterogeneity. This is consistent with the previous findings on the relationship between soil moisture and plant vigor by Dafflon et al. [24], using soil electrical conductivity from electrical resistivity tomography (ERT), and greenness index from pole-and UAS-mounted cameras.
Some of these spatio-temporal characteristics have contributed to the successful data integration. First, the persistent spatial heterogeneity allowed us to use one-time gI data to capture the spatial variability of vegetation vigor, and to use the spatially continuous tram data for quantifying the time-dependent correlations between gI and NDVI over time. In addition, the interlink between the spatial heterogeneity and co-variability among different properties is an important consideration. In general, data requirements and spatial heterogeneity are coupled in spatial statistics or geostatistics such that we need more datasets when the spatial heterogeneity is high [58]. During the early growing season, NDVI and NEEday are relatively spatially homogeneous, even though the data correlation is not significant (i.e., low Spearman correlation coefficients). In the peak season when the spatial heterogeneity is high, the data correlation is high and is able to capture the spatial heterogeneity. This system is suitable for the Kalman filter, which is able to integrate the two models of temporal evolution and data correlations. These two aspects-the data correlation is high to constrain the estimation when the heterogeneity is high-can compensate for one another to reduce estimation uncertainty. In this paper, we focused on the seasonal variability of NEEday.
Having the high-resolution NDVI and NEEday maps enabled us to compute the average and peak values over the growing season, and also the spatial average in each polygon geomorphic class. Our results show that the average NDVI and NEEday vary significantly as a function of the combination of polygon types and features. In particular, the average NEEday varies more than 150% between the highest regions (i.e., LCP troughs) and the lowest regions (i.e., HCP centers). However, the spatial average of NEEday and NDVI (both the average and peak) within each polygon type did not significantly differ, despite the fact that the LCP centers have significantly higher values than the HCP ones, and that the different polygon types were known to have distinct hydrological and ecosystem functioning [12,13,59]. This is because each polygon type includes high and low values associated with the features that are integrated during the spatial aggregation, and because the spatially aggregated effect depends on the spatial coverage of polygon features within each polygon type. This kind of effect is often called the Simpson's paradox [60], in which a trend or significant effect appears in several groups of data but disappears when these groups are aggregated (or vice versa). Typically, to obtain a representative value for each landscape class or zone, the sample average of point-scale measurements is calculated in each class [12,59]. This approach might, however, may miss a certain feature, or may not represent the spatial coverage of each feature, creating a bias in the spatially averaged estimates. At our site, for example, the HCPs tend to be smaller in terms of polygon size, and to have wider troughs where NDVI and NEEday are both high. This fact highlights the importance of estimating the spatial distribution of NDVI and NEE at high resolution over the landscape.
In addition, our results have some implications regarding to the impact of polygon degradation on carbon exchanges. The HCP centers have significantly lower NEEday than the LCP centers, which suggests that NDVI and NEE could decrease significantly at the centers of polygons when permafrost thaws and the LCPs transition to HCPs. At the same time, our results suggest that we need to consider microtopographic features (e.g., troughs, rims, centers) and their spatial coverage. During the polygon degradation, the size of polygons may not change, but the width of troughs may change when ice wedges melt and the polygon rims collapse [45]. To evaluate such an effect, it would be essential to have high-resolution images and DEM for delineating local-scale polygon features and changes. With our approach, it is possible to monitor the changes in carbon exchange associated with the polygon degradation over space and time, using tram data and infrequent (possibly annual or every a few years) high-resolution areal images.
In this study, we did not include the flux tower data into the integration, and used the tower flux data to validate the estimation results based on the integration of the chamber data, tram data and airborne image. It is possible to expand our framework to integrating the flux tower data so that the flux tower data provide additional constraints in the estimation, although it would increase the computational requirements significantly. Such an approach would enable us to include additional predictors such as temperature or photosynthetically available radiation (PAR) typically measured at the tower, and possibly to represent the diurnal variability. The diurnal variability can be captured by the tram measurements throughout the day, although coordinated chamber measurements are necessary. In addition, evaluating data over multiple years will be useful to provide the information on the relationship of NEE with changing air temperature and longer growing season.

Summary
In this study, we first developed the spatio-temporal data integration approach to estimate NDVI and NEE over ice-wedge polygonal tundra during the growing season. The datasets included the spatially extensive but temporally sparse high-resolution airborne RGB image and the temporally extensive but spatially sparse tram data that capture the spatial and temporal variability of NDVI along the 68-m transect. We first applied the Kalman filtering algorithm, which is a numerically efficient spatio-temporal estimation approach, to estimate daily NDVI over the site, and then computed NEEday based on the site-specific NDVI-NEEday relationship. We obtained the daily maps of daytime NDVI and NEEday at 0.5-m resolution over the site (750 m × 750 m). The weighted spatial average of NEEday over the flux tower footprint showed a reasonable agreement with the flux-tower-derived NEEday.
The estimated high-resolution NEEday maps were then used to evaluate the spatial variability of the average NDVI and NEEday during the growing season, depending on the combination of polygon types and features. We found that NDVI and NEEday varied significantly depending on the combination of the polygon types and features over the domain. However, the averaged NDVI and NEEday within each polygon type did not show a significant difference, even though the local effect was significantly different at the centers of different polygon types. This is because high and low NDVI and NEE features co-exist within each polygon type, and that they tend to cancel each out within each polygon type. The averaged NDVI and NEEday were then influenced by the spatial coverage of microtopographic features within each polygon type. In particular, the HCPs and LCPs had similar NDVI and NEEday values on average, even though the LCP centers have significantly higher values than HCPs. This is because the troughs have significantly higher productivity and CO 2 uptake than the centers of LCPs, and the HCP areas have a larger spatial coverage of troughs. Our results suggest the importance of considering microtopographic features for carbon exchange, and the need to integrate multi-type multiscale datasets to capture the microtopographic effects on the spatio-temporal dynamics of NDVI and NEE.