Technical note: A view from space on global flux towers by MODIS and Landsat: The FluxnetEO dataset

The eddy-covariance technique measures carbon, water, and energy fluxes between the land surface and the atmosphere at several hundreds of sites globally. Collections of standardised and homogenised flux estimates such as the LaThuile, Fluxnet2015, National Ecological Observatory Network (NEON), Integrated Carbon Observation System (ICOS), AsiaFlux, and Terrestrial Ecosystem Research Network (TERN) / OzFlux data sets are invaluable to study land surface processes and vegetation functioning at the ecosystem scale. Space-borne measurements give complementary information on the state of the 5 land surface in the surroundings of the towers. They aid the interpretation of the fluxes and support the training and validation of ecosystem models. However, insufficient quality, frequent and/or long gaps are recurrent problems in applying the remotely sensed data and may considerably affect the scientific conclusions drawn from them. Here, we describe a standardised procedure to extract, quality filter, and gap-fill Earth observation data from the MODIS instruments and the Landsat satellites. The methods consistently process surface reflectance in individual spectral bands, derived vegetation indices and land surface 10 temperature. A geometrical correction estimates the magnitude of land surface temperature as if seen from nadir or 40◦ offnadir. We offer to the community pre-processed Earth observation data in a radius of 2 km around 338 flux sites based on the MCD43A4/A2, MxD11A1 MODIS products and Landsat collection 1 Tier1 and Tier2 products. The data sets we provide can widely facilitate the integration of activities in the fields of eddy-covariance, remote sensing and modelling.

observational coverage with a high temporal overlap with most freely available EC records. Landsat measurements are of particular interest because they resolve small spatial details in pixels of 30 m size, but at the cost of missing out on short temporal features. The opposite is true for MODIS data products, which partly average over heterogeneous areas in spatially comparatively coarse pixels of several hundred meters. However, MODIS offers daily, partly even sub-daily temporal resolution. We process EO data sets of both surface reflectance and land surface temperature (LST) for a limited area around a given flux 45 site. For both the quality control and the gap-filling, the approaches aim to be generalisable across all sites without accounting for specific local conditions, yet flexible enough to accurately reproduce phenological behaviour and characteristic features such as disturbances or fast transitions in managed ecosystems. The procedure shall be as simple as possible, computationally efficient and not resort to additional data sources to facilitate a potential application to EO data at global scale.
Observation geometries need special attention as the MODIS instruments measure in a wide swath to obtain high temporal 50 coverage. They scan across their track from right to left with view zenith angles up to 65 degree from nadir. The wide range of viewing geometries leads to different fractions of surface types seen from one overpass to the next for a given site. In addition, vegetation structure and topography, together with the position of the sun relative to the sensors, cause variable shadowing specific information on the quality of the inversion of the bidirectional reflectance distribution function as well as snow cover, 85 platform information and land/water coverage in the scene (Schaaf and Wang, 2015a). pre-processed using the Landsat Ecosystem Disturbance Adaptive Processing System (LEDAPS, Schmidt et al., 2013) and the Landsat Surface Reflectance Code (LaSRC, https://landsat.usgs.gov/landsat-surface-reflectance-data-products) for atmo-100 spheric correction. The pixelQA layer contains information related to clouds, cloud shadows, snow, and ice and is useful for the quality control of the Landsat data (Zhu and Woodcock, 2012;Zhu et al., 2015). The services by Google Earth Engine (Gorelick et al., 2017) provided cutouts of the above mentioned products at the EC sites. Independently of the product and its spatial resolution, the cutout area was limited to a maximum distance of 2 km 105 between a given tower and the centre of a given satellite pixel. Downloading the EO data in tiff-format avoided intransparent reprojection of the data from sinusoidal to regular grid by Google Earth Engine, which would have been problematic for the quality flags in the MCD43A2 and MxD11A1 products. The Landsat data were already provided in regular grid by Google Earth Engine.

110
Data processing works separately for each pixel in a cutout (henceforth subpixel). We describe here the overall concept and rationale of the quality filter and the gap-filling, but report all technical details in the Appendix B.

Quality control and computation of spectral indices
Quality control of the MODIS reflectance-based vegetation indices focused on three aspects: good inversion quality of the 115 bidirectional reflectance distribution function as indicated by the BRDF_Albedo_Band_Quality_Bandx flags in the MCD43A2 product, snow-free conditions according to the Snow_BRDF_Albedo flag, and the omission of reflectance values that are affected by the presence of water in the field of view using the BRDF_Albedo_LandWaterType flag. For the selected data samples which passed those criteria we computed a large set of spectral vegetation indices (table 2). An additional check removed possible values of the vegetation indices outside their defined ranges. Some of the time series contained obvious out-120 lier values. We employed an empirical filter which largely removed those samples which had a particularly large difference to the median of their surrounding values in a temporal window (Papale et al., 2006, technical details on all filters in Appendix B).
In the Landsat data, the flag pixel_qa provided quality attributes (CFMask, Foga et al., 2017) and removed pixels that contained snow/ice, cloud, and/or cloud shadow using a binary flag of presence. Similar to the MODIS product, we computed 125 a series of spectral vegetation indices (table 2) using the good quality observations and removed possible values of the indices outside their defined ranges. A slightly modified filter removed possible outlier values also for the Landsat data (see details in Appendix B.)

Gap-filling
In the literature several gap-filling and smoothing approaches are available which work in one or more dimensions (e.g., Wang 130 et al., 2012; Kandasamy et al., 2013;v. Buttlar et al., 2014;Weiss et al., 2014;Yan and Roy, 2018;Zhang et al., 2021) or use fusion methods between sensors (Verger et al., 2011;Moreno-Martínez et al., 2020). They differ in their levels of sophistication and computational efforts. One of our requirements for the gap-filling approach was that is employs exclusively Appendix B). Consider all times with a snow flag larger than 0.1 or missing snow information as snow covered. The latter periods are included as the snow flag appears to systematically miss snow periods in higher latitudes in the beginning of 145 the winter. Still, frequent gaps with missing snow information also occur during the growing season. In order to avoid wrong filling with a constant value during the growing season this gap-fill step is not applied when the probability of snow cover is low, i.e. when the average seasonal cycle indicates typically snow-free conditions at a given time of the year, or when typically no snow occurs at all at a given site.
3. Subsequently, another moving median in windows of 40 days (4 months for Landsat) fills gaps shorter than 65 days (2 150 months for Landsat).
4. Compute the median seasonal cycle and use it to fill longer gaps by linearly scaling it to the time series in temporal windows. This windowed operation guarantees more flexibility to correctly represent inter-annual variations in the time series and might even partly account for changes in the shape of the seasonal cycle due to disturbances. It is, however, not suited to fill regularly recurring gaps at a certain time of the year, e.g. during rain seasons .

155
5. Fill the remaining gaps by piecewise cubic polynomial interpolation. Time series with less than 300 valid data points in the whole record after application of all the previous gap-filling steps will not be meaningful for analysis but are still filled by nearest neighbour interpolation.
6. Temporal operations cannot meaningfully fill gaps at the beginning and at the end of the record. Therefore the first/last valid data points are repeatedly appended at the beginning/end of the record. For Landsat, the number of available scenes is relatively heterogeneous across the globe (https://www.usgs.gov/media/images/ cumulative-number-scenes-landsat-archive) with some regions having a very good coverage (e.g., North America) while other regions are observed less frequently (e.g., Russia and Africa). Such differences in the availability of good quality data between sites strongly affect the quality of the gap-filling at site level. FluxnetEO therefore provides for each data layer a gap-fill flag 170 which describes whether and if so, how a certain data sample has been imputed which allows users to explore individual sites and use (parts of) the gap-filled data or resort to only using the high quality original data points.
3.2 Preprocessing of MODIS land surface temperature

Quality checks
The quality control of the MODIS LST did not use the flags provided in the MxD11A1 products, but focussed on the removal 175 of outlier values. Negative outlier values in LST might represent residual cloud contamination, whereas unusually high values might originate from undetected saturation in the level 1 data. Empirical quality checks followed the procedure for the MODIS reflectances, i.e. they discarded data points that deviated strongly from the median of their surrounding values in temporal windows of 30 days (Papale et al., 2006). An additional sanity check eliminated any daytime LST that was lower than the minimum of AQUA and TERRA nighttime LST for a given day.

Geometrical correction
For several applications, variable viewing geometries as inherent in the MODIS LST observations are not desirable. A geometrical correction approach developed by Ermida et al. (2018) accounted for directionality in LST retrievals due to vegetation structure and topographical effects. A parametric model estimates the magnitude of LST as if constantly observed from nadir or from an angle of 40 degrees between the sensor and the zenith above a given site. Ermida et al. (2018) derived the coeffi-185 cients for this geometrical model at a resolution of 0.05 degree. We followed the pragmatic approach of selecting the model coefficients for the correction from the pixel containing a given site, and acknowledge that we did not investigate to what extent the given site conditions represent the overall characteristics of the land surface in the allocated pixel. Further input to the geometrical model were the viewing azimuth angles, solar angles at the overpass time and estimates of daily potential radiation at the top of the atmosphere. The geometrical correction was applied to each subpixel in a cutout separately.

Gap-filling
Also for the gap-filling of LST several approaches are present in the literature (e.g., Gerber et al., 2018;Ghafarian Malamiri et al., 2018;Li et al., 2018;Dumitrescu et al., 2020). When using exclusively operations in time and no ancillary data to estimate invalid LST observations, one needs to take care to respect the shorter autocorrelation of LST compared to the reflectance-based indicators. According to Vinnikov et al. (2008) the weather-related component of clear-sky LST has an autocorrelation of about 195 3 days. The following sequence of steps filled the four MODIS LST data streams (for technical details refer to Appendix C): 1. Similar to the reflectances, a first step consisted in a temporal moving median to fill gaps, but in shorter windows of eight days. Li et al. (2018) and Crosson et al. (2012) and foresaw to use one of the four MODIS LST time series as a 'reference' to fill gaps in a second 'imputed' one. We computed a median seasonal cycle of the difference 3. In fully cloudy days without any valid LST observation, or in case a period has too few valid observations for a meaningful calibration of the linear model in the previous step, the gap-filling followed the same steps like for the reflectancebased spectral indices:

A second step was inspired by
Linearly scale the valid LST observations of each of the four data streams to their own median annual cycle in temporal 210 windows.
4. Interpolate the remaining gaps with cubic polynomials, or nearest neighbour in case of very low data availability (less than 300 valid data points in an entire time series).

5.
Missing values at the beginning and the end of the record cannot be meaningfully filled by temporal methods and are therefore simply repeated.

215
Steps 3-5 produced very smooth and therefore less realistic LST estimates than steps 1-2. Also, one needs to be aware that any LST estimate in data gaps from this procedure necessarily represents an LST estimate under clear sky conditions, which can be very different from the real LST under overcast skies (Ermida et al., 2019). This needs to be considered for a given application to prevent effects of clear-sky bias in the LST data sets on the results. Like for the vegetation indices, also the LST data layers have a gap-fill flag in FluxnetEO describing which data points are original and which gap-filling step filled the missing values.

Gap-statistics across indices
Data availability after quality screening is highly variable between sites and depends on the data stream ( Fig. 1). MODIS LST generally has less valid data points among the data sets than the reflectance-based indicators, and often less during daytime than nighttime. While the LST are instantaneous values, the reflectances represent averages over 16-day periods. A lower 225 number of good quality observations in indices that rely on band 6 relate to degraded detectors in AQUA MODIS band 6.
Large differences in the amount of good quality data between groups of plant functional types, especially for the reflectances, mirror general atmospheric conditions in different regions.

Temporal patterns of the gap-filled time series
We illustrate some characteristics of the MODIS time series in FluxnetEO using example sites. The Austrian site Neustift 230 (AT-Neu) was situated in a valley in the Alps and surrounded by grasslands which were typically mown three times a year (Wohlfahrt et al., 2008). According to their nature, the LST time series exhibit faster variability than the vegetation indices  years to produce smooth transitions between the good quality data and the gap-filled ones, the interpolation is able to preserve inter-annual variations in the EVI.
Missing LST values were estimated most reliably in the gap-filling steps 1-2 (moving median and scaled average shift to observations at other overpass times) because the typical short-term variability in the time series could be preserved. In the  is unrealistically low because the gap-filling needs to resort to filling by a median seasonal cycle of LST (obtained from those years in which the monsoon starts late) or by interpolation.

On the importance of spatial context
The type and distribution of the vegetation around a given EC measurement station are not necessarily homogeneous. Instead, clusters of different vegetation or land use types might prevail in different sections of the immediate surroundings of a site.

275
The area that a given flux measurement is representative of (the flux footprint, Schmid, 1997) changes rapidly with wind direction, turbulence conditions, atmospheric stability, and surface resistance (Schmid, 1997;Vesala et al., 2008;Chu et al., 2021). An exact match between the flux footprint and EO data (or a model grid cell) is challenging due to the often unknown or uncertain flux footprints and coarse spatial grid sizes. The scale mismatch is equally important in validation exercises for site-level measurements of surface reflectance (Román et al., 2009;Cescatti et al., 2012), site-level energy-balance closure 280 (Stoy et al., 2013) and model-data integration (Williams et al., 2009). The role that the scale-mismatch between site-level and EO data plays for ecosystem analyses clearly depends on the site and the application. Some applications try to account for the mismatch (Pacheco-Labrador et al., 2017;Wagle et al., 2020), others ignore it and use a custom area around each EC site.
Approaches to quantify and account for heterogeneity within a satellite pixel or a certain area around a given site do exist in the literature (Román et al., 2009;Chu et al., 2021;Duveiller et al., 2021), but seem less exploited. In this section, we present (pixels >=20 m). While the trees are evergreen, the herbaceous layer senesces in summer and re-greens in winter (Luo et al., 2018). The EO cutout of 2x2km 2 includes irrigated agricultural areas north of the flux footprint. These fields are barren in 295 winter and are covered with crops in summer. MODIS and Landsat EVI are strongly negatively correlated to GPP derived from EC in the pixels over agricultural areas, as are the anomalies of EVI and GPP (Fig. D1 a-d). Conversely, high positive correlations prevail across the remaining larger parts of the EO cutouts. Landsat EVI overlaid by the average flux footprint for two example months illustrates that the EC GPP is only representative of the tree-grass ecosystem (Fig. 5e, g). Hence, the spatial representativeness of EO data for EC fluxes might differ strongly depending on which satellite pixels are chosen 300 for the analysis. We computed the average EVI that is representative of the flux footprint (henceforth fpa for footprint area).
We compared it with an average EVI weighted with the probability density function of the flux footprint in order to take into account the decreasing influence of subpixels further away from the tower (henceforth fpw for weighted footprint area), as well as with two pragmatic approaches in case a flux footprint is unknown: an EVI average over all subpixels in the cutout with a radius of 2 km (henceforth fex for full extent) or only the single subpixel that contains the tower (cpx for center pixel). The 305 most noticeable difference between the time series for the different intersection methods is that the full extend (fex) in both while the footprint intersection methods (fpa and fpw) and the center pixel (cpx) EVI consistently indicate high greenness in the tree-grass ecosystem.
Gebesee, DE-Geb, is an agricultural site. The common approach in conducting EC measurements is to put the tower in a loca-310 tion where the land use is as homogeneous as possible, to be able attribute fluxes to a targeted ecosystem, e.g. a known crop type. In Gebesee, this was assured for most of the years in the long site history (e.g. Fig. 5h), but not from 2011-2013. In these years, the field was split into two different adjacent crop types that contributed to the measured fluxes (Fig. 5f), raising the risk for pitfalls in the analyses of the fluxes. Also, in situations/ years when the flux footprint represents a single field, additional potential difficulties originate from phenological differences between fields within the EO cutouts ( Fig. 5f,h) if not properly 315 matched. For example, the anomalies of both GPP and EVI are only highly correlated with each other in the immediate surroundings of the tower (Fig. D1g-h). Phenological heterogeneity between fields might explain why the EVI averaged over the full cutout (fex) is clearly different from the EVI in the footprint area (fpa, fpw) or the tower pixel (cpx) during the growing season maxima in 2015/16 (Fig. 5b,d). Also, consistently with the GPP, the EVI in the tower pixel indicates slightly later senescence in 2017 than averaged over the footprint area or the full cutout, highlighting considerable effects of a mismatch 320 between the flux footprint and the EO area.
Irrespective of the match between flux footprint and the area that the EVI is representative of, Fig. 5 illustrates the complimentarity between MODIS and Landsat in terms of resolution. Although Landsat offers high spatial detail, the temporal patterns that can be resolved with monthly averages are much coarser than the shorter variations that daily MODIS data can describe.
Depending on the application the user of FluxnetEO might choose one or the other.

325
RU-Zo2, the Zotino tall tower observatory ZOTTO, is located in the taiga-tundra transition zone. The landscape in the proximity of the EC station is a heterogeneous mix of forest, bogs and wetlands. At the tall tower, fluxes are measured at different heights above the canopy. The size of the flux footprint strongly increases with height and the fluxes at the highest level partly represent areas more than 2 km away from the site (Fig. 6b-d). Flux footprints of measurements closer to the canopy are usually 330 much smaller than the MODIS pixel size of 1 km for the LST, but the flux footprints of the higher measurement levels at RU-Zo2 partly integrate over multiple of such pixels. Size and direction of the footprint extents strongly vary over time (note that Fig. 6b-d represent three consecutive days), such that the vegetation types and surface conditions sampled do not only differ between measurement heights but also between days. We compare spaceborne LST AQUA day integrated over the flux footprint area (LST fpa ) with surface temperature inverted from long-wave outgoing radiation measured at the tower for clear-sky days 335 (Fig. 6a, see details about the methods in Appendix D). LST fpa of all three heights is consistently about 30% higher than the inverted surface temperature for most of the year, with a notably higher scatter under freezing conditions. The slope between LST fpa and surface temperature markedly decreases for the highest temperatures, which might indicate significant changes in surface emissivity during the brief peak growing season when vegetation extent is highest and the surface has drained from snow melt.
The proposed methods aim at assuring good quality and producing as reliable as possible gap-free estimates of EO-derived surface reflectance, vegetation indices, and LST for pixels around EC sites. Depending on the question/ application at hand, MODIS or Landsat EO data with their inherently very diverse spatial and temporal resolutions might be more suitable. The 355 requirements for the strictness of the quality checks and the sophistication of the gap-filling methods differ by use case. No approach can fit all requirements, but we expect FluxnetEO to offer many opportunities to advance our understanding of landatmosphere fluxes for individual sites, across regional networks and globally. It helps bridging the Fluxnet, remote sensing, and modelling communities, and facilitates consistent benchmarking of EO-based flux models of any kind. We anticipate that this will accelerate our ability to monitor and understand land-atmosphere fluxes across spatial and temporal scales. For the 360 future we plan to maintain, update and improve FluxnetEO. This will include extending the time series to most recent years, adding EC sites as measurements become available in one of the networks, improving the processing based on newly identified drawbacks and/ or user needs (e.g., Landsat sensors harmonisation), and updating to new EO data collections (e.g. Landsat collection 2). Importantly, forthcoming FluxnetEO versions shall more strongly facilitate complementary usage of multiple missions to exploit their synergy potential, so that future additions will include further EO products, for example the Sentinel eu/collections/tEAkpU6UduMMONrFyym5-tUW). Zipped folders package the data by continents and groups of countries. In the zip-370 directories, the files are organised by site and in two processing versions: One version contains spatially explicit data fields for each subpixel in the cutout of 2x2km 2 and is denoted by 'subpixel' in the file name. A second version is an average time series per site that represents the area within 1km radius of the site with the inverse distance to the tower as weight ('average_cutout'). In this version, at every time step all valid subpixels closer than 1km to the site are averaged after the quality checks, and the gap-filling procedure takes this average time series as input. The data fields contained in both processing versions are listed in table 2. Each data field has a complementary data layer with a 375 flag ('gapfilltype') indicating which data point is of original good quality or how a given point has been imputed in the gap-filling procedure.
The processing version 'average_cutout' has additional fields that indicate how many valid pixels within 1km of the tower contributed to the spatial average per time step ('N') and the spatial standard deviation of the vegetation index or LST for the given time step ('NSTD').  -15.4427, 167.192 ZA-Kru -25.0197, 31.4969 ZM-Mon -15.4378, 23.2528 In this section we provide all specific technical details necessary to reproduce our processing steps for the surface reflectance of MODIS and Landsat.
The quality control of the MODIS reflectance-based land surface indicators included the following steps: -Omission of the MCD43A2 BRDF_Albedo_Band_Quality_BandX flags ≥ 3 for each band to remove bad inversion 385 quality from the surface reflectances.
-The flag Snow_BRDF_Albedo eliminated pixels that contain snow. As the gap-filling procedure used the snow information, a spatially aggregated snow flag was needed for the processing version that averages valid data within 1 km of the tower. For this, we defined the aggregated snow flag as the fraction of subpixels in the cutout that are snow covered. If more than 50% of subpixels have missing snow information for a certain day, the aggregated snow flag is set to missing 390 as well.
-The presence of water in a scene seen by an optical sensor can strongly affect the observation. The BRDF_Albedo_LandWaterType flag allowed to filter for pixels exclusively on land (flag=1). This eliminated all data for many Swiss, Dutch, Italian and Finnish sites which are situated close to water bodies. Inclusion of ocean coastlines and lake shorelines (flag=2) and shallow inland water (flag=3) resulted in reasonable time series at most sites. This came at the cost of having few other 395 sites that were affected by the presence of water. As a trade-off between data availability and quality, we decided to include land-water flags 1-3.
-After the computation of the vegetation indices from the individual spectral bands, an additional check removed possible values of the spectral vegetation indices outside their defined ranges. An outlier filter compared each value to the median of all valid values in temporal windows of 30 days (Papale et al., 2006). A large difference of a given value to the median 400 of its surrounding values indicates a potential outlier. The threshold z as in Papale et al. (2006) was set to 2, and only a less conservative threshold of z=3 acted when more than 20 valid values were available in a given window.
The empirical outlier filter for Landsat slightly differed from the one for MODIS and removed observations in the five highest and lowest percentiles of the mean seasonal cycle of an index if they differed more than 75% from their surrounding 3-months moving window median. The second criterion was critical in order to preserve observations of disturbance events or 405 recovery dynamics.
Technical details for the gap-filling: 1. The first step is a moving window median to fill short non-snow related gaps. If the entire time series has less than 40% valid data, a given moving window contains both the actual values and the median seasonal cycle for the given time of  periods with missing snow information in the Snow_BRDF_Albedo flag needed special treatment. Some of these gaps appeared systematically in early winter in higher latitudes, so also times of missing snow information are considered as snow covered. However, also during the growing season long periods of missing snow information occur in several sites 415 globally. The following criteria check whether a period that is considered snow covered by high values or missing snow flags is filled with a constant baseline value or not: -If a given site has less than 60 days with valid snow coverage (i.e. Snow_BRDF_Albedo=1) in the total record, snow typically does not occur at the site. In this case the gap-filling procedure does not apply this gap-filling step at all for this site.

420
-The gap-filling with a constant value only addresses gaps with a minimum length of 20 consecutive days with snow flag missing or 1. This avoids filling very short intermittent snow periods or short gaps in snow information during the growing season.
-This gap-filling step does not consider gaps due to missing snow information if the median seasonal cycle of snow coverage indicates ≤ 5% of snow cover at the given time of the year and the difference between the fill value 425 and the median seasonal cycle is large (i.e. exceeds the 85 th percentile of the differences in times of missing snow information).
The constant baseline value that is used to fill snow periods in the time series for a site represents the 3 rd percentile of the median seasonal cycle of the spectral vegetation indices. If a given index typically has high values outside the growing season, the baseline value represents the 97 th percentile instead. However, if for a given winter the average over the last 430 5 valid data points at the end of the growing season or over the first 5 valid data points at the beginning of the next growing season is lower than the baseline value (higher than the baseline for indices which are typically high outside the growing season), the baseline takes the value of this average for the given winter (similar to Beck et al., 2007).
3. Linearly scale the median seasonal cycle to the time series to fill longer gaps . Calibration happens in moving temporal windows of 80 days, and application of the scaling in steps of 20 days.

435
Appendix C: Technical details about the processing of MODIS LST In this section we provide all specific technical details necessary to reproduce the processing steps for the MODIS LST.
The empirical filter to remove potential outlier values (Papale et al., 2006) followed the same procedure like for the vegetation indices, but used a constant z-value of 1.5 as it provided the best trade-off between filter success, wrong positives and wrong negatives.

440
Estimates of LST in data gaps originate from the following steps: dows of 8 days is not supplied by the median seasonal cycle in case of low data availability. The moving window median was not applied for windows with less than three valid values.
-Filling by linearly scaling the median seasonal shift between any two of the four MODIS LST time series to each other 445 (Crosson et al., 2012;Li et al., 2018). The following explains this gap-filling step for TERRA day as the 'imputed' time series: 1. Obtain the median seasonal cycle (MSC) of the shift between TERRA day and AQUA day : MSC( ∆(TERRA day , AQUA day ) ).
2. Linearly scale MSC( ∆(TERRA day , AQUA day ) ) to ∆(TERRA day , AQUA day ) in temporal windows of 80 days 450 (provided a minimum of 10 valid values in a given window). Apply the scaling in windows and steps of 20 days.
∆(TERRA day , AQUA day ) t=k:k+80 = f ( MSC( ∆(TERRA day , AQUA day ) ) t=k:k+80 ) ∆(TERRA day , AQUA day ))* t=k:k+20 = m · MSC( ∆(TERRA day , AQUA day ) ) t=k:k+20 + n. served to fill gaps in TERRA day , namely in the order of increasing standard deviation of the differences between valid TERRA day and each of the three estimated TERRA day *.
The procedure analogously filled AQUA day , TERRA night and AQUA night accordingly using valid observations of the 460 remaining three, respectively.
-Linearly scale the MSC of one LST time series to the actual time series in temporal windows. As in step 2, the calibration happened in temporal windows of 80 days, while the scaling was applied in windows of 20 days. Exemplarily for TERRA day : TERRA dayt=k:k+80 = f ( MSC( TERRA day ) t=k:k+80 ) TERRA day * t=k:k+20 = m · MSC( TERRA day ) t=k:k+20 + n 2020) using the R-code version (V1.41) of the FFP-tool. As a flux footprint for the intersection with EVI we define the area that contributes 80% to the flux footprint probability density function (80% isoline of the monthly/daily cumulative flux footprint for Landsat and MODIS, respectively).
Flux footprint calculation followed the same procedure for the three measurement heights at RU-Zo2. Surface temperature was inverted from long-wave outgoing radiation measured at a fixed height of 302 m using Stefan-Boltzmann law. As the inverted 475 surface temperature was compared to LST AQUA day , the average of half-hourly outgoing long-wave radiation for the nominal overpass time at 1.30pm ± 1.5 hours was taken. Surface emissivity is unknown and we assumed emissivity=1 throughout the year. Only days with good quality in both the LST and the long-wave outgoing radiation are used according to the following criteria: i) more than 90% of the EO cutout have valid (i.e. non-gapfilled) values which restricts the comparison to clear-sky conditions, and ii) at least 50% of the half-hourly long-wave fluxes in a given day are of good quality. A larger cutout of 480 5x5 km 2 was extracted for MODIS LST to fully cover also the extent of the flux footprint of the highest measurement level, but is used only for illustrative purposes and not in the data provided in the FluxnetEO collections.