Effects of Measurement Height and Low-Pass-Filtering Corrections on Eddy-Covariance Flux Measurements Over a Forest Clearing with Complex Vegetation

Flux measurements over heterogeneous surfaces with growing vegetation and a limited fetch are a difficult task, as measurement heights that are too high or too low above the canopy adversely affect results. The aim of this study is to assess implications from measurement height in regard to low-pass filtering, footprint representativeness, and energy balance closure for a clear-cut site with regrowing vegetation of varying height. For this, measurements from two open-path eddy-covariance systems at different heights are compared over the course of one growing season. Particular attention is paid to low-pass-filtering corrections, for which five different methods are compared. Results indicate significant differences between fluxes from the upper and lower systems, which likely result from footprint differences and an insufficient spectral correction for the lower system. Different low-pass-filtering corrections add an uncertainty of 3.4% (7.0%) to CO2 fluxes and 1.4% (3.0%) to H2O fluxes for the upper (lower) system, also leading to considerable differences in cumulative fluxes. Despite limitations in the analysis, which include the difficulty of applying a footprint model at this study site and the likely influence of advection on the total exchange, the analysis confirms that information about the choice of spectral correction method and measurement-height changes are critical for interpreting data at complex sites.


Introduction
With the establishment of the eddy-covariance technique over wide ranges of ecosystems and applications, more research is focused on sites with non-ideal heterogeneous characteristics B Oliver Reitz oliver.reitz@geo.rwth-aachen.de 1 Physical Geography and Climatology, Institute of Geography, RWTH Aachen University, 52062 Aachen, Germany 2 Agrosphere Institute, IBG-3, Forschungszentrum Jülich GmbH, 52428 Jülich, Germany (Griebel et al. 2020;Chu et al. 2021). This challenges basic assumptions of the method and can imply unknown modifications on measured fluxes as well as reducing the energy balance closure (Stoy et al. 2013). Furthermore, internal boundary layers form over surfaces with heterogeneous roughness or a limited fetch. It has been shown that they alter wind and friction velocity (u * ) profiles depending on the type of roughness transition (Jegede and Foken 1999;Dellwik and Jensen 2005). However, most footprint models, such as that after Kljun et al. (2015), do not consider these non-ideal conditions inducing horizontally heterogeneous flow.
A particular difficulty is the growth of vegetation, which demands a regular adjustment of the measurement height to ensure flux contributions from the same source area (Munger et al. 2012). Sensors too far above the canopy are susceptible to measured contributions from fluxes originating outside the area of interest when the fetch is limited (Gash 1986;Nicolini et al. 2017). In addition, steep roughness changes, such as forest edges, can induce recirculation areas behind the edge (Detto et al. 2008), further constraining the available fetch for measurements in forest clearings. On the other hand, multiple issues can result from measurements too close to the canopy. Measurements within the roughness layer may not be representative of the average ecosystem, rather sensors detect a near-field contribution of individual roughness elements leading to flux biases (Katul et al. 1999;Moureaux et al. 2012). Over inhomogeneous surfaces, a small source area resulting from a measurement height below the blending height can also induce a location bias, which is not representative of the average ecosystem flux (Schmid and Lloyd 1999). Lastly, spectral attenuation in the high frequency range, also called low-pass filtering (LPF), is expected to increase with a lower height of eddy-covariance sensors. This increase occurs because smaller eddies, which prevail closer to the ground, are more attenuated than larger eddies by individual LPF causes. For open-path systems, these causes mainly are sensor separation, time response, and path averaging (Burba 2013). A number of correction schemes exist to compensate for LPF, such as fully analytical methods modelling individual sources of attenuation (e.g., Moore 1986;Moncrieff et al. 1997;Massman 2000), in situ methods incorporating scalar spectra (e.g., Ibrom et al. 2007;Fratini et al. 2012), and fully experimental methods using the ratio of cospectral densities (e.g., Su et al. 2004;Polonik et al. 2019). Despite the availability of a variety of approaches and corresponding software packages, which in turn have different options, only a few studies have conducted a comparison of different LPF corrections (Fratini and Mauder 2014;Polonik et al. 2019), and no comparison of multiple corrections exists. Different LPF corrections are deemed suitable for specific set-ups. For example, Ibrom et al. (2007) proved a good performance of their method for measurements taken high above a rough forest surface, while Fratini et al. (2012) showed improvements of their method for measurements taken low over a smooth clover field. For measurements taken low over a regrowing clear-cut area, the contribution of high-frequency turbulence likely is important too, but in such a case the surface is comparatively rough. Thus it has to be investigated if the method of Fratini et al. (2012) also performs well compared to other corrections for such conditions.
Besides previous research at the deforested site of this study (Wiekenkamp et al. 2016(Wiekenkamp et al. , 2019Ney et al. 2019), observations at a single height focusing on carbon budgets have been conducted at wind-thrown sites, either for the growing season following a storm (Lindroth et al. 2009) or long term (Lindauer et al. 2014Matthews et al. 2017). Vickers and Mahrt (2006) investigated mean vertical motions above a forest clearing, indicating horizontal divergence. Peltola et al. (2015) analysed the spatial representativeness of CH 4 fluxes over extensive and homogeneous grassland, while Nicolini et al. (2017) measured fluxes at two heights above a fetch limited crop field. However, the influence of measurement height on fluxes over heterogeneous surfaces with limited fetch, such as forest clearings, is still not sufficiently known.
For this study, we added an eddy-covariance system at a second height of 5.4 m to an already existing one at 2.7 m above a deforested site, which has a relatively inhomogeneous surface created by undisturbed vegetation growth. We hypothesize that (i) flux data from the old and new measurement heights can seamlessly be used together to estimate the cumulative carbon uptake and evapotranspiration from the clear-cut area, and (ii) the choice of LPF correction methods available in EddyPro© significantly affects flux results. In this respect, we also evaluated the performance of each method for data with high expected flux loss.

Site Description
Measurements took place at the upper Wüstebach catchment, located in the Eifel National Park near the Belgian border (50.50305 N, 6.33596 E, 618 m elevation; see Fig. 1). The site is part of the TERENO (TERrestrial ENvironmental Observatories) Eifel/Lower Rhine Valley observatory, which is one of four observatories in Germany to analyse long-term impacts of climate and land-use changes (Zacharias et al. 2011). Mean annual precipitation is 1332 mm and mean annual temperature 7°C during the reference period 1981-2010 (Ney et al. 2019). Cambisols are the dominant soil type in the north-eastern part of the study area, whereas Gleysols and Histosols prevail in a boggy area in the southern part and near the stream . Elevation within the target area ranges between 596 m at the outflow of the Wüstebach stream in the north-west, and 628 m at the eastern edge, with an average slope of 4°.
The site mostly consisted of spruce monoculture (Picea abies and Picea sitchensis) until 2013, when an 8.6 ha area of it was cut to allow for natural succession towards a European beech forest. The only major exception to this were isolated alder stocks near the stream, which were not cut. The eddy-covariance station is located approximately in the centre of the clear-cut with the forest edge closest to the north and north-east, with a minimum distance of 72 m, and farthest to the west, with a maximum distance of about 292 m (see Fig. 1a). Only 3% of the original biomass remained on site (Baatz et al. 2015), mostly tree stumps, litter, and a few tree trunks. In 2020, the vegetation of the clear-cut area consisted of various grasses (e.g., Deschampsia spec., Molinia spec.), shrubs and bushes of different size (e.g., Cytisus scoparius, Calluna vulgaris, Epilobium angustifolium), and young trees (Sorbus aucuparia, Betula pubescens, Picea abies), some of which are typical pioneer species. In general, regrown vegetation inside the fence, which had been established against browsing damage, was denser than outside. After the 2020 growing season, young trees within the clear-cut had an average height of 1.60 ± 0.89 m. The spruce trees demarcating the forest edge had a uniform height of about 25 m and measured alder trees near the stream heights between 8.0 m and 18.3 m. These characteristics resulted in a very heterogeneous study site with vegetation of different height and composition and scattered coarse woody debris, which is expected to result in heterogeneous source and sink areas for CO 2 and energy fluxes.
The site heterogeneity can be further characterized by flow tilt angles, calculated as tan −1 (W /U ), where W is the vertical wind component and U is horizontal wind speed during neutral conditions (z/L < 0.1, where z is the measurement height and L is the Obukhov length) and U > 1 m s −1 . For flow tilt angles shown in Fig. 2a, we applied a yaw rotation on unrotated u and v wind components to include both horizontal wind components in U . Positive flow tilt angles prevailed from western wind directions and likely originated from the sloping terrain. On the other hand, negative flow tilt angles from the north and northeast possibly originated from the nearby forest edge. Figure 2b shows flow tilt angles after the application of a sector-wise planar fit rotation of wind components after Wilczak et al. (2001) for each 45°sector. Here, flow tilt angles were largely diminished, especially for the prevailing western wind directions. A significant influence of the alder trees is not evident, which might be attributed to the fact that the lower elevation next to the creek prevents the tree tops from protruding considerably above the canopy surface around the station and on the far side of the creek. Issues with the northern wind sector between about 325°and 025°a re indicated by large u * discrepancies between the upper and lower system from the north after planar fit rotation (Fig. 2c), indicating distortions from the nearby forest edge on the measurements at the upper height.

Eddy-Covariance Set-Up and Processing
Turbulent fluxes of latent heat (λE), sensible heat (H ), and CO 2 have been measured at the study site with an eddy-covariance system since 2013. In April 2020, a second eddycovariance system was established to replace the first one at its current height of 2.7 m, which in turn was moved to a new height of 5.4 m due to vegetation growth. The upper system consisted of a CSAT3 sonic anemometer (Campbell Scientific, Logan, Utah, USA) and a LI-7500 open-path gas analyser (LI-COR, Lincoln, Nebraska, USA). The lower system also consisted of a CSAT3, and a LI-7500RS open-path gas analyser, which features optical hardware improvements compared to the LI-7500. Both systems had an orientation of 224°, while the upper system had a sensor separation of 22 cm and the lower system of 19 cm to account for higher LPF. In addition, a net radiometer (NR01, Hukseflux Thermal Sensors, Delft, Netherlands) at 4.54 m and two heat-flux plates (HFP01, Hukseflux Thermal Sensors, Delft, Netherlands) at −8 cm were installed to provide 10-min averages of net radiation and soil heat flux. Measurements taken between 17 April and 30 September during the 2020 growing season were analysed for this work.
Raw data of wind components (u, v, w), sonic temperature (T s ), and H 2 O and CO 2 densities logged at 20 Hz were corrected and processed to 30-min fluxes using the software EddyPro© (v7.04, LI-COR, Lincoln, Nebraska, USA). A sectorial planar fit rotation for 45°sectors after Wilczak et al. (2001) was applied for tilt correction of an anemometer misalignment and to account for inclination of the ground (see Fig. 2). Time lags between the anemometer and gas analyser were compensated for with the covariance maximization with default method, which uses a default value if no covariance maximum can be attained within a time-lag window. A high-pass-filtering correction (Moncrieff et al. 2004) was applied to account for attenuation resulting from block averaging. As LPF is expected to have a stronger impact on the lower system, emphasis was put on LPF correction methods. Hence, all five methods implemented in EddyPro© were selected and compared. These are the corrections after Moncrieff et al. (1997), Massman (2000) and Massman (2001), Horst (1997), Ibrom et al. (2007), and Fratini et al. (2012) (hereafter referred to by the first authors' names). A short description of each method and their implementation in this study is given in the next section. Lastly, the density correction of Webb et al. (1980) was added to the fluxes and the 0-1-2 flagging policy after Mauder and Foken (2004) was applied. The latter includes spike removal, a steady state test, and a test on developed turbulence after Foken and Wichura (1996). Fluxes were further separated for daytime conditions based on sunrise and sunset times to exclude several potential problems at night, such as advection and drainage flows (Aubinet 2008). Besides that, u * filtering implemented in the REddyProc library (Wutzler et al. 2018) was applied to remove remaining low-turbulence data. Furthermore, data from the northern wind sector between 325°and 025°were fully excluded because of likely distortions from the nearby forest edge (see Fig. 2c). Finally, data were rejected for which the source area originated to less than 70% inside the target area (see Sect. 2.5). For further analysis, only such timestamps were considered, for which the respective flux had the highest quality (flag 0) and all further criteria were also met simultaneously at both systems.
Surface heat correction for the LI-7500 of the upper system after Burba et al. (2008) was not applied because the correction was intended for vertically adjusted sensors while the gas analysers had an inclination of 45°. Furthermore, errors from self-heating are expected to be significant during very cold conditions (< -10°C) whereas only data during the growing season were analysed here. Ney et al. (2019) previously compared annual sums of surface heat corrected and uncorrected net ecosystem exchange values at this site, and also opted for uncorrected fluxes.

Low-Pass-Filtering Correction
Out of the five methods applied here, the ones after Moncrieff and Massman are fully analytical, meaning that filtering is described as individual spectral transfer functions, which are deduced from a priori knowledge of the system's properties, such as sensor separation, the instruments' time responses and path lengths, atmospheric conditions, and site characteristics. Flux attenuation is then estimated using a cospectral model, i.e., after Kaimal et al. (1972), as a reference of ideal cospectra. The method after Horst is also based on an analytical formulation but is parametrized using an in situ assessment of the system's cut-off frequency with measured spectra. The methods after Ibrom and Fratini rely on an empirical determination of the cut-off frequency from the ratio of ensemble gas spectra to ensemble normalized temperature spectra as a proxy of ideal gas spectra. For the Ibrom method, the correction factor is then parametrized using the cut-off frequency, average wind speed, and atmospheric stratification. The Fratini method uses this parametrization in a slightly different way only for small fluxes. For large fluxes, the correction factor is calculated using the cut-off frequency and sensible heat cospectrum instead. The resulting correction factor of each method is then multiplied by the uncorrected flux to correct for spectral attenuation.
For the corrections after Horst, Ibrom and Fratini, binned (co)spectra were calculated for every 30 min using EddyPro©. They were filtered according to the statistical tests after Vickers and Mahrt (1997), the micrometeorological quality tests after Mauder and Foken (2004), and by friction velocity (u * > 0.2 m s −1 ) and flux magnitude (H > 20 W m −2 ; λE > 20 W m −2 ; CO 2 flux > 2 μmol m −2 s −1 ) before they were ensemble averaged for unstable stratifications (−650 m < L < 0; predefined in EddyPro©). For all correction methods and gas spectra, the frequency range for fitting the in situ transfer function was set from 0.005 to 2 Hz. As the methods after Ibrom and Fratini compensate for LPF on the scalar signal only and thus do not account for sensor separation, the correction after Horst and Lenschow (2009) was applied in addition. However, it was only applied for crosswind and vertical wind components, as along-wind sensor separation was already compensated by the time delay correction using covariance maximization. In addition, H cospectra were preliminarily corrected for small losses due to anemometer path averaging and time response before using them for the Fratini method.

Energy Balance
The energy balance closure was calculated for each 30-min interval to estimate the performance of different spectral corrections and to compare the flux results of the upper and lower system. Ideally, the sum of sensible and latent heat fluxes measured by an eddy-covariance system should equal the available energy, that is net radiation minus ground heat flux and energy stored in the air and biomass (Wilson et al. 2002). Hence, the equation is where R n is the net radiation, G is the ground heat flux, P is the energy used for photosynthesis, and S the change of energy stored in the air and vegetation below the eddy-covariance measurements. As measurements were taken relatively close to the surface, S was neglected for this study. The energy P was also not measured but can be considered small compared to the other terms . The terms H , λE, R n , were measured directly, and G was assessed by correcting the soil heat flux for the estimated change in heat stored between the soil surface and the heat-flux plate according to Graf et al. (2020). Two energy balance parameters were calculated on a 30-min basis: i) the energy balance ratio (EBR) as the sum of turbulent fluxes divided by the available energy and ii) the energy balance closure (EBC) as the regression between the sum of turbulent fluxes and available energy. Here, a reduced major axis regression was used instead of an ordinary regression. In this way it is possible to handle likely random errors of available energy by evaluating the slope as the geometric mean of an ordinary regression and one with switched dependent and independent variables (Wilson et al. 2002). For net radiation, the maximum error was estimated at about 25 W m −2 ).

Footprint Estimation
Prior to footprint determination, roughness length z 0 and displacement height d were first estimated for each wind direction quadrant using wind velocities from the two systems during neutral conditions (z/L < 0.1). This was done by solving the logarithmic law after z 0 and d according to the Integrated Surface Flux System Guide (UCAR/NCAR 1990) where k is the von Kármán constant (0.41), U 1 and U 2 the wind speeds at measurement height z 1 and z 2 , and u * the mean of both heights. Values for d ranged from 1.07 m (north-east), 0.61 m (south-east), 0.53 m (south-west) to 0.51 m (north-west), and for z 0 from 0.29 m (north-east), 0.16 m (south-east), 0.24 m (south-west) to 0.18 m (north-west). The relatively small values for d compared to z 0 generally match the patchy structure of the study area well.
The two-dimensional footprint model of Kljun et al. (2015) was then applied to estimate the footprint for both heights for every 30 min, as well as a footprint climatology over the whole timeframe. The planetary boundary-layer height was hereby derived according to Appendix B in the respective paper. To analyse modelled source area differences between the two systems, footprint rasters for individual 30 min at the lower system were subtracted from the ones at the upper system in order to achieve patterns displaying if a pixel was more important for the upper or for the lower system. Two datasets were then created, for which negative (positive) pixels were set to 0 and remaining values converted to absolute values to show only pixels that were more important for the upper (lower) system. Lastly, the datasets were averaged over all timesteps to two raster images.

Spectral Analysis
Ensemble cospectra for the upper and lower system were calculated after the time delay correction was applied and are displayed in Fig. 3 for unstable conditions. Figure 3a and 3c show a clear attenuation of w C and w q cospectra at high frequencies (where C and q are the CO 2 and HO 2 mixing ratios) compared to w T s cospectra at both systems. The w T s cospectra can be used as reference cospectra because the sonic temperature is considered as an unfiltered scalar although it is also affected to a small extent by LPF due to path averaging and limited time response of the anemometer (e.g., Ibrom et al. 2007). Hence, the ratio of gas cospectra divided by w T s cospectra gives an experimental transfer function describing the spectral loss of CO 2 and H 2 O fluxes. Both cospectra, w C and w q , start diverging from w T s cospectra already at lower frequencies for the lower system compared to the upper system (Fig. 3b, d), resulting in a larger frequency loss and a higher demand for correction. The integral of w C cospectra in the inertial subrange (vertical lines in Fig. 3a, c) is 67% of the w T s cospectra integral for the upper system. For the lower system, this share is only 60%. Furthermore, ensemble sonic temperature spectra have a maximum density at 0.014 Hz for the upper system and at 0.02 Hz for the lower system, demonstrating a shift to higher frequencies for turbulent fluxes at the lower measurement height and thus a higher susceptibility to LPF. Likewise, the infinite impulse response filter cut-off frequency after Ibrom et al. (2007) is 1.1 Hz for the upper system and 1.0 Hz for the lower system.

Correction Factors
Correction factors of the five applied spectral correction methods for CO 2 fluxes of the upper and lower system are shown in Fig. 4. The correction factors across different spectral correction methods were similar with slightly higher values for the Fratini method. Correction factors for CO 2 fluxes were on average smaller for the upper system than for the lower system, with 1.06 (Moncrieff), 1.05 (Massman), 1.06 (Horst), 1.05 (Ibrom), and 1.08 (Fratini) for the upper system compared to 1.11 (Moncrieff), 1.07 (Massman), 1.10 (Horst), 1.09 (Ibrom), and 1.15 (Fratini) for the lower system. Correction factors for H 2 O were almost identical to those for CO 2 and are thus not displayed separately. The outliers of high correction factors from the Moncrieff and Massman methods are associated with low U , for which the high-pass-filtering correction increased correction factors. For the Fratini method, not constraining correction factors to the bounds of Eq. 9 in Ibrom et al. (2007) but depending on the stochastic nature of turbulence by incorporating H cospectra may have led to the outliers. Furthermore, U is considered an important factor for spectral attenuation as high wind speeds favour high-frequency eddies and therefore correction factors might be expected to increase with U (Moncrieff et al. 1997). The dependency of the CO 2 correction factors on U is shown in Fig. 5 for the different LPF methods and for both systems. The correction factors after Moncrieff, Horst, Ibrom, and Fratini slightly raised with increasing U , more pronounced for the lower system than for the upper one. The correction factors after Massman, however, do not show any dependence on high U at all.

Flux Results
Out of 4729 possible daytime 30-min intervals, 1992 (42%) CO 2 flux data points passed the quality tests at both heights simultaneously (λE: 47%; H : 47%). A substantial part of missing data does not result from quality control but is caused by a data acquisition failure of the upper system during 20 days in July. Table 1 lists mean values of LPF-corrected and uncorrected fluxes for those selected 30-min intervals. It can be noted that higher H and especially λE were measured at the upper system compared to the lower one. CO 2 , however, were slightly more negative at the lower system, especially after corrections were applied. As expected, uncorrected H values differ only slightly from corrected ones, whereas the CO 2 flux and λE differ more strongly, especially at the lower system. A reduced major axis regression between 30-min flux values of the upper and lower system corrected after Moncrieff gives slopes of 1.14 (R 2 = 0.71), 1.23 (R 2 = 0.93) and 1.03 (R 2 = 0.97) the CO 2 flux, λE, and H respectively. For CO 2 fluxes, the two least correlated LPF corrections were Massman against Fratini at the lower system with a slope of 0.86 and an R 2 of 0.95. The average Bowen ratio is 0.83 for the upper system and 0.98 for the lower system. The uncertainty induced from the choice of LPF corrections was calculated as the standard deviation between fluxes of the five LPF corrections averaged over all 30-min intervals. For CO 2 fluxes of the upper system, the average flux with uncertainty of the five LPF corrections was -4.75 ± 0.16 μmol m −2 s −1 (3.4%), for the lower system it was -5.63 ± 0.35 μmol m −2 s −1 (7.0%). For λE, these values were 111.33 ± 1.51 W m −2 (1.4%) for the upper system and 94.50 ± 2.91 W m −2 (3.0%) for the lower system (the percentages represent the size of this uncertainty compared to the flux, averaged over all time steps). Likewise, the uncertainty from the two measurement heights was on average 24.8% of CO 2 fluxes and 9.7% of λE, averaged over all LPF corrections. As a comparison between the two most contrasting LPF corrections, CO 2 (λE) fluxes corrected after Fratini were on average 7.4% (2.9%) higher than after Massman at the upper system and 16.2% (8.1%) higher at the lower system. In contrast, CO 2 fluxes averaged over all LPF corrections were 18.3% higher at the lower system, whereas λE was 17.8% higher at the upper system. For comparison, random errors were estimated according to Finkelstein and Sims (2001) on a 30-min basis and were on average 34.2% (26.3%) of CO 2 fluxes of the upper (lower) system and 16.2% (10.6%) of λE. Figure 6 shows the frequency of the CO 2 flux and λE corrected after Moncrieff separated for 10°wind direction sectors. Wind generally prevailed from western directions, whereas wind from southern and especially northern directions was less frequent. For western wind directions, a slight clockwise wind direction shift from the lower to the upper system is also noticeable. At both heights, the highest average CO 2 uptake was detected from wind directions between 225°and 270°, with a CO 2 flux of -5.85 (-6.43) μmol m −2 s −1 at the upper (lower) system. The smallest CO 2 uptake was recorded from 090-135°(-3.30 μmol m −2 s −1 ) at the upper system and from 180-225°(-3.59 μmol m −2 s −1 ) at the lower system. The highest CO 2 flux differences between the two systems were recorded from 045-090°, both on average (1.14 μmol m −2 s −1 ) and summed up (8.17 g C m −2 ). For λE, the highest average flux was recorded from 135-180°(131.0 W m −2 ) at the upper system, and from 180-225°(107.3 W m −2 ) at the lower system. The highest λE differences between the two systems occurred on average from 135-180°(131.0 W m −2 ), and summed up from 225-270°( 16.36 mm), as it was a more frequent wind direction.
Fluxes of both daytime and night-time conditions are visualized as cumulative fluxes in Fig. 7 for the corrections after Massman and Fratini, as examples of an analytical and an in situ method, as well as without correction. Cumulated H 2 O fluxes were 23% higher at the upper system than at the lower one without LPF correction, but this discrepancy was smaller after corrections were applied (14% higher after Fratini). For cumulative CO 2 sequestration, in contrast, the lower system yielded larger cumulative fluxes (12% without correction), and corrections even increased this discrepancy (28% after Fratini). The correction after Fratini produced larger fluxes than the correction after Massman of both CO 2 (upper system: 8%, lower system: 19%) and H 2 O (upper system: 3%, lower system: 8%). Relations between the upper and lower system and between correction schemes remained consistent over the vegetation period.

Energy Balance Closure
An ideal EBC would be represented by a slope of 1 and an intercept of 0 from a linear regression of the sum of turbulent fluxes versus available energy, and an ideal EBR would be 1. Table 2 demonstrates that EBC and EBR were generally better at the upper system than at the lower system with a higher EBR of about 0.09, though at both heights the sum of turbulent energy fluxes was lower than the available energy. All correction methods improved the EBC as well as the EBR and had relatively similar results, with Fratini performing slightly better especially for the lower system (EBR of 0.81). However, R 2 of the reduced major axis regression was slightly higher at the lower system than at the upper one.

Footprint
The footprint climatology estimations are displayed in Fig. 8 and show footprints for the upper system extending about 2.5 times as far from the tower than for the lower system. The Fig. 7 Cumulative non gap-filled evapotranspiration (a) and CO 2 fluxes (b) including both daytime and nighttime situations that passed quality control. Up and Lo stand for the upper and lower system, respectively 90% cumulative footprint contour line contains forest outside the target area, whereas the 90% line of the lower system is still within the target area. Only 0.8% of all 30-min intervals of the lower system originated to less than 70% within the target area, whereas for the upper system that value was 4.2%. When considering daytime conditions only, these values dropped to 0.1% for the lower system and 0.2% for the upper system. A 70% threshold was used for flux filtering to be in line with a previous study at the research site (Ney et al. 2019), discarding all affected values from the further analysis. The general shape of the footprints extends in the east-west direction, with long upwind distances resulting from prevailing winds from western directions.
Wind direction changes between the upper and lower system are expressed in the averaged footprint differences between the two systems (Fig. 9). The upper system recorded more wind from north-western and north-eastern directions, while the lower system recorded more wind from western and south-eastern directions. The relatively large wind frequency differences between the heights from north-western wind directions, however, did not yield to analogous λE or CO 2 flux differences from these directions. The modelled results further indicate that source area differences between the two systems mostly originated from within the fence not farther than about 130 m from the tower, while the area in the direct vicinity (< 20 m) around the tower was more important for the lower system.

Assessment of Measurement Height
The higher daytime CO 2 uptake observed at the lower system despite lower evapotranspiration is a counterintuitive result that requires clarification. One possible explanation is that within the source area of the upper system more evaporation without accompanied photosynthesis occurs. Such an area could be the wetter and temporarily flooded Gleysols and Histosols in the Fig. 9 Relative importance of pixels more important for the upper system (a) and pixels more important for the lower system (b), without observations from north. Please note the different colour scales for the subplots, which were expedient to visualize the relatively smaller importance of individual pixels in (a) compared to (b). The coordinate origin is at the eddy-covariance station and the white line delineates the fence for orientation southern part of the clear-cut and near the stream to the west. The soil moisture measurements of Wiekenkamp et al. (2016) show that soil moisture was about 30% higher there compared to the direct vicinity of the tower in summer 2014. Furthermore, Graf et al. (2020) showed that peatlands responded with a disproportionately low ratio of CO 2 uptake to evapotranspiration compared to other ecosystems during drought conditions. For the study site, this was the case in the previous two years and may have affected flux in 2020. Accordingly, the highest mean λE and CO 2 flux were recorded from different wind directions, and the highest λE differences between the two systems were recorded from south-south-east on average, while cumulative differences were largest from west-south-west due to more observations from there. The ratio between CO 2 uptake and evapotranspiration was also the lowest from 135-180°compared to other wind directions at both systems. However, it remains not fully explained why the south-south-east sector stands out compared to south-western wind directions with similar or even wetter soil properties (Wiekenkamp et al. 2016). Figures 8 and 9 also indicate that source area differences between the two heights mostly originated from within the fence, thus not clearly demonstrating a strong influence of the boggy area outside of it. The lower system, on the other hand, could have a location bias of young, fast-growing trees located in the direct vicinity of the tower (see Fig. 1b), which could have increased CO 2 uptake there. The energy balance could not be closed with either LPF correction for both systems but was generally better for the upper system. A number of studies demonstrated that even with very carefully applied eddy-covariance set-ups, the sum of turbulent energy fluxes remained below the amount of available energy (Foken et al. 2010;Stoy et al. 2013). It is assumed that this results from low frequency eddies not detectable by eddy-covariance systems because of a limited averaging period (Foken 2008). In addition, closure of the energy balance cannot be expected for a heterogeneous exchange surface inducing advection . Hence, the lower sums of H + λE compared to the available energy are in line with expectations based on previous research. Since advection could be present even in rigorously filtered data, it could have both increased or decreased the EBR at both systems. Vickers and Mahrt (2006) showed that a mass continuity approach indicated longterm sinking motions above a forest clearing. Tilt corrections such as planar fit, however, remove the mean vertical motion, hence partially not taking into account vertical advection in the flux averaging period. Such long-term sinking motions above the study area may be induced by the rough-to-smooth surface change or by drainage flows following the sloped terrain (Lee 1998). On the other hand, Vickers and Mahrt (2006) also pointed out that vertical advection of CO 2 based on mass continuity was a large term of net ecosystem exchange mainly on weak mixing nights, which were excluded for this analysis altogether. Besides that, the EBR discrepancies between the upper and lower system can have multiple causes. G can differ within the footprint of each system from the measurements beneath the tower as soil properties are not uniform throughout the deforested area . However, it is expected that average G of the lower system's smaller footprint is more similar to the measured G than that of the upper system. The same applies to possible differences in net radiation, in particular due to different surface albedos. The energy stored in the air and biomass was investigated to be negligible for similar vegetation heights ). If anything, the error induced by disregarding the energy storage in the air should be higher for the upper system, where the air column beneath is larger compared to the radiometer. Hence, EBC might be poorer at the lower system partially because high frequency attenuation is not fully compensated by spectral corrections, as indicated by considerably lower λE at the lower system in Table 1. This raises the question of whether the CO 2 flux at the lower system is likewise underestimated, given the spectral similarity between CO 2 and H 2 O fluxes. Thus, insufficient LPF corrections at the lower system may also have contributed to the counterintuitive flux results described above.
Spectral corrections might be insufficient for the lower system because sensors are not placed high enough in the inertial sublayer. Moore (1986) stated that for his analytical correction the measurement height above d should be at least 10 times the sensor separation. For the lower system this means 2.58 m above ground, which was barely met in our case. Measurements in the roughness layer can yield the CO 2 flux and λE representing only local disturbances and thus being spatially variable within the same ecosystem (Katul et al. 1999).
However, a precise definition of the roughness layer height and thus an appropriate measurement height is still lacking. For structurally complex ecosystems, Munger et al. (2012) where h c is the average canopy height. Since rough estimates of the average h c are generally accepted (Rebmann et al. 2018), we calculated it as d/0.67 according to the EddyPro© manual, which results in z ≈ 2.48 m above ground. For shrublands, however, Munger et al. (2012) recommended a fixed height of about 6 m, which in our case was barely satisfied by the upper system. Nicolini et al. (2017) accomplished feasible measurements as low as 0.9 m above d, though over a homogeneous surface. Although these recommended heights can only be seen as very rough estimates, they indicate that the lower system might be at best at the lower end of the suitable range of z and will be in the roughness layer with further expected vegetation growth.
The footprint estimates revealed that the source area of the upper system was to a large extent within the target area and therefore only few observations were removed. However, during 2014-2017 the 90% cumulative footprint of the lower system had a maximum distance of about 200 m from the tower (Ney et al. 2019). In 2020, this distance decreased to 123 m, while the 90% footprint of the upper system had a maximum distance of about 311 m. This result indicates that the source area of both systems differed from previous observations, but with further vegetation growth it is expected that the upper system's source area will approximate that of previous measurements by the lower system. The footprint model of Kljun et al. (2015) assumes horizontal homogeneity of the flow and thus has limited applicability to the study area. The complex flow over the forest edge particularly cannot be resolved, for which large-eddy simulations or, as a less computationally intensive solution, turbulence closure models such as SCADIS would be more suited (Sogachev and Lloyd 2004). This model was also able to indicate source hotspots in contrast to analytical footprint models in heterogeneous areas (Sogachev and Dellwik 2017). A recirculation area behind the edge inducing downward flows can be expected for a distance of 2-5 times the forest canopy height (Detto et al. 2008), corresponding to a distance between 50 and 125 m at the study site. This is a problem for northern wind directions where the forest edge is within this distance, and distortion of the mean flow is indicated by much higher u * at the upper system (see Fig. 2c). On the other hand, for the prevailing western wind directions, such edge turbulence effects were not detected. Roughness changes were also roughly taken into account for footprint modelling by including z 0 for each wind direction quadrant. Hence, the footprint results might be useful for a first approximation of the source area and for testing spatial representativeness of the fluxes.
Despite these general considerations, the presented results strongly speak against the first hypothesis. The large differences between fluxes of the two heights (see Table 1 and Fig. 7) prevent a seamless use of data from both time series and likely result from a different source area within the heterogeneous clear-cut area and insufficient LPF corrections for the lower system. Instead, in any future analysis of CO 2 fluxes at the clear-cut, the period with two simultaneous measurements heights can be used to estimate the uncertainty from measurement height choice, which can then be compared to long-term trends or differences between sites.

Spectral Corrections
The higher average correction factors for the lower system throughout all methods are in line with the higher spectral attenuation observed there compared to the upper system (see Fig. 3c, d). This observed shift to higher frequencies with a lower sensor height coincides well with other experiments and well-known theoretical considerations (e.g., Moncrieff et al. 1997;Foken et al. 2012;Zhao et al. 2019). The correction factors after Moncrieff and Massman show a clear dependence on U because specific quantities of transfer functions are defined as functions of U there. However, correction factors actually decrease initially with increasing U , since attenuation dominates in the low frequency range due to block averaging at U < 0.5 m s −1 but becomes less important with increasing U in unstable conditions. An insensitivity of correction factors to U can be observed for the Massman method because for open-path systems, time constant equivalents from path averaging and sensor separation decrease with increasing U , and thus were assumed to compensate the shift to high frequency eddies (Massman 2000). These comparatively small correction factors at higher wind speeds resulted in slightly smaller λE and CO 2 flux at both systems for the Massman method. Polonik et al. (2019) concluded that the Fratini correction is not well-suited for open-path analysers because it accounts only for scalar attenuation, as it does not consider sensor separation, and therefore produced smaller fluxes than Massman. However, with the additional correction after Horst and Lenschow (2009), this limitation was not an issue for our analysis. The fluxes of the lower system corrected after Fratini had a higher magnitude and better energy balance closure compared to Massman or other methods, confirming its applicability to low measurement heights, even for a comparatively rough surface. Polonik et al. (2019) did not apply the correction after Horst and Lenschow (2009) because it produced unrealistically high correction factors in stable conditions, but in our case this correction increased the correction factor only by 0.07 for the upper system and 0.13 for the lower system during stable conditions. Nonetheless, in a few cases it added large values up to 0.7 to the correction factor. In unstable conditions, the maximum value added was 0.07 for the upper system and 0.24 for the lower system. Fratini and Mauder (2014) found a difference of about 3% in λE and the CO 2 flux caused by the use of spectral corrections either after Moore (1986) or Horst (1997), which contributed most to discrepancies between flux processing in EddyPro© and TK3. In our analysis, the highest differences (16.2%) were found between CO 2 fluxes at the lower system corrected after Fratini and Massman. In contrast, Rannik et al. (2020) assessed that differences in fluxes from the choice of coordinate rotation were less than 10%. Nevertheless, it should be kept in mind when comparing different spectral correction schemes that spectral corrections are not the last step in the processing chain of EddyPro©, but density-correction terms (Webb et al. 1980) are further added, which in addition can be implemented differently in other software (Fratini and Mauder 2014). It is also important to assess the importance of sources of uncertainty, such as measurement height and LPF correction, against the magnitude of real fluxes between sites or years that are the target of past and future studies on carbon budgets of forests and clear-cuts. For example, for annual net ecosystem exchange over the first four years after deforestation, Ney et al. (2019) found a source-towards-neutral change of 0.439 kg C m −2 , and differences of more than 0.600 kg C m −2 compared to the surrounding spruce forest. The largest differences of daytime cumulative growing season C uptake resulting from combinations of LPF correction and measurement height in our study, 0.650 kg C m −2 between Massman of the upper system and Fratini of the lower system, would not change these results fundamentally but account for a non-negligible additional relative uncertainty.
The uncertainty resulting from the choice of LPF correction can be subsumed under systematic errors associated with data processing in the classification scheme of Mauder et al. (2013). Other sources of uncertainty include systematic errors from instrumental calibration and random errors due to changes in footprint, instrumental noise, or the stochastic nature of turbulence. Stochastic errors estimated according to Finkelstein and Sims (2001) were considerably larger than the differences induced by the choice of a LPF correction method on a 30-min basis. Over longer time periods, however, random errors are cancelled out, whereas systematic differences from LPF corrections add up, as illustrated in Fig. 7. These discrepancies of different LPF corrections were stronger at the lower system, where LPF and concurrent correction factors were higher than at the upper system. Therefore, the results support the hypothesized importance of the choice of LPF correction, although flux differences between the two measurement heights were larger than even between the two most contrasting LPF corrections.

Conclusion
We compared turbulent flux measurements at two heights above a clear-cut site, demonstrating the trade-offs that have to be considered when choosing the measurement height above a fetch-limited heterogeneous surface. Major limitations of these results include potential advection biasing the EBR at both heights and the limited applicability of the Kljun et al. (2015) footprint model to a study site with heterogeneous flow. However, the footprint model has shown a limited utility for estimating the influence of source and sink heterogeneities within the clearing. The upper system, with its larger footprint, is more influenced by the forest and edge turbulence effects from the northern sector, while the lower system likely lacks representativeness of the clearing and is susceptible to higher LPF. These effects resulted in significant flux discrepancies between the two heights, which oppose the first hypothesis that a seamless use of the data from both time series is acceptable. We also evaluated different LPF correction schemes. The differences between the methods after Moncrieff, Massman, Horst, Ibrom, and Fratini induced a systematic uncertainty to the fluxes, which was stronger for the lower system (CO 2 : 7.0%, H 2 O: 3.0%) than for the upper system (CO 2 : 3.4%, H 2 O: 1.4%). The flux discrepancies of the different correction methods added up over time and hence support the second hypothesis. Compared to other methods, the Fratini approach yielded higher fluxes and a better energy balance closure for the lower system. Hence, our analysis confirms that for long-term single-point flux observations above forest clearings, information about changes of measurement height are critical for interpreting the data, and that it is also important to consider the spectral correction method when comparing fluxes between sites.