1 Introduction

The energy balance for land-surface cover follows the first law of thermodynamics, where the net radiation (\({R}_{n}\)), the soil heat flux (G), and the change of heat storage (S) within the system is partitioned into turbulent fluxes of latent (\(LE\)) and sensible heat (\(H\)) with the atmosphere, i.e.

$${R}_{n}-G- S= H+ LE.$$
(1)

Given that the incoming radiation remains the same over a larger area, the available energy (left-hand side of Eq. 1) of different land-surface covers with scanty soil and vegetation is mainly controlled by their albedos and surface temperatures. This implies that heterogeneous land surfaces consisting of a mix of different land-surface cover potentially experience substantially different microclimates (Rotach and Calanca 2002). To measure the exchange of heat, momentum, and trace gases between the atmosphere and Earth’s surface at a micrometeorological scale, the eddy-covariance method (Baldocchi et al. 1988; Aubinet et al. 1999) has shown to be the most accurate and reliable technique (Finkelstein and Sims 2001). The technique provides a considerable contribution to the study of environmental, biological, and climatic controls between the atmosphere and different terrestrial ecosystems (Baldocchi et al. 2001). However, an independent assessment of the turbulent eddy-covariance flux estimates is important for the evaluation of their reliability and validity. The energy balance closure (EBC) is widely used and considered as a standard method for this purpose (e.g., Anderson et al. 1984; Verma et al. 1986; Hollinger et al. 1999; Foken et al. 2010; McGloin et al. 2018). Most eddy-covariance sites experience a non-closure of the energy balance with larger available energy than turbulent energy. The lack of EBC has been heavily studied and reported in the last 25 years, as shown in the review papers of Foken (2008) and Mauder et al. (2020). Even for flat and homogeneous ideal surfaces, the energy balance is rarely closed (Stannard et al. 1994; Xin et al. 2018). Stoy et al. (2013) found, among other things, a general tendency of increasing energy imbalance with increasing landscape heterogeneity in vegetation cover and topography.

At eddy-covariance sites with heterogeneous surface cover, one of the additional challenges that arises is the mismatch of footprints of the different sensors used to measure the different energy fluxes. Typically, the horizontal footprint scale ranges from 0.1 m for the ground heat flux, to 10 m for net radiation, and up to several hundreds of metres for the turbulent fluxes (Mauder et al. 2020). Uncritical comparison of the available and turbulent energy where spatial variation of surface cover exists within the different sensors’ footprints could therefore potentially aggravate the EBC, or artificially improve it. At eddy-covariance sites characterized by considerable heterogeneity, net radiation measurements should therefore be corrected for the surface-cover variability to improve the compliance with the surface of the turbulent energy exchange (Schmid 1997). As such, fixed flux towers and regular remote-sensing satellites fail to meet the needs of fine-spatial radiation measurements.

The objective of this study is to develop and test a methodology (Fig. 1) for providing surface-cover-corrected net radiation at a heterogeneous eddy-covariance site. Since a tower-based radiometer only provides a one-point measurement of outgoing radiation, we hypothesize that a surface-cover-corrected net radiation serves as a more reliable estimate of the available energy and therefore plays a considerable role in explaining the lack of energy balance when surface heterogeneity is present. The study was conducted at an alpine eddy-covariance site in mountainous Norway with both complex and heterogeneous terrain and surface cover. The site is not ideal for eddy-covariance measurements, however, it is important for improving the knowledge of energy exchange processes in complex mountainous ecosystems, which have been rarely investigated and studied (McGloin et al. 2018).

Fig. 1
figure 1

Key steps of methodology for providing surface-cover-corrected net radiation matching the footprint of turbulent fluxes at a heterogeneous eddy-covariance site

Full size image

2 Materials and Methods

2.1 Study Area

The study site (60°36′N 7°30′E) is situated in an alpine valley at an elevation of around 1200 m above sea level (Fig. 2a.) near the Finse Alpine Research Center in southern central Norway. The valley extends in the north-west–south-east direction, where Hardangerjøkulen Glacier (~ 1850 m above sea level) is located south-west of the valley. There is a confluence of cool glacial meltwater and warmer non-glacial streams from lake Finsevatnet in the river Utsekveikja, which runs through the study site. The climate is arctic and features maritime influences with relatively mild winters, while the summers are cool. This is a result of the transition zone between the easterly continental and westerly oceanic climate, where warm, moist westerly winds from the Atlantic Ocean dominate (Leinaas and Schumacher 2005). The annual mean (1991–2020) temperature is −1.1 °C and annual average precipitation is 967 mm. The typical low summer temperature is somewhat explained by the influence of the glacier (Leinaas and Schumacher 2005), where cold air tends to sink down to the valley. The surface cover is dominated by lichen heaths on wind-exposed ridges, while narrow zones of dwarf shrubs and mountain heather are found on the leesides. Further downslope, the snow cover increases, and snow-beds with bare rocks and mosses dominate until the topography flattens out towards the wetlands and the river. In the flat areas, water accumulates to generate mires and small ponds. Along the river and streams there are narrow bands of dwarf willows.

Fig. 2
figure 2

a Location of the study site in south-west Norway (black square); b The mobile radiation tower and typical surface cover at the study site

Full size image

2.2 Flux Tower Instrumentation

The exchange of sensible and latent heat between the surface cover and the atmosphere was measured using the eddy-covariance method (Baldocchi et al. 1988; Aubinet et al. 1999). The three velocity components and the speed of sound were measured by a three-dimensional sonic anemometer (CSAT3, Campbell Scientific, U.S.A.) installed at the northern end of a horizontal boom at 4.4-m height above ground. Next to this, an enclosed path infrared gas analyzer (LI7200RS, LiCor, U.S.A.) with a heated intake tube of 0.71-m length was installed to measure the concentration of water vapour (and carbon dioxide). The eddy-covariance system sampled with a frequency of 10 Hz and was connected to a flash drive, on which the data were stored. The EddyPro v6.2.0 (LiCor Inc, 2016) software was used to process the raw data. First, a double-rotation method was applied, so that the assumption of negligible average vertical velocity component was fulfilled. After extracting the turbulent fluctuations from the raw signal, 30-min fluxes of sensible and latent heat were estimated from the covariance between the turbulent fluctuations of the vertical velocity component and the temperature and water vapour concentration, respectively. Spectral corrections were done to compensate for both high- and low-frequency loss (Moncrieff et al. 1997, 2004). We assessed the impact of high-frequency corrections by comparing the fluxes to calculations using the alternative in situ correction of Fratini et al. (2012), which is based on the combination of a direct approach and the analytical formulation of Ibrom et al. (2007). This correction considers the effect of relative humidity on the signal dampening at high frequencies. During unstable conditions, the in situ method gave 5–10% larger values for LE (and no effect on H), but it also introduced unrealistically large variability (often exceeding a value of 5 compared to uncorrected fluxes), especially during stable conditions, typically during the nights. To get the most homogeneous and comparable time series, it was therefore decided to use the simpler, but more robust, analytical method for high-frequency spectral corrections. To assure the quality of raw data, statistical tests suggested by Vickers and Mahrt (1997) were applied. Furthermore, we assured the fulfilment of the assumption of steady state and fully developed turbulence following Foken and Wichura (1996). All the 30-min estimates without the highest measure of quality (i.e., turbulent flux estimates with quality flag 1 or 2, as given by EddyPro) were discarded. Also, the 30-min estimates of the turbulent fluxes where the mean wind speed was less than 2 m s−1 (corresponding to a friction velocity of about 0.2 m s−1) were discarded to avoid totally windless periods. Most of the 30-min observations had ratios close to 1:10 between the friction velocity and mean wind speed, which therefore fulfilled this filter criterion. Filtering the turbulent data according to these procedures, resulted in a data coverage of 59% for the time during the study period.

One-point tower-based net radiation (\({R}_{n}\)) was measured by a four-component radiometer (CNR4, Kipp & Zonen, The Netherlands), which was installed at the southern end of the horizontal boom at 4.4-m height above ground. Since the upward-looking pyrgeometer was defective during the field survey, the incoming longwave radiation used to calculate \({R}_{n}\) was measured by another CNR4-radiometer located 250 m south-east of the eddy-covariance tower. Soil heat fluxes (G) were measured in proximity to the tower by two soil-heat-flux plates (HFP01, Hukseflux, The Netherlands) installed approximately 0.1 m below surface and 0.5 m apart from each other. The measurements were sampled at 0.2 Hz and stored on a datalogger (CR6, Campbell Scientific, U.S.A.), before being averaged to 30-min fluxes. See Sect. 3.2 for further discussion of challenges and limitations related to the flux measurements.

2.3 Methodology for Providing Surface-Cover-Corrected Net Radiation

2.3.1 High-Resolution Vegetation Mapping

The vegetation of the study site was mapped in 2017 by using the Nature in Norway (NiN) system. The NiN system is a full coverage system that describes nature variation in Norway. The mapping at Finse was done at a spatial scale of 1:5000, where the (recommended) smallest distance of the cross section of a mapping unit was 7.5 m. However, for nature types occurring at spatial resolutions smaller than the mapping scale, the mapping unit was described as a mix (either as a regular or non-regular mix) of the relevant nature types (Bryn et al. 2018). The high-resolution vegetation mapping made it possible to characterize the footprint of the flux tower in terms of vegetation. This information is important for understanding which types of vegetation and surface covers contribute the most to the measurements of the turbulent fluxes. In total, the study site consisted of 32 different nature types (Bryn and Horvath 2020). Nature types with similar vegetation structure, vegetation function, surface brightness, and approximately same access and ability to store and capture water were categorized into the same surface type. This assured that the surface types consisted of nature types having approximately the same albedo and surface characteristics controlling the longwave radiation balance. The categorization resulted in 11 surface types. However, one of these consisted of nature types strongly affected by human disturbances (buildings and gravel roads covering an area that contributed on average about 1.4% of the eddy-covariance flux signal) and was therefore not of interest. The area of each surface type was calculated to detect the percentage contribution of the study site. This revealed that three of the surface types only accounted for 1.07%, 0.21%, and 1.74% of the total mapped area of the study site. Due to their insignificant contributions, these surface types were excluded for the purpose of providing in situ surface-cover-specific net radiation measurements. The remaining seven surface types for which net radiation was measured are shown in Fig. 3. Their nature types are described in Table 1.

Fig. 3
figure 3

The eddy-covariance site with surface cover classified into surface types as given by the colour coding in the legend. Each surface type consists of nature types as given in Table 1. The black crosses show the location of the measuring points where net radiation was measured by the mobile radiation tower. The black circle displays the position of the eddy-covariance tower. The two shaded areas illustrate the difference of surface cover contributing to two of the 80% 30-min footprints of the turbulent fluxes for wind directions from north-west and east, respectively

Full size image
Table 1 Classification of the seven surface types based on the NiN nature types with their descriptions
Full size table

2.3.2 Surface-Cover-Specific Net Radiation

The in situ measurements of net radiation of the different surface types were performed during a field survey from 18 to 25 August 2018. A mobile radiation tower (Fig. 2b) was specially built to allow necessary equipment to be easily moved around between the different surface types. The mobile radiation tower consisted of a four-component net radiometer (CNR1, Kipp & Zonen, The Netherlands), which was connected to a battery-powered datalogger (CR6, Campbell Scientific, U.S.A.) with a logging frequency of 0.2 Hz. The pyrgeometers (spectral range of 4500–42,000 nm) had 150° field of view, while the pyranometers (spectral range 305–2800 nm) had 180° field of view. The manufacturer of the net radiometer suggested that the radiometer should be mounted at a height of at least 1.5 m above ground to minimize shading effects of the instruments on the surface. However, this height was found to correspond to a spatial scale larger than most of the extent of the patches of the different surface types. A mounting height of 1.1 m was therefore chosen as a compromise to ensure homogeneous footprints within each surface type.

Two to three patches with the largest extent within each surface type were chosen as suitable measuring points. The tower was placed at a measuring point for 4 h at a time before it was moved to the next measuring point. This procedure was performed during the daytime from 0800 to 2000 LT (local time = UTC + 2 h). During night-time, from 2000 to 0800 LT, the tower was placed at the same measuring point. Following this schedule, it was possible to measure net radiation within each of the seven different surface types evenly distributed throughout the whole field survey. This procedure also ensured that all surface types were represented by at least one measurement from all times of the day. The radiometer was placed parallel to the surface by the same reference to gravity for each measuring point. A compass was used to ensure that the radiometer always pointed to the south.

2.3.3 Predicting Time Series of Net Radiation

Continuous time series of net radiation for each surface type were provided by fitting linear regression models. Linear regression models have been widely used to predict net radiation based on radiative and other meteorological parameters (e.g., Iziomon et al. 2000; Kjaersgaard et al. 2007). However, the predictions of various regression models are dependent on the local fitting of the regression coefficients, which implies that regression models are dependent on the surface type and its radiative conditions, thus only valid for a restricted area. The model form used for predicting cover-specific net radiation (\({R}_{n}^{*})\) for each of the seven individual surface types was

$${R}_{n}^{*}=\alpha +{\beta }_{1} {R}_{n}+{\beta }_{2} {SW}^{\downarrow *}+{\beta }_{x}{ LW}^{\downarrow *} + \varepsilon $$
(2)

The left-hand side of the equation refers to net radiation obtained by the mobile radiation tower. Here \({R}_{n}\) represents the one-point net radiation as measured by the eddy-covariance tower, with model parameter \({\beta }_{1}\), \({SW}^{\downarrow *}\) and \({LW}^{\downarrow *}\) represent the surface-cover-specific incoming shortwave and incoming longwave radiation, respectively, as measured by the mobile radiation tower, α, \({\beta }_{2}\), and \({\beta }_{3}\) are thus the local fitted cover-specific regression coefficients for each individual surface type, and act as surface-dependent corrections of the general tower-based net radiation, and \(\varepsilon \) represents the error term, which was assumed to be uncorrelated and normally distributed with expectation zero and constant variance. The outgoing longwave radiation for each specific surface type is largely constrained by the surface temperature and necessarily not the incoming longwave radiation. An accurate prediction of the longwave net radiation might therefore be deficient if only relying on the incoming longwave radiation as an explanatory variable. By including the tower-based \({R}_{n}\) as an explanatory variable, we account for this effect by contributing with information about the relationship between the net longwave radiation and incoming longwave radiation.

All regression coefficients in Eq. 2 are assumed to be time independent, which is a simplification. The parameter \({\beta }_{2}\) is directly related to the surface albedo. Prominent temporal variations with larger albedo during large solar zenith angles and small albedo around solar noon, at least during clear-sky conditions with dominant direct radiation, should be expected. This diurnal time-dependency can be accounted for by including the interaction effect between the solar zenith angle and the incoming shortwave radiation as an explanatory variable. However, due to the prevailing cloudy weather with mostly diffuse radiation during our field survey (see Sect. 3.1), this effect was minimal and therefore omitted. This was substantiated by a pre-analysis, which revealed a Pearson correlation of 0.99 between \({SW}^{\downarrow *}\) and its product with the solar zenith angle.

Continuous time series of net radiation were predicted for each surface type by using the individual fitted regression coefficients and the continuous measurements of net radiation and incoming shortwave and longwave radiation as measured by the eddy-covariance flux tower. The accuracy of all the linear regression models was assessed by the adjusted coefficient of determination (\({R}_{adj}^{2}\)), the root-mean square (RMS) and the root-mean-square error (RMSE). Since all measurements were used as training data for fitting the linear regression models, leave-one-out cross validation (LOOCV) was conducted to assess the linear regression models’ suitability as predictive models. Their predictive power was reported in terms of root-mean-square-error-prediction (RMSEP) and the coefficient of determination for prediction (\({R}_{pred}^{2})\).

To test whether there was a significant effect of the surface characteristics on the net radiation between the different surface types, an analysis of covariance (ANCOVA) with \({R}_{n}\), \({SW}^{\downarrow *}\) and \({LW}^{\downarrow *}\) as predictor variables was performed. In addition, the interaction effects between measurement day and \({R}_{n}\), \({SW}^{\downarrow *}\) and \({LW}^{\downarrow *}\) were included in the ANCOVA model. This was done to be able to control for the effect of the different atmospheric illumination conditions and the between-day variation in the incoming radiation fluxes present during the field survey. Twenty-one pairwise comparisons were performed by using the package emmeans (Lenth et al. 2018) in the R software (R Core Team 2013) with test statistic set to the 5% level of significance. To compensate for the increased probability of erroneous rejection of the hypothesis of no difference between the surface types when performing several tests on the same data, the Bonferroni correction (Miller 1981) was applied for the multiple comparisons.

2.3.4 Footprint Analysis

The footprint of the turbulent fluxes may be visualized as the horizontal source area of the upwind direction of the surface containing effective sources and sinks contributing to the measurement point (Kljun et al. 2002). For a given measurement height of the eddy-covariance system, the footprint is dependent on the wind speed, wind direction, turbulence intensity, and atmospheric stability. The footprint of the eddy-covariance system was estimated following Kljun et al. (2015). The measurement height above ground, displacement height, boundary-layer height, and percentage source area of the footprint were set to 4.4 m, 0 m, 2000 m, and 80%, respectively. Additional input parameters were derived from the eddy-covariance data from the flux tower. Footprints were estimated for each averaged 30 min of input data. Within each 30-min estimated footprint, the percentage contributions, \({p}_{i}\), from each different surface type were calculated by dividing their respective areas with the total area of the footprint. By assuming a uniform weighting within the turbulent footprint, the weighted, surface-cover-corrected net radiation (\({\tilde{R }}_{n}\)) was found as

$${\tilde{R }}_{n} = {\sum }_{i=1}^{7}({p}_{i} {R}_{n, i})$$
(3)

where \({R}_{n, i}\) (predicted net radiation) now is \({R}_{n}^{*}\) for each of the \(i=1,\dots ,7\) different surface-cover types. Thus, \({\tilde{R }}_{n}\) represents the corrected value, which accounts for the surface heterogeneity within the footprint of the turbulent fluxes.

2.4 Assessment of Methodology

To test the effect of the methodology, we computed the EBC with both \({R}_{n}\) and \({\tilde{R }}_{n}\), one at a time, and compared the results. As such, we did not intend to close the energy balance, but assessed the effect of providing a more reliable estimate of net radiation for the heterogeneous surface cover. No measurements of the change of heat storage (S) were provided for this study, meaning that the EBC was assessed without it. The EBC was derived with three different methods. Firstly, the ordinary least squares (OLS) regression between the turbulent fluxes (H + LE) and the available energy (\({R}_{n}\) – G), was calculated. The OLS method assumes that there are no random errors in the independent variables of \({R}_{n}\) and G, which is a simplification (Meek et al. 1998). When neglecting this, the regression produces a downward bias in the slope coefficient and a flatter curve is obtained. To eliminate the effect of this error, the reduced major axis (RMA) regression between the turbulent fluxes and the available energy was derived as a second method. A lack of energy balance is often more pronounced at 30-min time scales compared to diurnal or seasonal scales. This is typically caused by systematic effects of overestimated positive daytime fluxes and underestimated negative night-time fluxes (Blanken et al. 1998). To account for missing estimates of energy stored during morning and early midday and energy released in the afternoon and evening, the energy balance ratio (EBR) was applied as a third assessment. The EBR is the ratio of the sum of turbulent fluxes to the sum of available energy (see, e.g., Wilson et al. 2002), and requires continuous measurement series within the averaging period. Therefore, to provide continuous data material of turbulent fluxes during the study period, discarded 30-min estimates of the turbulent fluxes were gap-filled separately for H and LE by using a random forest regression with 11 predictors (the predictor variables were incoming shortwave and longwave radiation, surface temperature, growing degree days, soil and air temperature, albedo, soil moisture, vapour pressure deficit, snow depth, and year). The model used 2000 decision trees and was trained on all available flux measurements (i.e., also including times before and after our field survey). Westerly footprints of the flux tower were used for the gap-filling because this was the dominating wind direction during the field survey. The models had coefficients of determination of 0.91 for the sensible heat flux and 0.88 for the latent heat flux. The EBR was calculated for the period from 19 to 25 August, starting and ending at 0830 LT.

3 Results and Discussion

3.1 Surface-Cover-Specific Net Radiation

Harsh weather conditions with cloudy weather and rapidly changing cloud cover prevailed during the field survey. In addition, two smaller storms occurred during the period. This resulted in sudden shifts in net radiation during the daytime, revealed as the spikes in the diurnal distribution of net radiation (Figs. 4 and 6). There were only a few periods of several hours with continuous clear-sky and sunny weather during the daytime. Therefore, larger differences in surface-cover-specific net radiation than reported here, especially for persistent, sunny daytime conditions, are expected to exist. Accumulated precipitation during the period was 46.3 mm. The predominant wind direction was west-north-west (Fig. 5).

Fig. 4
figure 4

Surface-cover-specific net radiation of the seven surface types, and the difference between the one-point tower-based net radiation (\({R}_{n}\)) and surface-cover-corrected net radiation (\({\tilde{R }}_{n}\))

Full size image
Fig. 5
figure 5

Wind rose for the period from 18 to 25 August 2018 based on all observations. The colours indicate the wind speed, while the size of the sectors represents the relative frequency of the wind directions

Full size image

Visual inspection of error diagnostic plots of the linear regression models revealed no serious deviations from the model assumptions. A subsequent analysis for errors of both predictor and response variables revealed neither any extreme outliers. The linear regression models fitted for the surface types water, mountain heathlands, and moss heathlands had each one or two outliers for the response variable, detected as studentized residuals with absolute values greater than three. Additionally, one observation of the explanatory variables in the linear regression model fitted for water was identified as an outlier (Cook’s distance of 1.89). However, since these outliers occurred during sudden cloud shifts, they were not removed from the data material. All local fitted cover-specific regression coefficients (\(\alpha \), \({\beta }_{2}\) and \({\beta }_{3}\)) (Table 2) were statistically significant for all linear regression models at 95% level. This indicated a significant linear relationship between both shortwave and longwave radiation and the net radiation for all surface types and substantiates the importance of correcting for local surface characteristics. Statistical assessment of the linear regression models, of both the fitting and prediction process, showed a high quality for all surface types (Table 3). The small RMSE and RMSEP substantiate the linear regression models’ good fit to the data and their strong predictive power for the different surface types. All linear regression models explained more than 99% of the variation in net radiation, also when performing the LOOCV (Table 3). Nevertheless, it should be recognized that coefficients of determination involve the sum of squares of deviations from the mean. A good fitting of large extreme observation of net radiation during midday might therefore inflate the coefficients of determination. However, deviations between surface-cover-specific observations of net radiation and predicted surface-cover-specific net radiation (Fig. 7 in Appendix 1) were found to be small, and the largest deviations occurred during midday with large incoming shortwave radiation. Based on these findings, the magnitude of the coefficients of determination reported here should therefore not be considered as unreliable. The errors introduced through the least squares regression in each individual cover-specific linear regression model were systematically smaller than the difference between tower-based net radiation and cover-specific net radiation for each individual surface type (Fig. 8 in Appendix 1). This illustrated the large accuracy of the linear regression models and their good ability to predict and detect cover-specific differences in net radiation during both shortwave- and longwave-dominated time periods. It is important to emphasize that if applying the methodology for surface-cover-correcting net radiation at another eddy-covariance site, the regression coefficients (Table 2) require to be rederived.

Table 2 Number of 30-min observations (N) of net radiation for each surface type, and regression coefficients for the linear regression models for each surface type
Full size table
Table 3 Root-mean-square (RMS) of net radiation for each surface type. Root-mean-square error (RMSE) and adjusted coefficient of determination (\({R}_{adj}^{2})\) for linear regression models fitted for each surface type, and root-mean-square prediction (RMSP), root-mean-square error prediction (RMSEP) and coefficient of determination for LOOCV (\({R}_{pred}^{2})\) for continuous predicted time series of net radiation for each surface type. Mean and standard deviation (Sd) of percentage contribution (\({p}_{i}\)) of the footprints, which were based on all observations
Full size table

The surface-cover-specific net radiation was not corrected for sloping terrain. To account for the full effect of the complex variation of both terrain slope and orientation, a detailed fine-spatial consideration of terrain slope and aspect is required. In open sloping terrain it has been found necessary to either provide slope-parallel measurements (Matzinger et al. 2003; Serrano-Ortiz et al. 2015) or correct horizontally measured radiation (Whiteman et al. 1989; Ramtvedt et al. 2021) to represent the actual radiation received by the surface. A failure to do so influences the energy balance by net radiation being out of phase with the turbulent energy (Hammerle et al. 2007; Hiller et al. 2008; Serrano-Ortiz et al. 2015; Wohlfahrt et al. 2016). This is particularly important for homogeneous, equally sloping surfaces. However, for heterogenous surfaces with varying sloping terrain with different orientation, the systematic effect on the net radiation and the albedo is less significant. At our study site, the terrain variability was expected to both reduce the diurnal time dependency of the albedo and minimize the importance of slope-correcting the net radiation. Likewise did the prevailing overcast conditions during the field campaign weaken the strong directional dependency of the incoming shortwave radiation. Both mechanisms substantiate that the simplification of expressing the albedo in the linear regression models with time constant regression coefficients did not critically impair the accuracy of the prediction modelling.

The ANCOVA revealed that there was a statistically significant difference between 13 of the 21 pairwise comparisons of surface types. Net radiation of lichen heathlands was found to be statistically different (at 95% level of significance) from all other surface types. During the daytime the lichen heathlands had systematically the lowest net radiation throughout the whole period (Fig. 4). This was reasonable considering their large albedo (see, e.g., Aartsma et al. 2020), as indicated by \({\beta }_{2}\) (which represents 1—albedo) in Table 2. During night-time, lichen heathlands were detected with the lowest net radiation among the vegetation surfaces (Fig. 4). Since lichens are spectrally similar to other vegetation for thermal infrared wavelengths (Abbot et al. 2013), this was likely explained by larger longwave emissions during night, which was reinforced by full-saturated lichen cover (Kershaw 2010) due to humid and rainy weather during the field survey. Net radiation of water and boulder fields was found to be statistically different (at 95% level of significance) from all other surface types except from each other and flood plains. Visual inspection (Fig. 4) revealed that water had substantially larger net radiation during daytime, and smaller net radiation during night-time compared to the other surface types. The daytime findings were reasonable due to the low albedo of water, as illustrated by \({\beta }_{2}\) in Table 2. Since the incoming longwave radiation from the atmosphere were assumed to be homogenous for the eddy-covariance site, the smaller net radiation of water during night-time was likely caused by larger emissivity and/or higher surface temperatures. According to Brewster (1992) the thermal emissivity (for surface temperatures of 300 K) is 0.92–0.96 for vegetation, ≈ 0.96 for water and 0.88–0.92 for rocks. In addition, the water’s large heat capacity dampens the diurnal temperature variation, thus increasing night-time surface temperatures when compared with the other surface types.

3.2 Assessment of Methodology

The small data coverage after discarding low quality observations of turbulent fluxes illustrates the challenging conditions of providing high-quality eddy-covariance measurements at the study site. Averaged diurnal fluxes for the period are shown in Fig. 6. During night-time the residual in the energy balance was systematically more negative when calculated with \({\tilde{R }}_{n}\) than when calculated with \({R}_{n}\). During daytime no difference was found (Fig. 6). The reported mean and standard deviation of \({p}_{i}\) in Table 3 were based on all observations. Water and fens were the surface types that contributed the most on average to \({\tilde{R }}_{n}\) throughout the period (Table 3). Three percent of the average footprint consisted of surface types not analyzed in the study (coloured grey in Fig. 3). A large variation in 30-min \({p}_{i}\) existed within each surface type, especially for water and fens (Table 3). This is reasonable due to the pronounced fine-spatial heterogeneity and the large differences in spatial extent of the 30-min footprints, as illustrated by the footprints in Fig. 3. However, despite the large difference in 30-min \({p}_{i}\), no clear difference in the EBC between wind directions could be detected (not shown here). Since the prevailing wind directions at the study site generally are north-west and south-east (approximately equally represented), as observed during this field survey, it should not be expected that larger surface variability itself effectively influence \({\tilde{R }}_{n}\) or the EBC. Insignificant differences (maximal difference of 2 percent for mean of \({p}_{i}\) and 1 percent for standard deviation of \({p}_{i}\)) between \({p}_{i}\)’s when computed with all observations and only high-quality observations were detected. This indicated that the large variability found in 30-min footprints was likely affected by the wind direction and to a lesser extent the conditions of wind speed, turbulence intensity, and atmospheric stability.

Fig. 6
figure 6

Averaged diurnal turbulent fluxes (H and LE), soil heat flux (\(G\)), one-point tower-based net radiation (\({R}_{n}\)), and surface-cover-corrected net radiation (\({\tilde{R }}_{n})\) for only high-quality observations. The shaded regions display the residuals in energy balance when calculated with \({R}_{n}\)(blue) and \({\tilde{R }}_{n}\)(orange)

Full size image

The OLS regression gave larger energy imbalance between the turbulent fluxes and available energy, than compared to the RMA regression (Table 4). The ranking of these closures was in good agreement with former results (Wilson et al. 2002; Li et al. 2005). The smaller energy imbalance found when allowing for diurnal effects was explained by the cancelling effect of energy surplus (storing) during daytime and energy deficit (releasing) during night-time. Especially during daytime, the momentary outgoing energy exchange of the turbulent fluxes was smaller than the net radiation. Lacking estimates of energy storage at a 30-min resolution produces a regression slope less than one. Therefore, longer averaging periods, such as a diurnal scale, or use of EBR as EBC assessment, should be considered appropriate to account for missing 30-min energy storage estimates (Leuning et al. 2012). No improvement was found in the EBC when correcting the mismatching flux footprints of turbulent fluxes and net radiation when assessed with the regression methods (Table 4). When assessed with the EBR, the energy balance improved only by 8% when correcting for surface heterogeneity. It is important to emphasize that comparisons between the regression assessments and the EBR and the effect of correcting for surface heterogeneity when assessed with the EBR should be made with caution, because differences could be devoted to the gap-filling itself. However, EBR based on data without gap-filling gave exactly the same closures as those reported in Table 4. The gap-filling was therefore excluded as cause to the energy balance improvement. Based on the findings, no or only a small effect seems to appear when providing surface-cover-corrected net radiation at heterogeneous eddy-covariance sites. Nevertheless, only a small dataset served as basis for this demonstration, and the results might not necessarily be unconditionally generalizable. Especially during daytime with cloudless conditions and large fractions of direct incoming shortwave radiation, larger differences in surface-cover-specific net radiation might give detectable differences between \({\tilde{R }}_{n}\) and \({R}_{n}\).

Table 4 Assessment of EBC with net radiation (\({R}_{n}\)) provided by the eddy-covariance tower and with surface-cover-corrected net radiation (\({\tilde{R }}_{n}\)), for ordinary least squares (OLS) regression, for reduced major axis (RMA) regression and energy balance ratio (\(\text{EBR}\)). The intercepts have unit W m−2
Full size table

At our eddy-covariance site it is likely other site-specific challenges than mismatching footprints that play considerable roles in explaining the energy imbalance. Both the hilly terrain with topographical slopes and the contrast between the land surface and the water surface in the river and the lake are potential mechanisms in generating advection as breezes and drainage flows. Substantially different thermal properties and differences in diurnal temperature rates between the water and the land surface questions the assumption of fully developed turbulence in the atmospheric surface layer (Blanford et al. 1991). These secondary circulations and internal boundary layers could decouple the eddy-covariance system from the surface (Aubinet 2008). This underestimates the turbulent fluxes. Since our eddy-covariance site is placed in a valley, cold and dense air tend to sink downslope by gravity (Stull 1988), which creates a shallow sublayer of cold air. The glacier located just a few kilometres south-west of the study site is expected to strengthen this effect due to katabatic flow from the glacier. It is known that the study site’s summer temperature is influenced by cold air from the glacier (Leinaas and Schumacher 2005). During wintertime, measurement of vertical temperature profiles in the surface layer at the study site indicated warmer layers at 10–30 m above ground level, which persisted even during well-mixed conditions (Pirk et al. 2019). This indicates possible biases in the eddy-covariance measurements. Additionally, it is expected that the heterogenous surface cover disturb the low-frequency covariance spectrum (Panin et al. 1998). It should be recognized that the likely presence of these mechanisms systematically lowers the turbulent transport, which might explain the large energy imbalance (Table 4) reported in this study. However, an unconventional eddy-covariance calculation using the ogive optimization method (Desjardins et al. 1989; Oncley et al. 1990; Sievers et al. 2015) to identify and separate out low-frequency flux contributions did not indicate an improvement in the EBC at our study site (Pirk et al. 2019). To understand the full effect of potential vertical divergence, and detect possibly biased data, large-eddy simulations with in situ eddy-covariance measurements and simultaneous observations at two heights of the fluxes are needed.

4 Conclusions

We have demonstrated a methodology for providing surface-cover-corrected net radiation matching the footprint of the turbulent fluxes at a heterogeneous eddy-covariance site. The key steps were to perform a high-resolution vegetation mapping for which surface-cover-specific in situ net radiation were measured, before continuous time series were predicted and area-weighted according to the turbulent footprint of the eddy-covariance tower. Statistical testing revealed significant differences in surface-cover-specific net radiation for 13 of 21 pairwise comparisons. This indicated the importance of accounting for surface heterogeneity and not relying on only the tower-based one-point measurement of net radiation. Because of cloudy weather at daytime during the field survey, it was not possible to provide and examine the full effect of differences in net radiation caused by differences in albedo between the surface types. Therefore, larger differences in net radiation than presented here are expected to occur between the surface types.

In contrast to our hypothesis, applying this methodology only gave a minor improvement of the energy imbalance. When assessed with the EBR the energy imbalance decreased by 8%, but no difference was detected when assessed with regression methods (OLS and RMA). The difference between the assessment methods was reasonably explained by lacking estimates of energy storage and diurnal sources and sinks. Based on our findings, providing surface-cover-corrected net radiation does not explain the energy imbalance at heterogenous eddy-covariance sites. However, only a small dataset served as basis for this demonstration. Further testing and application is required to understand the full effect of surface-cover-correcting mismatching footprints of turbulent fluxes and net radiation at heterogeneous eddy-covariance sites. Despite this, the study serves as a contribution to the knowledge about energy exchange processes in rarely investigated and reported ecosystems.