#### 2.1. Study Site and Instrumentation

The Ulu Peninsula, covering 552 km

^{2} (see [

26]), is the largest ice-free area in the AP region. JRI lies to the south-east of the tip of the AP (

Figure 1a), with a prevailing easterly flow often modified by the AP orography in relation to the southern and south-western winds [

27]. Consequently, the local climate on JRI is characterized by a 4–5 °C lower mean air temperature compared to the north-western AP [

28]. The study site was situated at 10 meters above sea level (m a.s.l.) on a Holocene marine terrace in the northern part of the Ulu Peninsula (

Figure 1). The Johann Gregor Mendel Station is located approximately 100 m to the north-east of the meteorological tower and a sonic anemometer, while the seashore is about 150 m to the north. The area is composed of dry sandy ground (about 10% moisture) underlain by continuous permafrost with an active layer of thickness between 50 and 65 cm [

23].

The measurements were carried out between 23 January and 21 March, 2018 (58 days). Details on the instruments used for measurements are in

Table 1 and they were all temporarily installed at one site. For the latent heat flux calculation [

29], we used the measurements of 2-m air temperature, relative humidity and surface ground temperature measurements. As only relative humidity was measured, 2-m specific humidity was calculated with the use of air temperature and atmospheric pressure data. Furthermore, vertical kinematic eddy heat flux, Monin–Obukhov length and virtual temperature from sonic anemometer were also included in the latent heat flux calculation. The ground surface temperature was defined as the ground temperature at a depth of 2 cm.

#### 2.2. Data Processing and Analysis

The surface energy budget of an ice-free surface can be written as:

where Q

_{net} is the net radiation, Q

_{G} is the ground heat flux, Q

_{H} is the sensible heat flux, Q

_{E} is the latent heat flux and C is the residual. The Q

_{net} has the following components:

where Q

_{Snet} is the net shortwave flux consisting of incoming (Q

_{S↓}) and outgoing (Q

_{S↑}) shortwave radiation, while Q

_{Lnet} is the net longwave flux, which can be divided into incoming (Q

_{L↓}) and outgoing (Q

_{L↑}) longwave radiation.

All of the measured data were processed as 30-min averages. The Q

_{S↓} measurements were further corrected with the use of a more precise CM11 pyranometer (Kipp and Zonen). Simple linear regression was used for adding the relative deviation to CNR4 measurement [

30]. Surface albedo was calculated daily as the ratio of Q

_{S↑} and Q

_{S↓} in the 3 h closest to the maximum of Q

_{S↓} in order to avoid measurement unreliability under large solar zenith angles and low flux magnitudes.

Ground heat flux measurement was corrected using the calorimetric method [

31,

32] as follows:

where G

_{i–1} is the heat flux density at the top of a layer, here at a depth of 2 cm, δz

_{i} is the layer thickness (3 cm), C

_{i} is the volumetric heat capacity for the layer (1.17 MJ m

^{−3} K

^{−1} defined for the area by Hrbáček et al. [

23]), ∂T

_{i}/∂t is the rate of a 30-min change in the mean layer temperature and G

_{i} is the heat flux density at the bottom of the layer at a depth of 5 cm.

The sensible heat flux was calculated from sonic anemometer using the eddy covariance technique (e.g., [

33]). However, it needs to be noted that when the term “sensible heat flux” is used in this study, it is the buoyancy flux, which was calculated instead of true sensible heat flux due to the absence of fast-response moisture measurements [

34]. Sensible heat fluxes of 30-min average were calculated only if more than 54% of the measurement were available. Consequently, <2% of sensible heat flux data were missing. Latent heat was calculated as in Langer et al. [

29]. Briefly, latent heat flux was estimated as:

where ρ

_{air} is the density of air, L

_{lg} is the latent heat of vaporization of water, r

_{a} is atmospheric resistance, q is specific humidity, z

_{m} is measurement height and z

_{0} is roughness length. Roughness length was originally chosen to be 10

^{−3} m over bare ground and 10

^{−4} m over a snow-covered surface [

35,

36], but due to unrealistically high values of latent heat flux, they were both lowered by one order to 10

^{−4} and 10

^{−5} m, respectively. The atmospheric resistance r

_{a} was calculated as in formula D4 in Langer et al. 2011 [

29]. Specific humidity at the surface q(z

_{0}) was derived from the surface temperature using the Magnus formula [

37].

Latent heat flux measurements were only available for 92% of the study period. For all surface energy budget components, daily and weekly means were calculated if <20% were missing, hence the loss of daily mean fluxes on 8 February and 12–14 March, 2018 (

Figure 2).

The data were mostly analysed and compared as daily means; 30-min averages were only examined in

Table A1 and day-long case studies in order to determine the influence of snow and cloudiness on the surface energy budget (

Section 3.4). For the correlation analysis (

Section 3.3), Spearman’s rank correlation coefficient was used due to non-normal distributions of some data. Simple linear regression was employed for calculating regression function in

Section 3.3. The following convention was used: the net radiation components were considered positive if energy transport was directed towards the ground surface; the turbulent and ground heat fluxes were denoted as positive if energy transport was directed away from the ground surface.

For the case studies, it was necessary to determine a cloudy and a clear-sky day. We only had three categories: clear sky, cloudy and overcast. Cloudiness was derived from a combination of hemispherical photos and Q_{S↓} measurements. We calculated the mean daily ratio of Q_{S↓} to extra-terrestrial radiation (simply approximated from solar constant and solar zenith angle) in hourly intervals. We then estimated that if the daily mean ratio was above 0.6, the sky was (mostly) clear. If the ratio was below 0.4, the day was overcast. Ratios between 0.4 and 0.6 meant a cloudy day.

#### 2.3. Measurement Accuracy

The residual term C in our study was 36.1 ± 104.0 W m

^{−2} (

Table A1), indicating that errors prevented surface energy budget closure. Unlike the measurements from KGI, where C was not determined but was contained within Q

_{H} [

17] or Q

_{G} [

13], we tried to limit the uncertainty by calculating unmeasured Q

_{E}. Non-zero C has been almost unavoidable in surface energy budget studies from polar regions [

29,

38], even though errors should have been cancelled out unless there was a systematic error. In our case, 35% of Q

_{net} was allocated to neither Q

_{G} nor the turbulent fluxes due to a combination of multiple measurement errors.

Firstly, the measurement accuracy of the CNR4 for longwave radiation of ±10% might have increased the values of Q

_{net} (

Table 2). Moreover, Westermann et al. [

38] pointed out that uncertainty about albedo under conditions with patches of snow and bare-ground could be up to 0.1, which would further lower or raise Q

_{net} by 3.4%.

Secondly, the eddy covariance technique requires stationarity. However, complex orography around the measuring site was prone to the development of local circulation patterns [

36]. That could have led to systematic underestimation of Q

_{H} up to 25% [

39], while the use of buoyancy flux as Q

_{H} means that another 10% might need to be added to Q

_{H} [

34].

The error in Q

_{E} was difficult to assess because no comparable measurements existed. Nevertheless, Q

_{E} calculating depended strongly on the correct choice of z

_{0}. In this study, using the values for snow-covered surfaces and bare ground according to current knowledge [

35,

36] could lead to an increase of Q

_{E} up to 33%. On the other hand, as both Prosek et al. [

17] and Choi et al. [

13] reported mean Q

_{E} ~ Q

_{H}, it is likely that the error is smaller. Even though the maximum uncertainty in Q

_{G} estimation caused by instrument accuracy could have been up to ±20%, the mean value of Q

_{G} was below 1 W m

^{−2}. As a result, the error in Q

_{G} would not contribute significantly to C.

Finally, individual energy budget components were measured on different spatial scales. Q

_{G} was derived from a point measurement and radiation sensors had a footprint area of several square meters. For the sonic anemometer at a height of 2 m, 60% of the flux might have come from an area with a radius of 75 m or larger [

40]. During high wind speeds, the turbulent eddies (and derived energy fluxes) measured by the sonic anemometer were coming from the Bohemian Stream Catchment (see

Figure 1 and

Section 3.1), where the ground was waterlogged or moister than close to the meteorological tower. Consequently, more energy was allocated to Q

_{E} in the area from which the turbulence was coming, influencing the residual term C.

Considering the above-mentioned errors, the magnitude of C could be explained by the uncertainty inherited from methods of measurements, the landscape diversity around the measuring site and the choice of parameters in the calculation. Essential for future experiment concept and better accuracy of eddy-covariance method could be two-level humidity measurement or krypton hygrometer for measuring water vapour fluctuations with a high-frequency sampling rate.