Evaluating multi-year, multi-site data on the energy balance closure of eddy-covariance flux measurements at cropland sites in southwest Germany

The energy balance of eddy covariance (EC) measurements is typically not closed resulting in one of the main challenges in evaluating and interpreting EC flux data. Energy balance closure (EBC) is crucial for validating and improving regional and global climate models. To investigate the nature of the gap in EBC for agro-ecosystems, we analysed EC 20 measurements from two climatically contrasting regions (Kraichgau (KR) and Swabian Jura (SJ)) in southwest Germany. Data were taken at six fully equipped EC sites from 2010 to 2017. The gap in EBC was quantified by ordinary linear regression, relating the energy balance ratio (EBR) calculated as the quotient of turbulent fluxes and available energy, to the residual energy term. In order to examine potential reasons for differences in EBC, we compared the EBC under varying environmental conditions and investigated a wide range of possible controls. Overall, the variation in EBC was found to be higher during 25 winter than summer. Moreover, we determined that site had a statistically significant effect on EBC, but neither did crop nor region (KR vs SJ). The time-variable footprints of all EC stations were estimated based on data measured in 2015, complemented by micro-topographic analyses along the prevailing wind direction. The smallest mean annual energy balance gap was 17 % in KR and 13 % in SJ. Highest EBRs were mostly found for winds from the prevailing wind direction. The spread of EBR distinctly narrowed under unstable atmospheric conditions, strong buoyancy, and high friction velocities. 30 Smaller footprint areas led to better EBC due to increasing homogeneity. Flow distortions caused by the back head of the anemometer, negatively affected EBC during corresponding wind conditions.


Introduction
Studying turbulent exchange at the land surface is important for assessing water cycling, plant growth and carbon fluxes of ecosystems and for enhancing soil-crop, climate and weather models. Currently, the best technique to determine these fluxes is the eddy covariance (EC) method. It is considered the most direct and accurate measurement of turbulent fluxes in the soilplant-atmosphere system (Baldocchi et al., 2001;Burba, 2013). In EC flux data, the measured available energy (incoming net 5 radiation minus ground heat flux) is generally higher than the sum of turbulent exchange fluxes (latent and sensible heat).
Accordingly, either the turbulent fluxes are incompletely captured or the measured available energy is positively biased. This gap in energy balance closure (EBC) is a long-standing problem in EC measurements and one of the most frequently discussed concerns in micrometeorological research (Foken, 2008a).
Globally, a large number of research sites has been established to, inter alia, study reasons for the energy imbalance. This 10 includes the FLUXNET network with more than 500 EC towers around the world (Wilson et al., 2002), and the AmeriFlux network operating in North, Central and South America (Peng et al., 2017). An example is an Energy Balance Experiment which was conducted to determine the reasons for the energy imbalance of EC measurements over irrigated cotton fields. The results showed that the net radiation differed by up to 10 W m −2 across a single field . In a review on EBC, Foken (2008b) summarized the most important factors for the energy imbalance as related to; measurement errors of the 15 energy balance components, incorrect sensor configurations, influences of heterogeneous canopy height, unconsidered energy storage terms in the soil-plant-atmosphere system, inadequate time averaging intervals, and long-wave eddies (mesoscale circulations) (Foken, 2008b;Jacobs et al., 2008;Wilson et al., 2002). Additionally, the energy transport with near-surface secondary circulations (large eddies) cannot be measured with a single EC station (Cava et al., 2008;Foken, 2008b;Xu et al., 2017). 20 In parts, this may be rectified by altering the time averaging intervals. For example, Kidston et al. (2010) found that at a forest site, EBC peaked at 90 % when applying a 240 min averaging interval. At a boreal forest site, Sánchez et al. (2010) applied a range of different time averaging intervals and found that increasing the interval from the traditional 30 min to one day improved the EBC from 75 to 100 %. However, this picture is inconclusive, since Oncley et al. (2007) found at an irrigated cotton field that increasing the time averaging interval up to 4 h did not result in a higher contribution of turbulent fluxes, and 25 that the contribution of low turbulent fluxes was less than 10 W m −2 . In most cases, the standard 30-min averaging period is proven to be the best compromise for simultaneously capturing most of the turbulent fluxes while fulfilling the precondition of stationarity (Charuchittipan et al., 2014;Masseroni et al., 2014;Sun et al., 2006).
The influence of site characteristics (e.g. vegetation type, canopy height, terrain) on the EBC has been studied extensively. Wilson et al. (2002) reported no clear differences in EBC between flat and sloped terrain sites across 22 research sites. A 30 comparison of two different agroecosystems in China, a degraded grassland and a maize cropland, showed similar EBC of about 80 % (Du et al., 2014). The comparison of a mature boreal jack pine forest and a jack pine clear-cut site by Kidston et al. (2010) revealed that depending on the surface characteristics the loss of low frequencies can contribute significantly to the energy imbalance. Canopy height may impact EBC although Wilson et al. (2002) found in their study that vegetation height did not control EBC. However, consideration of the stored energy in the soil-plant-atmosphere zone can noticeably improve EBC (Jacobs et al., 2008;Meyers and Hollinger, 2004;Zeri and Sá, 2010). Meyers and Hollinger (2004) compared the energy stored or released by CO2 exchange and crop enthalpy change and showed that in their study maize stored more energy than soybean crops. Additionally, Eshonkulov et al. (2019) reported that mean EBC was improved from 78 % to 87 % when minor 5 energy storage and flux terms were taken into account during the main vegetation period. Lastly, the mismatch of measurement scales is also considered to be a reason for the energy imbalance (Sánchez et al., 2010;Xu et al., 2017).
During the last decade, the identification of the contributing source area (footprint) and evaluation of the representativeness of the EC flux data for the field of interest has found an increased interest (Göckede et al., 2006;Kljun et al., 2004;Schmid, 2002). Knowledge about the footprint is important to clarify whether the EC station measures local or nonlocal energy fluxes 10 (Eugster and Merbold, 2015;Pirk et al., 2017). Currently, there is a variety of models in use to estimate footprint areas. While most analytical footprint models assume a homogeneous flux source area, footprint calculations for heterogeneous sites require greater computational effort and detailed information on surface characteristics (Mauder et al., 2013). Despite the existing methods and studies, Stoy et al. (2013) concluded that the relationship between footprint and EBC in agricultural cropland has not been sufficiently studied. 15 Therefore, the presented study evaluates the energy balance at agricultural crop lands. The analyses are based on EC measurements conducted from 2010 to 2017 at six fully equipped sites in two climatically different regions of southwest Germany. We hypothesized that multi-year, multi-site observations will provide new insights into the nature of the energy imbalance of EC flux measurements. The objectives of this study are to evaluate if crop type, site characteristics, wind direction, atmospheric conditions, and footprint area act as controls on the EBC. 20

Site description
The study sites in the Kraichgau region (KR) were located at "Katharinentalerhof" characterised by mostly flat terrain and are located approximately 4 km north of the city Pforzheim (48.92° N, 8.70° E). Three EC stations (EC1, EC2, and EC3) were 25 installed at adjacent fields with the respective areas of 14.9, 23.6 and 15.8 ha (Fig. 1). The prevailing wind direction is west.
Approximately 500 m to the South of the experimental fields a former landfill site is located whose maximum elevation is about 41 m above its surroundings. KR is one of the warmest regions in Germany with a mean temperature and annual precipitation of 9.4 °C and 890 mm during 1981−2010 (Meteorological station Pforzheim-Ispringen, German Weather Service, located about 3 km from the research sites). The soils of this region developed from deep loess layers overlying a shell 30 limestone. Detailed information about meteorological and soil conditions can be found in in Table 1 and in Imukova et al. (2016), Ingwersen et al. (2015), and Wizemann et al. (2014).
Due to its higher elevation, the Swabian Jura region (SJ) is characterized by a colder and harsher climate compared to KR. The prevailing wind direction is southwest to west. Mean temperature and annual precipitation were 7.5 °C and 1042 mm during 1981−2010 (Meteorological station Merklingen, German Weather Service, about 2 km from the research sites). Information 5 about meteorological and soil variables is given in Table 1. Accordingly, crops are generally sowed and harvested later than in KR. SJ is the largest contiguous karst landscape in Germany, with generally rather shallow soils. The study sites are located close to the town of Merklingen (Fig. 1). The areas of the three research fields at EC4, EC5, and EC6 were 8.7, 16.7 and 13.4 ha, respectively. While EC4 and EC5 were adjacent fields, EC6 was situated 1.5 km to the north (Fig.1).
The crop rotation in SJ was more diverse than in KR and the most frequently grown crops were winter wheat and silage maize 10 (Table 2). In 50 % of the site years in KR, winter wheat was cultivated; this value was only 25 % in SJ. KR showed a lower variety in cultivated crops to SJ, with three in KR and six in SJ. At all sites, farmers frequently grew cover crops between winter and summer crops. These were mainly mustard, phacelia or multi-species mixtures.

Eddy covariance measurements 15
One EC station was installed at the center of each field site in spring 2009 (Ingwersen et al., 2011;Wizemann et al., 2014).
All stations were equipped with a fast-response CO2/H2O open-path infrared gas analyzer (LI-7500; LI-COR Biosciences, Lincoln, NE, USA) and a three-axis ultrasonic anemometer (CSAT3; Campbell Scientific Inc., Logan, UT, USA). The raw data of the gas analyzer and sonic anemometer were recorded at 10 Hz and stored on a CR3000 data logger (Campbell Scientific Inc., Logan, UT, USA). In early 2009, the CSAT3 orientation at EC1 and EC3 was 230°, and at EC2 255°. In late April 2010, 20 the orientation was changed to 170° and varied over the subsequent years between 160° and 190°, ensuring that winds from the prevailing wind direction (west and east) enter the anemometer from the side. In SJ, the mean CSAT3 orientation was 220°±15 from late March 2010 until the end of 2017. The gas analyzers were factory-calibrated biannually. Sensor heights were adjusted accounting for increasing canopy heights, particularly during the vegetation periods of maize. This ensured that the distance between sensors and canopy was roughly 2-3 m. Maximal sensor heights in KR and SJ were 6.00 m at EC2 (2014) 25 and 4.80 m at EC6 (2010), and the minimal sensor height was approximately 2 m in both regions. Each EC system was powered by two 12 V batteries (each 240 Ah) charged by four 20 W solar panels. To enable continuous EC measurements during winter, portable fuel cell systems (Efoy Pro 800 Duo, FSC Energy AG, Brunnthal-Nord, Germany) were installed in autumn 2015, one at EC2 and one at EC6. At the others stations the LI-7500 was shut down during the winter, mostly from late November to mid-March. 30 Net radiation was measured using a 4-component radiometer (NR01, Hukseflux Thermal Sensors, Delft, The Netherlands).
The radiometers were placed above the cropped field area in close proximity to the EC stations. Air temperature and relative humidity were measured at a height of 2 m at each EC station using a temperature and relative humidity probe (HMP45C, Vaisala Inc, Helsinki, Finland). Soil temperature was measured at the depths of 0.02, 0.06, 0. 15,0.30 and 0.45 m (107 Thermistor probe,Campbell Scientific Inc.,Logan,UT,UK). To measure the soil heat flux near the EC stations, three heat flux plates (HFP01, Hukseflux Thermal sensors, Delft, The Netherlands) were installed at a depth of 0.08 m. The soil volumetric water content at 0.05, 0.15, 0.30, 0.45 and 0.75 m depth was monitored with frequency-domain reflectometry 5 sensors (CS616, Campbell Scientific Inc., Logan, UT, USA). In the shallow soil at EC6, however, soil variables could be measured only down to 0.3 m. Data from thermistor (0.02 m and 0.06 m) and water content sensors (0.05 m) were used to calculate the soil heat storage between the soil heat flux plates and the ground surface (Eshonkulov et al., 2019;Wizemann et al., 2014). Precipitation was measured with a 0.2 mm tipping bucket rain gauge (ARG 100, Environmental measurements Ltd., North Shields, UK) which was installed 1 m above ground. The rain gauges were recalibrated once per year. 10

Data processing and quality control
High-frequency raw data from 2010 to 2017 were processed with 30 min averaging interval using the software package TK3.1 (Mauder et al., 2013). The following settings were used to compute latent and sensible heat flux: spike detection (Vickers & Mahrt, 1997), planar fit coordinate rotation (Wilczak et al., 2001), correction of spectral loss (Moore, 1986), sonic virtual 15 temperature conversion into actual temperature (Schotanus et al., 1983) and correction for density fluctuations (Webb et al., 1980). Additionally, the raw data of 2015 were processed with the software Eddypro ® (Version 6.2.1, LI-COR Inc., 2012) to obtain input parameters (Obukhov length, standard deviation of lateral velocity fluctuations after rotation, friction velocity, mean wind speed and direction) for deriving flux source area (footprint). Data processing and correction in EddyPro ® were conducted with the same settings as in TK 3.1. Both programs yield comparable results (Fratini and Mauder, 2014). 20 The nine flag system after Foken et al. (2004) was used as the quality criterion. For the evaluation, we used only data with quality flags 1-3 as suggested by Mauder and Foken (2011). Moderate (flags 4-6) and poor quality (flags 7-9) data were discarded. In a second step, a median filter was applied for additional de-spiking of half-hourly fluxes. The filter removes all fluxes exceeding 5 times the median of the previous three days (Demyan et al., 2016). No gap-filling was performed in this study. 25

Energy balance closure of eddy covariance measurements
In the ideal case, the surface energy balance obeys the following equation: where Rn is the incoming net radiation, LE is the latent heat flux, H is the sensible heat flux (both positive upwards) and G is 30 the ground heat flux (positive downwards). All components are expressed in W m −2 . Note that in Eq. 1 minor flux terms such as energy storage in the canopy or energy conversion by photosynthesis are neglected. All available filtered half-hourly flux data of the four terms in eq. 1 were used to calculate the EBC. Three measures were used to evaluate the EBC. Firstly, we determined the slope and intercept from ordinary linear regression (OLR) of turbulent fluxes (H + LE) against available energy (Rn − G). In the ideal case of a fully closed energy balance, the slope and intercept of the linear regression are equal to one and zero, respectively (Ping et al., 2011;Wilson et al., 2002). In this study, we considered also the intercept (W m −2 ) of the OLR 5 in evaluating EBC as suggested by Franssen et al. (2010).
Secondly, by the energy balance ratio (EBR) calculated as: 10 Thirdly, by comparing the energy balance residual (Res; W m −2 ) given by:

Atmospheric conditions
As a proxy for the role of shear and buoyancy in the production or consumption of turbulent kinetic energy, we used the friction 15 velocity, u * (m s −1 ), and the kinematic virtual temperature flux, respectively. The latter is the covariance (w'Tv') between vertical wind speed (w) and virtual temperature (Tv). As the virtual temperature can be replaced by the sonic temperature (Ts) with negligible loss of accuracy (Kaimal and Gaynor, 1991), we computed the virtual temperature flux from the covariance (w'Ts') between w and Ts.
The relationship between atmospheric stability and the EBC was examined using the dimensionless atmospheric stability 20 parameter ζ defined by Stull (1988): where zm (m) is the measurement height of the sonic anemometer and L (m) is the Obukhov length. The stability parameter 25 expresses the relative roles of shear and buoyancy. Using ζ, the stability of the atmosphere can be divided into three classes (Franssen et al., 2010): stable (ζ ≥ 0.1), neutral (−0.1 < ζ < 0.1) and unstable (ζ ≤ -0.1).

Footprint analyses and micro-topography
To determine the relationship between the contributing source area of turbulent fluxes and the EBC, we performed footprint analyses. We used the flux footprint prediction (FFP) online tool of a simple two-dimensional parametrization (Kljun et al. 2015, http://geography.swansea.ac.uk/nkljun/ffp/www/). The footprint parametrization uses the Lagrangian stochastic particle dispersion model (Kljun et al., 2002). As input parameters to the model, we used displacement height, zd (m), mean wind speed 5 (m s −1 ), Obukhov length (m), standard deviation of horizontal wind speed (m s −1 ), friction velocity, u * (m s -1 ), wind direction (°), and measurement height above the ground surface, zm (m), which was calculated by 10 where is the height of the sonic anemometer and the gas analyser, and is calculated by where (m) is the time variable canopy height because of crop growth. This was accounted for by biweekly measurements, 15 and considered in TK3.1 for the respective two-week periods. Data for footprint analyses were constrained to u* > 0.1 m s −1 and ζ ≥ −15.5.
Additionally, the micro-topography of the EC sites was determined along a transect in the prevailing wind direction (Fig. 1).
About every two meters, the elevation of the fields above mean sea level was measured with a differential global positioning system (Altus APS 3M, Septentrio, Belgium). 20

Statistical analyses
For the statistical analyses, we used all available data on energy fluxes from the onset of measurements (late March/early April) until harvest. In the case of maize, however, full data for the calculation of energy balances was generally available from May. Autocorrelation of the data was tested using the Durbin-Watson test (Faraway, 2014). Analysis of variance 25 (ANOVA) was used to test for significant effects of region, site, year and crops on EBC, for which linear mixed models were defined (Piepho et al., 2004). The data were assumed to be normally distributed but heteroscedastic due to the different years.
We based these assumptions on graphical residual analyses. Generally, the factors of interest were defined as fixed and interaction terms were considered. Remaining factors not included in the ANOVA were defined as random. Multiple contrast tests (Bretz et al., 2011) were performed to identify significant differences between the different factor levels. Unless indicated 30 otherwise, the significance level was set to = 0.05. Calculations were done using the statistical software R (R Core Team, 2014) and packages multcomp for simultaneous tests of linear mixed models (Hothorn et al., 2017), nlme for fitting and comparing the models (Pinheiro et al., 2016), gplots for creating plots (Gregory et al., 2009), and gdata for importing input data from MS Excel formatted files (Gregory et al., 2017).

Kraichgau
At the KR sites, the annual mean air temperature ranged between 8.4 °C in 2010 and 10.9 °C in 2014. The overall average was 9.8 °C (Fig. 2a), which is 0.4 °C higher than the 30-y climatological mean (1981−2010) measured at the meteorological station Pforzheim-Inspringen. The lowest and highest monthly mean temperature was −3.2 °C in February 2012 and 21.1 °C in July 2015. The mean annual precipitation was 796 mm, which is 93 mm lower than measured in Pforzheim-Inspringen. In 2013, 10 the wettest year within the 8 year period, total precipitation amounted to 973 mm. The lowest annual precipitation (629 mm) was measured in 2015 ( Fig. 2b). Fig. 3 shows the height transects along the prevailing wind direction. At the KR sites, the mean slopes along the transects were 0.4 %, <0.01 % and 0.3 % at EC1, EC2 and EC3, respectively. The micro-relief of station EC1, located on a micro-bank, fluctuates more strongly than that of EC2. The immediate surroundings of EC2 are very homogeneous in elevation. Station 15 EC3 was positioned in a micro-depression. Overall, the three transects show that the KR fields can be regarded as flat, with EC2 being the flattest.

Swabian Jura
The mean temperature in SJ (7.4 °C) was 2.4 °C lower than in KR varying from 5.9 °C in 2010 to 8.5 °C in 2015 (Fig. 2c). As 20 in KR, the lowest and highest mean monthly temperatures were recorded in February 2012 (−6.6 °C) and July 2015 (18.6 °C).
The mean annual precipitation was 874 mm. As in KR, 2015 was the year with the lowest precipitation. Highest total rainfall was measured in 2017, not in 2013 as in KR. November 2011 was the month with the lowest monthly cumulative precipitation (5 mm), and July 2014 that with the highest (187 mm) (Fig. 2d).
In SJ, only EC4 is relatively flat (Fig. 3). Its topography is comparable with that of EC1 in KR. The elevation along the transect 25 at EC5 gently increases, with a mean slope of 0.6 % from SE to NW. Station EC5 itself is situated in a local micro-depression.
The topography of station EC6 differs considerably from that of the other fields. The station is positioned on the top of a ridge.
Whereas in NW direction the terrain drops with a mean rate of 3.7 m per 100 m, in SE direction the terrain is nearly flat (slope = 0.3 %).

Energy partitioning at the land surface
The energy partitioning at the canopy surface of different crop stands is shown, by way of example, for the vegetation period of 2016 (Fig. 4). In that year, five different crops (winter rapeseed (WR), spring barley (SB), winter wheat (WW), silage (SM) and grain maize (GM) were grown at the EC sites. From April to June, most of the net radiation was transformed into latent 5 heat at the crop stands, except for SM at EC5. The daytime Bowen ratio was lowest for WW and WR, with 0.14 and 0.13, respectively. Also GM, SB, and SM at EC6 led to daytime Bowen ratios distinctly below unity (about 0.21). Only silage maize at EC5 had a Bowen ratio of about unity, which indicates that the available energy was partitioned into latent and sensible heat in similar proportions. At the WW, SB, and WR sites and years, the ground and sensible heat fluxes were nearly the same and showed a similar diurnal course. At the maize stands, the ground heat flux tended to be higher than the sensible heat flux during 10 the morning hours, while in the afternoon the order switched and more sensible heat than ground heat was formed. At all sites, the measured energy residual was similar to the sensible heat fluxes, ranging from 23 W m −2 at EC3 to 44 W m −2 at EC1. The daily net radiation was 149, 133, 134, 130, 138, and 164 W m −2 at EC1 to EC6, respectively. The mean daily LE ranged from 54 W m −2 at EC5 to 94 W m −2 at EC3.
For July to September, the strongest shift in energy partitioning occurred at the WR site. Over noon, the Bowen ratio was in 15 the range of unity, and sometimes the half-hourly sensible heat flux was even higher than the latent heat flux. A similar shift was observed at the WW site, but weaker than at the WR site. At the GM, SM and SB sites the largest difference compared with the period April to June was the ratio between sensible and ground heat flux. From July to September the sensible heat was about twice the ground heat flux. The mean net radiation ranged from 125 W m −2 at EC5 to 176 W m −2 at EC2, and LE varied from 78 W m −2 at EC4 to 89 W m −2 at EC1. The residual energy for this period was 23, 28, 22, 24 and 16 W m −2 at the 20 sites EC1, EC2, EC3, EC4, and EC6, respectively. Note: EC5 data are missing due to a damage in the sonic anemometer and gas analyser.

Energy balance closure
The mean EBR over the 48 site-years was 0.75, corresponding to a mean energy residual of 41.6 W m −2 (Tab.  In general, the EC method performed best (EBC was highest) over the vegetation period from April to August. The highest EBC was usually found during July and August, and distinctly declined over autumn and winter. At station EC6, the SJ station equipped with a fuel cell system, for example, the EBC declined to 42 % in January 2016 and 23 % in December 2017. A low EBC was usually associated with a larger variation (see winter months; Fig. 5).

Effect of region, station, year and crop
The statistical analyses showed that the EBC over the main vegetation period from early April until harvest did not differ between the two regions (Fig. 6a). The EBC at stations EC2 and EC4 was significantly higher (p < 0.001; pprobability level) than at the other stations (Fig. 6b). The lowest spread in values was observed at station EC4. In 2013 and 2014, EBC was lower (p < 0.001) than in the other six years (Fig. 6c). The crops had no significant effect on mean EBC (Fig. 6d). EBC over winter 10 rapeseed showed the highest variation in comparison to the other four crops, varying between 57 and 88 %.

Effect of wind speed and direction
Typical for the mid latitudes, the KR sites' prevailing wind direction was from west to east. The fraction of WSW to WNW (240° -300°) winds was 43.2, 36.8 and 33.7 % at EC1, EC2 and EC3, respectively (Fig. 7). The highest wind speeds were 15 also measured within these wind direction sectors. Wind blowing from north-and southward directions was rarely measured (< 10 %). While at EC1 the wind speed averaged 2.9 m s −1 , at EC2 and EC3 the values were 2.4 and 1.9 m s −1 , respectively.
Moreover, high wind speeds (> 6 m s −1 ) clearly decreased from the most westerly station (EC1) to the most easterly station (EC3). At EC1, EC2 and EC3, the share of these high wind speeds was 4.6, 3.3 and 0.2 %, respectively.
The distribution of EBR as a function of wind direction is shown in Fig. 8. The EBR was averaged for 30° wind sectors over all available daytime data (global radiation >10 W m −2 ). At the KR sites, the highest EBR was achieved when the wind blew 25 from westerly directions, which is the prevailing wind direction. Wind from northern and southern directions was related to lower EBR. Particularly wind from the south was associated with EBR below 0.6 at all three stations. This phenomenon was most pronounced at station EC1. Also at the SJ sites, highest EBR was usually coincided with wind from the prevailing direction. One exception is the high EBR at EC4 for the SSE sector. At EC4, for six out of the twelve wind sectors, EBR was above 0.8. In contrast, at EC5 and EC6 the EBR exceeded 0.8 in only three wind sectors. At all stations, EBR was lowest (< 0.6) for winds from north-east (Fig. 9).

Effect of atmospheric conditions
This section evaluates the EBR as a function of buoyancy, shear and atmospheric stability. For this purpose, we plotted EBR 5 against the kinematic virtual temperature flux (w'Tv'; proxy for buoyancy), friction velocity (u * , proxy for shear) and the stability parameter ζ (Fig. 10). Again, only half-hourly daytime fluxes (Rs > 10 W m −2 ) were evaluated. The plot EBR versus w'Tv' shows a vast scatter at weak buoyancy. Here, the EBR ranges from plus four to minus four. The scatter decreases substantially as the modulus of w'Tv' increases. Note that w'Tv' < −0.15 K m s −1 were measured only at stations EC2 and EC4.
Plotting EBR against friction velocity also reveals a large scatter, which narrows as friction velocity (shear) increases. The 10 scatter, however, does not narrow as much as for increasing buoyancy. During neutral or stable atmospheric conditions, the EBR showed a large spread (Fig. 10). In contrast, this range distinctly declined when the stability parameter reached strongly negative values, indicative for highly unstable conditions. EBR above unity or below zero was rarely observed under these conditions.
From the total dataset, only 7 % of daytime measurements were made under stable conditions, 34 % under unstable conditions 15 and 59 % under neutral conditions (Table 4). During unstable conditions, EBC and EBR at all sites were slightly lower compared to neutral conditions. During neutral conditions, however, the standard deviation (SD) of EBR was about twice as high as under unstable conditions. During stable conditions, EBC and EBR were systematically lower than at unstable and neutral conditions. At EC4, for example, the EBR under neutral and unstable conditions was 0.78 and 0.85, respectively. Under stable conditions, the value declined to 0.41. Moreover, the huge spread in EBR under stable conditions is underlined by its 20 high SD, which is about three times the mean value.

EBC and energy balance components
From July to September, daily mean Rn varied between 125 and 176 W m −2 . Similar ranges of Rn were observed over maize in Livraga, Italy (Masseroni et al., 2014). The latent and sensible heat fluxes varied strongly over the observational period. In early-covering crops (winter rapeseed, winter wheat, winter barley), LE was about two to three times higher than H in the 5 period AMJ (April-May-June), while in the period JAS (July-August-September) LE and H were in a similar range (Fig. 4, EC3-WW, EC4-SB). The period JAS is when ripening and harvest of cereals and winter rapeseed occurs, as well as postharvest management such as tillage and seeding of cover crops or winter rapeseed. During AMJ, the patterns of LE and H at EC5 and EC6 differed, even though the maize was sown on similar days of the year (May 07 at EC5 and May 03 at EC6). This can be explained by the substantially higher leaf area index at EC6 (0.74±0.15) compared to EC5 (0.35 ±0.06), measured on 10 June 22.
The mean EBR of the 48 site-years was 0.75 (Table 3). In comparison, Wilson et al. (2002) reported an EBR of 0.84, on average, for the 50 analyzed FLUXNET site-years, ranging from 0.34 to 1.69. In three agricultural and one industrial site in South Korea, the mean value varied between 0.46 and 0.83 (Kim et al., 2014). Majozi et al. (2017) found a mean EBR of 0.93 at a semi-arid savannah site in South Africa, over a period of 15 years. 15 The slopes of the OLR and EBR differed by a maximum of 5 %, which is consistent with previously published data. Such a small difference points at a high reliability of the presented EC measurements (Wilson et al., 2002). The highest annual EBC occurred at EC4 (87 % in 2010), the second highest at EC2 (83 % in 2016), the lowest at EC5 (62 % in 2016). The lowest EBC was observed mainly in the cold, non-growing season, which may be attributed to insufficient thermally and mechanically induced turbulence (Franssen et al., 2010) as well as to freezing (Varmaghani et al., 2016). 20 The incomplete EBC in our dataset has several potential explanations. One is related to the neglected minor storage terms (Eshonkulov et al., 2019;Masseroni et al., 2014;Meyers and Hollinger, 2004). Importantly, considering minor storage terms is not straightforward because they are not measured when conventional EC equipment is used. Only the energy fixed and released by photosynthesis and respiration can be directly derived from EC data because the net CO2 flux is generally measured. Considering minor storage terms in calculating the EBC at a maize field improved the mean value from 87 to 91 % 25 (Xu et al., 2017) and from 81 to 86 % (Masseroni et al., 2014). Eshonkulov et al. (2019) demonstrated that the contribution of minor storage and flux terms over winter wheat in southwest Germany was largest during the main vegetation period in May.
During this month the minor terms helped to close the energy balance by an additional 7−8 %.

The effect of meteorological conditions and surface-layer turbulent parameters
In both KR and SJ, EBR was highest for winds blowing from the prevailing wind direction. These winds were associated with high wind speeds favouring well-developed turbulent conditions. This is consistent with other studies. Xin et al. (2018), for example, also found that winds with high speeds blowing from the prevailing direction yielded consistently higher EBC compared to other directions. Kim et al. (2014), for example, grouped EBR into two different categories, with lower (< 0.75) 5 and higher EBR (> 0.75) and for their four research sites observed that EBR was higher at high wind speed. Fig. 10 shows that the spread of EBR distinctly narrowed at high friction velocities (u* ≥ 0.5). Prior studies have noted the importance of u* on the EBC. Anderson and Wang (2014) found that, under these conditions, EBC was closed on days with continuous turbulence. Results of hourly daytime EBR and u* showed a strong relationship at our sites (Fig. 10). This is consistent with other studies carried out in selected croplands such as irrigated sugarcane (Anderson and Wang, 2014), maize 10 plantations (Masseroni et al., 2014) and rice fields (Kim et al., 2014). Sánchez et al. (2010) also reported that EBR was > 0.90 when high friction velocities prevailed (> 0.8 m s −1 ) at a boreal forest site in Finland. Mauder et al. (2013) investigated EBC at the TERENO site Lackenberg (Germany) and found that it was almost closed. They explain this result by the very good turbulent mixing and the high homogeneity at this site. This confirms that, at high u*, the production of high-frequency fluxes is elevated (Fratini and Mauder, 2014). 15 At our study sites, neutral conditions dominated (~ 60 %), followed by unstable conditions (~ 34 %) and by stable conditions (6 %) (Table 4). Importantly, average EBR changed from 0.67 (± 0.32) to 0.72 (± 0.69) and 0.41 (± 1.33) during unstable, neutral, and stable conditions, respectively (SD in brackets). Under stable conditions, the EBR was lowest and had the largest variation. Averaged over all EC stations, the slope of OLR under neutral conditions was slightly higher than under unstable conditions. This is also evident in the mean and variance of the calculated energy residuals. The average residuals under stable, 20 neutral and unstable conditions were 9.7 (± 31.5), 47.5 (± 58.2) and 81.5 (± 62.0) W m −2 , respectively. The coefficient of variation was highest under stable conditions and decreased over neutral to unstable conditions. This result differs from previous studies. Mauder et al. (2010) reported a residual energy close to zero for a cropland in Ontario, Canada, under stable conditions, peaking at 150 W m −2 under neutral, and decreasing to 100 W m −2 under unstable conditions. The scatter of EBR versus buoyancy flux at EC2 and EC4, the two stations with the highest EBC, differed from those of the 25 other stations (Fig. 10). At these two sites, strong negative buoyancy fluxes below −0.15 K m s −1 were recorded. This means that the atmosphere was not heated by the land surface, but that the land surface was significantly heated by the atmosphere.
Such a situation points to a stable boundary layer (SBL). Lan et al. (2018) report that they measured the highest buoyancy fluxes under a weak SBL with strong surface shear. They argue that the strong mechanical shear produced at the ground favours the development of turbulent eddies with larger scales that enhance vertical mixing of momentum and heat transporting the 30 aloft warm air downward and the surface cold air upward. Moreover, the mechanical mixing weakens the magnitude of mean temperature gradient and allows turbulent eddies with larger vertical scales to develop. Conversely, under a SBL weak winds occur near the surface and turbulent eddies are depressed and detached from the boundary leading to suppressed vertical mixing.
Several studies recommended considering secondary circulations to achieve a better EBC (Foken et al., 2010;Kidston et al., 2010;Mauder et al., 2010). Those studies postulate that heterogeneity-induced and buoyancy-driven quasi-stationary circulations are probably the dominant processes behind underestimated energy fluxes. The studies that suggested the use of 5 an averaging period higher than 30 minutes usually refer to unstable conditions. These studies suggested that averaging periods of 2-4 h are often needed to statistically resolve the largest convective turbulent eddies or also non-stationary mesoscale motions that sometimes can modulate turbulent fluxes (Mahrt, 1998). Larger averaging improved short term EBC during the diurnal hours in the Salentum peninsula of Apulia, Italy (Cava et al., 2008). In considering secondary circulations, different time averaging intervals instead of the standard 30-min period can be used. Although a 60-min interval might be suitable for 10 capturing the major turbulent fluxes (Kilinc et al., 2012), in most cases the standard 30-minute period is still sufficient (Kidston et al., 2010). The classical averaging period of 30 minutes can be a proper choice for unstable or neutral conditions. Shorter averaging period is suitable for capturing energy fluxes in very stable conditions (Sun et al., 2012;Vickers and Mahrt, 2006).
Finding an optimum averaging period is a very complex and nearly impossible task. This is because atmospheric turbulence changes irregularly and there is no clear-cut "switch" in time. Therefore, the averaging time could be modified during raw data 15 processing. In practice, however, this is unlikely because it drastically increases the complexity of data processing (Lenschow et al., 1994). Moreover, the sources of secondary circulations are unclear, and they are most probably not well linked with the locally measured available energy. Accordingly, excluding secondary circulations in EC measurements can be locally meaningful. Recently, a new method, known as ogive optimization, was proposed by Sievers et al. (2015). The method enables separating low-frequency influences from vertical turbulent fluxes to isolate the local exchange processes of interest. 20 Although EC measurements contain uncaptured energy components, the flux data are used, among others, to evaluate models and interpret simulation results. In such studies, EC flux data are usually post-closed, i.e., the measured turbulent fluxes are adjusted so as to close the energy balance (Ingwersen et al., 2015). The standard approach is the Bowen-ratio post-closure method (Twine et al., 2000). It assumes that the missing energy has the same Bowen ratio as the measured turbulent fluxes.
This approach, however, may introduce a systematic bias to simulated surface energy fluxes (Chen and Li, 2012). Analyses of 25 the energy imbalance by Ingwersen et al. (2011) showed that soil water contents simulated by a land surface model agreed better with measurements when the residual was fully assigned to H. As discussed by Charuchittipan et al. (2014), secondary near-surface circulations attributed to low frequencies mainly transport sensible heat. Therefore, they proposed a new alternative energy balance correction method they termed the Buoyancy flux ratio. At very large Bowen ratios (> 10), the Bowen ratio post-closure and buoyancy flux correction methods yield similar results. At Bowen ratios ranging from 0.1 to 0.2, 30 which is typical for croplands during the main growing period, the Buoyancy flux ratio method assigns most of the energy residual (> 50 %) to the sensible heat flux. The Bowen ratio method, in contrast, distributes most of it (> 90 %) to latent heat.
As long as the composition of the residual remains unknown, it is important to communicate the possible error in EC flux data, for example with the post-closure method uncertainty band (PUB) (Ingwersen et al., 2015). Working with only one postclosure method may result in serious misinterpretations in model-data comparisons.

The effect of the instrumental setup
At the SJ sites, we found particularly low EBR in the wind sector 0−90°. The CSAT3 sensor was oriented mostly to 225°, so 5 that the sector 30−90° was located behind the anemometer head. To substantiate the idea that the anemometer negatively influences EC measurement quality, and taking the data from EC4 as an example, we recalculated EBR across all years,

Relationships between EBC and footprint
Accurate measurements of energy balance components are important to achieve a good EBC. In this context, one key requirement is that the EC station be located in a place that represents the fluxes from the area of interest (Burba and Anderson, 2010). According to those authors, the terrain must be horizontal and uniform. Three parameters are needed in footprint analysis: measurement height, surface roughness and atmospheric stability. When turbulent fluxes originate from a horizontal 20 and homogeneous surface, the footprint depends solely on the distance between the location of the measurement point and the emission element. We found a distinct tendency that the smaller the footprint, the higher the EBC. We express two explanations. First, the smaller the footprint, the higher the chance that the assumption of a homogeneous source area is fulfilled. Second, the smaller the footprint, the better the scale match between the measurement of available energy and turbulent fluxes. Alfieri and Blanken (2012) found that variations of surface energy fluxes over tens of meters ranged from 30 25 to 40 W m −2 using single-point (immobile) and mobile EC towers at a uniform site (Colorado, USA). They concluded that a single-point EC tower cannot capture all the relevant energy fluxes because they vary spatially. Our results confirm that if the footprint is small, the EBC from EC measurements is better, which can be interpreted as a reduction in the variation of surface energy fluxes.
Many studies claimed that surface heterogeneity is a potential reason for the energy imbalance (Stoy et al., 2013;Xu et al., 30 2017). The latter authors reported that EBC decreased with increasing surface heterogeneity. The degree of heterogeneity was derived from high-resolution remote sensing images and land surface temperatures. To handle this effect, some authors recommend using direction-specific coefficients that indicate the degree of heterogeneity. For example, Panin et al. (1998) introduced a heterogeneity factor that comprises surface parameters such as roughness, radiation and the thermal humidity of the internal boundary layer. That factor can be used for data interpretation. Nonetheless, deploying heterogeneity factors still does not explain how the residual energy is composed. The lack of EBC at the KR sites might partly reflect katabatic advection 5 (Heinesch et al., 2008;Kutsch et al., 2008), which results from stable atmospheric conditions and occurs especially in hillslope areas (Loescher et al., 2006;Mauder et al., 2010). Moreover, complex topography can induce advective fluxes (Feigenwinter et al., 2008;Rebmann et al., 2010). The former landfill site located about 500 m south of the fields in KR (Fig.1) might have been responsible for advective fluxes, since it's elevation is approximately 41 m higher than the study sites. Moreover, the topography could also affect EBC. The elevation transects show that the immediate terrain surrounding the stations EC3, EC5 10 and EC6 is not totally flat (Fig. 4). This is a well-known problem for micro-meteorological field measurements (Wilczak et al., 2001). At EC1, EC2, and EC4, however, the terrain can be considered flat.

Conclusions
We evaluated the EBC of long-term EC measurements at six different cropland sites in two contrasting environmental regions 15 in southwest Germany. EBC depended on how well thermally and mechanically induced turbulence was developed. On average, 25 % of the available energy was not detected by our EC stations, with the lowest annual imbalance (energy residual) of 17 % in KR and 13 % in SJ. This range of EBC is common in cropland, and such recovery rates must be accepted in heterogeneous landscapes. We interpret the range of the highest mean annual EBC (83 % resp. 87 %) as the upper detection limit of the EC method at our sites and settings. During winter months and under stable atmospheric condition, EBC was 20 problematic. EBC was negatively affected by: i) stable atmospheric conditions, ii) non-horizontal or heterogeneous source area, iii) larger obstacles in the landscape, i.e. the former landfill site that may have induced adjective flux components, and iv) flow distortions of winds first travelled past the back head of the anemometer, which reduces wind speed and corrupts the spectral characteristics of turbulence at specific wind directions. EBC was positively affected as the footprint area decreased, probably because this tends to decrease the heterogeneity of the source area and improves the match of available energy 25 measured locally with the mean available energy in the footprint. We thank Dr. Paul Stoy for handling the manuscript, one anonymous reviewer and Marcelo Zeri for helpful and constructive 5 comments.
Hermann (EC5) and Mr. Reichart (EC6) in SJ for the permission to conduct measurements in their fields. The authors would also like to thank the technical staff Benedikt Prechter, Felix Baur, Christian Schade, and Thomas Schreiber.