Georg Lackner
email: georg.lackner@mailbox.org

Daniel F Nadeau

Florent Domine

Annie-Claude Parent

Gonzalo Leonardini

Aaron Boone

François Anctil

Vincent Fortin

The effect of soil on the summertime surface energy budget of a humid Subarctic tundra in northern Quebec, Canada

Rising temperatures in the southern Arctic region are leading to shrub expansion and permafrost degradation. The objective of this study is to analyse the surface energy budget (SEB) of a subarctic shrub tundra site that is subject to these changes, on the east coast of Hudson Bay in eastern Canada.

We focus on the turbulent heat fluxes, as they have been poorly quantified in this region. This study is based on data collected by a flux tower using the eddy covariance approach and focused on snow-free periods. Furthermore, we compare our results with those from 6 Fluxnet sites in the Arctic region and analyse the performance of two land surface models, SVS and ISBA, in simulating soil moisture and turbulent heat fluxes. We found that 22.4% of the net radiation was converted into latent heat flux at our site, 33% was used for sensible heat flux and about 15% for ground heat flux. These results were surprising considering our site was by far the wettest site among those studied, and most of the net radiation at the other Arctic sites was consumed by the latent heat flux. We attribute this behaviour to the high hydraulic conductivity of the soil (littoral and intertidal sediments), typical of what is found in the coastal regions of the eastern Canadian Arctic. Land surface models overestimated the surface water content of those soils but were able to accurately simulate the turbulent heat flux, particularly the sensible heat flux and to a lesser extent, the latent heat flux.

3

Significance statement. Isostatic uplift after the last deglaciation led to the emergence of beaches, which represent a large area in the Canadian Arctic. We studied the surface energy budget of such a beach emerged 6000 years ago in northeastern Canada. Results suggest that evaporation is up to 20% less than at previously studied sites which is explained by sandy soils which retain little moisture despite abundant precipitation. Deployed numerical models showed difficulties in simulating the soil conditions but proved successful in simulating the surface energy budget after manual adjustments of the soil conditions. These simulation difficulties probably apply to other parts of the Canadian Arctic, possibly leading to errors in meteorological and climate forecasting.

Introduction

The increase in global air temperature, which is twice as pronounced in polar regions than in all other parts of the world [START_REF] Chylek | Arctic air temperature change amplification and the atlantic multidecadal oscillation[END_REF], is changing the distribution of vegetation zones (Myers-Smith and Hik 2018). This change is particularly notable in the forest-tundra ecotone [START_REF] Payette | The Subarctic forest-tundra: The structure of a biome in a changing climate[END_REF]) at the interface between the boreal forest and the Arctic shrub tundra. This region is characterized by a mosaic of forest patches, usually restricted to humid, wind-sheltered areas, surrounded by shrub tundra that occupies the well-drained and wind-exposed sites. Rising temperatures favor the progressive expansion of these forest patches, which lead to the northward expansion of the boreal forest [START_REF] Tape | The evidence for shrub expansion in northern alaska and the pan-arctic[END_REF]).

The disappearance of shrub tundra at the southern Arctic border is leading to changes on multiple levels. Heat and moisture exchanges at the interface between the land surface and the atmosphere are being altered. As such, the surface energy budget (SEB) as a whole is also changing due to shifts in leaf area index, land surface albedo, and evapotranspiration. These changes then create a feedback affecting the climate at multiple scales by modifying ecosystem processes and weather patterns [START_REF] Baldocchi | How will land use affect air temperature in the surface boundary layer? Lessons learned from a comparative study on the energy balance of an oak savanna and annual grassland in california, usa[END_REF][START_REF] Pielke | Land use/land cover changes and climate: modeling analysis and observational evidence[END_REF]). Due to the circumpolar extent of the forest-tundra ecotone, its potential impact on the climate demands our attention.

Moreover, the surface energy budget influences the fate of the permafrost as it controls how much heat is stored in or released from the ground. Permafrost in this region has become discontinuous to sporadic and is degrading quickly [START_REF] Lemieux | Groundwater occurrence in cold environments: Examples from Nunavik, Canada[END_REF]. In order to model carbon emissions due to permafrost thaw [START_REF] Schuur | Climate change and the permafrost carbon feedback[END_REF], information on the partitioning of the net radiation into sensible, latent, and soil heat fluxes is key. Thus, an improved understanding of the energy balance can ultimately lead to more accurate predictions of permafrost thawing and the resulting greenhouse gas emissions.

Due to the harsh conditions and limited accessibility in the Subarctic region, few in-situ observations have been conducted. Available datasets come from sites located in central and western North America, Russia and Europe [START_REF] Lafleur | Energy balance differences and hydrologic impacts across the northern treeline[END_REF][START_REF] Eugster | Land-atmosphere energy exchange in Arctic tundra and boreal forest: Available data and feedbacks to climate[END_REF][START_REF] Langer | The surface energy balance of a polygonal tundra site in northern Siberia -part 1: Spring to fall[END_REF][START_REF] Stiegler | Tundra permafrost thaw causes significant shifts in energy partitioning[END_REF]). The topography of these sites often consists of plains or gentle slopes, which tend to reduce drainage efficiency and provide wetter soils. The latent heat flux at these sites tends to dominate the SEB over the other terms. For example, [START_REF] Eugster | Land-atmosphere energy exchange in Arctic tundra and boreal forest: Available data and feedbacks to climate[END_REF] reported latent heat flux to net radiation ratios in the range of 0.44 to 0.72 at several low Arctic sites in Alaska, and values between 0.48 to 0.65 for treeline shrub tundra across central Canada, Alaska, and Russia.

One aspect that has so far been poorly documented is the impact of soil properties on the SEB of subarctic sites, as the vegetation tended to be considered the dominant factor on the issue. In the eastern region of the Canadian Subarctic, for instance, soils are often composed of coarse quaternary deposits such as glacial till and marine sediments [START_REF] Lemieux | Groundwater occurrence in cold environments: Examples from Nunavik, Canada[END_REF]). These soils are usually very well-drained even on small inclines, limiting water availability and making this region distinct from those of previous studies.

In addition to the quantification of surface fluxes, a pressing issue is whether current land surface models (LSMs) such as ISBA [START_REF] Noilhan | A simple parameterization of land surface processes for meteorological models[END_REF] and SVS [START_REF] Husain | The multibudget soil, vegetation, and snow (SVS) scheme for land surface parameterization: Offline warm season evaluation[END_REF][START_REF] Alavi | Warm season evaluation of soil moisture prediction in the soil, vegetation, and snow (SVS) scheme[END_REF]) are able to adequately simulate the SEB in Subarctic regions. LSMs are used in a variety of applications, such as numerical weather prediction models and terrestrial ecosystem models.

Although the surface of the forest-tundra environment is snow-covered for most of the year, the adequate modeling of the SEB during snow-free period is equally important because this period regulates the growth of vegetation and to a larger extent, the water cycle. Studies have highlighted difficulties in modeling dry sites, stating that the models capture energy limitations better than water limitations [START_REF] Engstrom | Modeling evapotranspiration in Arctic coastal plain ecosystems using a modified BIOME-BGC model[END_REF]), thus leading to increased errors for simulated latent heat fluxes.

In this study, we use eddy covariance data to analyse the Arctic surface energy balance and examine the performance of two LSMs: SVS and ISBA. SVS will soon replace the old version of ISBA and be used as an operational weather-forecasting model in Canada. This LSM has been proven to be competent at simulating latent and sensible heat fluxes at six sites in arid, Mediterranean, and tropical climates [START_REF] Leonardini | Evaluation of the soil, vegetation, and snow (svs) land surface model for the simulation of surface energy fluxes and soil moisture under snow-free conditions[END_REF]). This makes it particularly interesting to compare whether SVS performs better in simulating northern Canadian regions than the current version of ISBA. This study presents a comprehensive analysis of observed and modeled surface energy fluxes during three consecutive snow-free periods in a forest-tundra ecotone site in northern Quebec in eastern Canada. The objectives of this study are to (i) use eddy covariance observations to identify key factors that influence the energy partitioning at the surface, (ii) determine the interannual variability of the partitioning, (iii) contrast the observations with those from other Subarctic sites across the globe, and (iv) assess the capability of the current LSMs, ISBA and SVS, to simulate surface heat fluxes in this environment.

Methods

a. Study Site

The study site is located in the Tasiapik valley close to the community of Umiujaq, on the eastern shore of Hudson Bay in Quebec, Canada (56 • 33'31"N, 76 • 28'56"W; Umiujaq in Figure 1a). The valley is 4.5 km long and 1.3 km wide and surrounded by cuestas (steep ridges) with a maximum height of nearly 400 m. The lower part of the valley borders the Tasiujaq Lake at an elevation of 0 m and is connected to Hudson Bay by an inlet. The vegetation cover in the valley is typical of a forest-tundra ecotone, with shrub and some lichen tundra in its upper region and forest tundra in the lower region. While the dwarf birch Betula glandulosa is the dominant vegetation type in the shrub tundra, the prevailing tree species in the forest-tundra are black and white spruces (Picea mariana [Mill.] BSP and Picea glauca [Moench] Voss) [START_REF] Paradis | Greater effect of increasing shrub height on winter versus summer soil temperature[END_REF]). Following Gagnon et al. (2019), vegetation in the vicinity of the flux tower could be classified as one of three types: pure lichen cover, small shrubs with lichen understory and a low LAI (presumably <0.5) and small shrubs with moss understory. Lichens covers 20 to 30% of the surface while small shrubs with differing understory make up the rest (See Figure A3).

In the valley, permafrost is discontinuous to sporadic [START_REF] Lemieux | Groundwater occurrence in cold environments: Examples from Nunavik, Canada[END_REF]) with rapidly degrading permafrost mounds that are ≈ 100 m in diameter [START_REF] Fortier | Fast permafrost degradation near Umiujaq in Nunavik (Canada) since 1957 assessed from time-lapse aerial and satellite photographs[END_REF]. Close to the study site, permafrost is only present to the south at a distance of about 20 m. The active layer thickness is about 1.5 to 2 m. This active layer depth together with the small inclination of a few degrees of the surface presumably provides a rather undisturbed drainage of precipitation.

There are no complete long-term meteorological records available for Umiujaq. Instead, recent observations from the research station in the Tasiapik valley were used to estimate climatological conditions starting in 2012 for air temperature and from 2015 for precipitation. Umiujaq exhibits a Subarctic climate with a mean annual temperature of -4.0 • C. Their winters are long and cold and their summers are cool. There is usually continuous snow cover from late October to early June with a maximum thickness between 0.6 m and 1.0 m. Monthly temperatures and precipitation are shown in Figure 2.

Weather patterns are strongly influenced by Hudson Bay to the west and whether or not it is covered with ice. Hudson Bay usually remains frozen until mid-June, and then stores heat in summer that is released in the fall. This results in relatively lower air temperatures in summer and higher air temperatures in the fall. The average total annual precipitation is about 800 -1000 mm, out of which about 50% falls as snow and the precipitation cycle shows a strong correlation with the freezing of Hudson Bay. After the freeze-up in mid to late December, precipitation remains low from January to May/June. Precipitation subsequently increases and peaks in late summer and autumn. Fog is common throughout the summer [START_REF] Robichaud | The Weather of Atlantic Canada and Eastern Quebec[END_REF]. Advection fog forms when warmer air moves over the cold Hudson Bay, mostly in July and August.

The geology of the valley is thoroughly described in [START_REF] Lemieux | Groundwater dynamics within a watershed in the discontinuous permafrost zone near Umiujaq (Nunavik, Canada)[END_REF]. The bedrock is covered by a 10-to 50-m thick layer of moraine deposits and fluvioglacial sediments. Above the fluvioglacial sediments, there is a layer of marine sediments (mainly silt) that is ≈ 30-m thick, but it does not cover the entire valley. The topsoil layer is sand (littoral and intertidal sediments, see

Figure 3), with some limited exceptions. [START_REF] Gagnon | The carbon sink due to shrub growth on Arctic tundra: a case study in a carbon-poor soil in eastern Canada[END_REF] found that the soil covered with lichen is made up almost entirely of sand (90% for the first 5 cm and almost 100% for deeper layers up to 30 cm). They reported a soil organic content of about 1.4 kg m -2 under lichen and 4.3 kg m -2 for medium height shrubs tundra with lichen understory. Most of the soil organic content is found in the first 2 cm for the lichen covered areas and 5 cm for medium shrub covered areas.

b. Instrumental Setup

In June 2017, a 10-m tower was equipped with a sonic anemometer and a CO 2 /H 2 O gas analyzer that was installed 4.2 m above the ground (IRGASON, Campbell Scientific, USA), facing eastwards towards a plateau with a small slope of about 3 • (Figure 1 b). The tower holding the instrument and a steep ridge are to the west, in the opposite direction.

The tower also features general meteorological sensors, such as temperature and relative humidity (model HMP45, Vaisala, Finland) and wind direction/speed (model 05103, R.M. Young, USA) probes. Some 10 m west of the flux tower, another meteorological station was installed in 2012, which also measures basic environmental variables such as air temperature and wind speed. It is equipped with a 4-component radiometer (CNR4, Kipp and Zonen, The Netherlands). In the proximity of the two stations, soil temperature and water content (5TM, Meter Group, USA) are measured at several locations and under various vegetation types directly at the tower and several meters distance, at 5 depths from 5 cm to 50 cm, in order to capture spatial variability.

Precipitation is also measured next to the tower with a T200B precipitation gauge (GEONOR, USA) equipped with a single Alter shield. Raw precipitation observations were corrected for wind-induced undercatch following [START_REF] Kochendorfer | Testing and development of transfer functions for weighing precipitation gauges in WMO-SPICE[END_REF]. Flux data was recorded using a CR3000 datalogger while the other data was recorded using a CR1000 datalogger (both Campbell Scientific, USA). The site was visited twice a year for maintenance and data collection.

c. Data Processing

Raw 10-Hz eddy-covariance data were processed using the EddyPro ® (version 7.0.3; Li-COR Biosciences, USA) software package. This software includes multiple corrections and QA/QC procedures such as turbulent fluctuation detrending based on a running mean, covariance maximization, density fluctuation compensation (Webb et al. 1980), and analytic correction of high-pass and low-pass filtering effects [START_REF] Moncrieff | A system to measure surface fluxes of momentum, sensible heat, water vapour and carbon dioxide[END_REF]). Turbulence data were rotated into a streamline coordinate system using the sector-wise planar-fit method [START_REF] Wilczak | Sonic anemometer tilt correction algorithms[END_REF]). The random uncertainty quantification approach from [START_REF] Finkelstein | Sampling error in eddy correlation flux measurements[END_REF] was used to detect outliers, spikes, and other artifacts. The 0-1-2 quality scheme from [START_REF] Mauder | A strategy for quality and uncertainty assessment of long-term eddy-covariance measurements[END_REF] was applied to the data, and only segments that were flagged as 0 and 1 were kept in the study. These values correspond to good-and medium-quality observations, which are both suitable for the analysis of seasonal energy budgets.

Ground heat flux was computed using a time series of a measured soil temperature profile ( 1 , ), ( 2 , ) and ( 3 , ) at three depths 1 < 2 < 3 and the one-dimensional heat equation

( , ) = 2 ( , ) 2 , (1)
where is the thermal diffusivity [START_REF] Westermann | The annual surface energy budget of a high-arctic permafrost site on svalbard, norway[END_REF]). As the thermal diffusivity is defined by the thermal conductivity and the heat capacity ℎ by = / ℎ we calculated through and ℎ using empirical equations which are based on the soil density and water content [START_REF] Chen | Thermal conductivity of sands[END_REF][START_REF] Bittelli | [END_REF]. Then, the heat equation was solved numerically in Python using two of the temperature sensors in the soil ( ( 1 , ) and ( 3 , )) as Dirichlet boundary conditions. The result is a temperature distribution in the considered spatial domain. It was checked against the third soil temperature measurement ( 2 , ) for correctness (R 2 = 0.99, RMSE = 0.21 • C). Then, the temperature gradient at the upper boundary 1 was calculated. The ground heat flux at depth 1 is then given by the Fourier equation:

1 = ( , ) = 1 (2)
The heat storage change in the soil is accounted for by taking the change of the mean temperature of the soil above 1 . The mean temperature is obtained by extrapolating the temperature distribution given by the heat equation. The final ground heat flux is than given by = 1 + . This method was chosen as measurements under different vegetation types were available making it possible to analyze the spatial variability. Furthermore, problems have been reported of alternative methods such as using heat flux plates (Ochsner et al. 2006).

Post-processing was also necessary to sort out the remaining outliers and to fill data gaps. This was achieved using the program PyFluxPro [START_REF] Isaac | Ozflux data: network integration from collection to curation[END_REF]). This program has six processing levels, uses EddyPro output files as input and produces a continuous time series for all fluxes. For the first three processing levels, data are read, quality controlled, and finally auxiliary measurements are merged when gaps are present. The quality control includes (i) range checks based on userdefined limits, (ii) spike detection, (iii) manual removal for specific dates and (iv) data rejection based on other variables. For the data rejection, CO 2 and H 2 O signal strengths from the IRGA as well as internal error codes from both the sonic anemometer and the IRGA were used to remove erroneous flux data. For the fourth processing level, meteorological variables were gap filled with ERA5 data, a reanalysis product from the European Centre for Medium-Range Weather Forecast.

These data provide hourly estimates for various meteorological and soil variables beginning in 1979 at a spatial resolution of 30 km (https://www.ecmwf.int/en/forecasts/datasets/ reanalysis-datasets/era5). Each variable was bias-corrected using a linear fit between ERA5 and flux tower observations during periods when both were available.

Finally, the fluxes were gap-filled using a self-organizing linear output map (SOLO), a type of artificial neural network (ANN) (see [START_REF] Hsu | Self-organizing linear output map (SOLO): An artificial neural network suitable for hydrologic modeling and analysis[END_REF] and [START_REF] Abramowitz | Towards a benchmark for land surface models[END_REF]). SOLO requires a set of environmental drivers such as air temperature, radiation and humidity, as well as the fluxes as inputs. SOLO first constructs relationships between the environmental drivers by applying an ANN-equivalent of a principal component analysis. It then uses an ANN-equivalent of a multiple linear regression to make connections between the drivers and the fluxes. ANN together with marginal distribution sampling (MDS, [START_REF] Reichstein | On the separation of net ecosystem exchange into assimilation and ecosystem respiration: review and improved algorithm[END_REF]) were shown to be the best choices for gap-filling flux data [START_REF] Moffat | Comprehensive comparison of gap-filling techniques for eddy covariance net carbon fluxes[END_REF]).

On average, 56 This study focused on the summer period, defined here as snow-free periods. It started on the date that snow disappeared and ended when the first snowfall occurred in order to avoid any effects associated with the presence of snow on the ground. In order to facilitate inter-annual comparisons, we applied these criteria to the three measurement years and defined common start and end dates:

July 1 to September 5 ℎ (66 days).

d. Comparison Sites

To give the energy balance measurements from the Umiujaq site context, data from six other flux tower sites across the Arctic were selected from the FLUXNET [START_REF] Pastorello | The FLUXNET2015 dataset and the ONEFlux processing pipeline for eddy covariance data[END_REF]) database.

These sites were located in Canada, the United States, and Russia and were selected within latitudes ranging from 56 • N to 70 • N (Figure 4). FLUXNET is a global network of micrometeorological tower sites that collects eddy covariance observations. We have chosen sites that show some similarities concerning the vegetation (shrubland) and are not considered as permanent wetland.

The comparison sites cover a variety of biomes located in Arctic and Subarctic regions (Figure 4). Their characteristics are detailed in Table 1. Mean annual temperatures range from -14. Vegetation covers at all the sites, except for RU-COK, is classified as open shrubland. Those sites are covered by all evergreen or deciduous woody vegetation less than 2 m tall at coverage between 10% and 60%. The vegetation at RU-COK is closed shrubland, given that the canopy cover exceeds 60%.

Data processing was similar to that described in section 2c. Quality assurance and quality control methods follow those in [START_REF] Pastorello | Observational data patterns for time series data quality assessment[END_REF], where meteorological variables were gap-filled using ERA-Interim reanalysis data [START_REF] Vuichard | Filling the gaps in meteorological continuous data measured at FLUXNET sites with ERA-Interim reanalysis[END_REF], and flux data were gap-filled with marginal distribution sampling [START_REF] Reichstein | On the separation of net ecosystem exchange into assimilation and ecosystem respiration: review and improved algorithm[END_REF]).

As for most of the eddy covariance sites worldwide, all studied sites experienced non-closure of the energy budget. Since not all terms of the energy budget were available for each site and given the ongoing debate on how to best distribute the residual energy budget [START_REF] De Roo | A semi-empirical model of the energy balance closure in the surface layer[END_REF], no corrections were applied to any of the energy budgets.

e. Model Description

To simulate the SEB, two LSMs were used on offline mode: ISBA and SVS. ISBA was included in SURFEX v8. Both models require meteorological data as inputs. These data can be obtained from an atmospheric model when operated in coupled mode, or from observations in stand-alone mode (offline).

In this study, the stand-alone option was used to ensure the quality of the driving data and to avoid errors that may be introduced by the atmospheric model itself.

1) ISBA-MEB ISBA [START_REF] Noilhan | A simple parameterization of land surface processes for meteorological models[END_REF] simulates all water and energy exchanges between the atmosphere, vegetation, and soil. The water and energy budgets of soil are computed by solving the one-dimensional Fourier law and a mixed-form of the Richards equation explicitly [START_REF] Boone | The influence of the inclusion of soil freezing on simulations by a soil-vegetation-atmosphere transfer scheme[END_REF][START_REF] Decharme | Local evaluation of the interaction between soil biosphere atmosphere soil multilayer diffusion scheme using four pedotransfer functions[END_REF]). The soil is therefore divided into layers of increasing depth down to 12 m. Each layer has specific hydraulic and thermal properties, which are either calculated internally based on the soil texture or forced externally. The phase change between the solid and liquid phases of water in the soil is computed using the Gibbs free-energy method [START_REF] Boone | The influence of the inclusion of soil freezing on simulations by a soil-vegetation-atmosphere transfer scheme[END_REF][START_REF] Decharme | Impacts of snow and organic soils parameterization on northern Eurasian soil temperature profiles simulated by the ISBA land surface model[END_REF]). This method calculates the temperature that governs the phase change based on the soil texture. During any phase-change, the total amount of water in each layer remains unchanged. Therefore when the ice content increases, the water content decreases by the same amount, and vice versa.

Vegetation characteristics can be automatically selected from 19 available categories using the site coordinates and ECOclimap (https://opensource.umr-cnrm.fr/projects/ ecoclimap-sg/wiki), or they can be specified by the user (vegetation type, LAI, canopy height and albedo).

In order to specifically include vegetation in ISBA, Boone et al. (2017) recently added what is called the multi-energy balance (MEB) parametrization. MEB largely follows the representation of vegetation detailed in the Rossby Centre Regional Climate Model [START_REF] Samuelsson | The Rossby Centre Regional Climate Model RCA3: Model description and performance[END_REF]).

The most important additions (relative to the standard version of ISBA) are the presence of canopy turbulence, a new radiation transfer scheme, and a litter layer. Within MEB, the sensible and latent heat fluxes between multiple components are calculated separately and then combined to form the final heat flux. These components are the vegetation, the bare ground, a ground-based snowpack, and the canopy air space.

2) SVS

The Soil, Vegetation and Snow scheme (SVS) [START_REF] Husain | The multibudget soil, vegetation, and snow (SVS) scheme for land surface parameterization: Offline warm season evaluation[END_REF][START_REF] Alavi | Warm season evaluation of soil moisture prediction in the soil, vegetation, and snow (SVS) scheme[END_REF]) is partly based on an early version of ISBA (Bélair et al. 2003a,b), in which a force-restore method is used for the energy budget, as well as a single layer interaction between the canopy and the atmosphere.

A notable improvement of this model is the different land surface tiling and vertical layering, as well as having separate energy budgets for bare soil and (low and high) vegetation tiles. Another difference between the two LSMs concerns soil hydrology. The earlier version of ISBA monitored soil moisture at a superficial layer and the rooting-depth layer. In SVS, a certain soil water content is attributed to each of the N soil layers and its evolution is modeled based on Darcian flow.

However, in SVS soil temperature is only defined at the soil surface and the associated variation is determined using the single-layer approach mentioned above.

3) M

Both ISBA and SVS require the same input data: air temperature, specific humidity, wind speed, incoming shortwave and longwave radiation, atmospheric pressure, and (solid and liquid) precipitation rates. Observational data for these variables are available from 2012, except for atmospheric pressure, for which the data start from June 2017 and precipitation where measurements started in 2015. Before this date, ERA5 atmospheric pressure and precipitation was used. A simple linear regression was conducted between the measured data and the ERA5 data in order to correct for the atmospheric pressure data from ERA5 and make it consistent with in situ observations. A fixed threshold of 0 • C was applied to separate the solid and liquid phases of precipitation. This simple threshold was implemented because we were only studying the summer period, and thus a more elaborate phase partitioning strategy was unnecessary.

4) M S

Both models were set up in a similar fashion to maximize the comparability of the results.

However, we have opted to use the soil discretization that is most commonly used for each model and which happen to be different from each other. For ISBA, the soil column has a depth of 12 m and is divided into 20 layers of increasing depth. The soil column for SVS is 10-m-deep and has 7 layers, also of increasing depth.

Models were initialized with a spin-up to ensure an equilibrium of soil moisture and temperature.

We used 5 complete years of observations (2012-2017) as the spin-up, which was sufficient for both models to reach an equilibrium. We analyzed the melt and onset of the snow cover in the model to make sure that snow did not interfere with our study period. We used a time step of 1 hour for both models.

The soil composition that was used for both models was 95% sand and 5% silt [START_REF] Gagnon | The carbon sink due to shrub growth on Arctic tundra: a case study in a carbon-poor soil in eastern Canada[END_REF]). Hydraulic and thermal parameters were for the most part calculated by the models themselves, based on the soil composition provided. A modified version of both ISBA and SVS models were also tested where hydraulic parameters were adjusted. The adjusted parameters were the saturated soil water content and the field capacity, both of which were based on soil water content observations at several depths. For both parameters, a profile that decreased with increasing soil depth was selected. The saturated soil water capacity was set to 0.3 m 3 m -3 at the surface and 0.08 m 3 m -3 at depths of 1 m and more. The field capacity at the surface was 0.18 m 3 m -3 in both models and was decreased to a minimal value of 0.05 m 3 m -3 at soil layers deeper than 1 m. For comparison, we also ran both models without changing any of the parameters. This is referred to as the default version of the model.

5) M E

Only a subset of the whole dataset was used for model evaluation. We selected data with energy budget closures (EBC) between 0.8 and 1.2, where EBC is defined as:

= + + . (3)
In equation 3, and are the turbulent fluxes of sensible and latent heat, is the ground heat flux and is the net radiation. We applied this procedure to minimize possible measurement errors and assure comparability between the models, as their energy budget closures are equal to 1 by definition. In order to evaluate and visualize the performance of ISBA and SVS in simulating and , we used Taylor diagrams [START_REF] Taylor | Summarizing multiple aspects of model performance in a single diagram[END_REF], which can effectively illustrate several performance metrics on the same plot. The normalized standard deviation for a dataset from its mean is given by the radial distance from the origin at (0|0). The angle between the point of interest, and the abscissa defines the correlation of the model with the observations. The normalized RMSD of the model is shown by the distance to the point indicated as Observation at (1|0). The standard deviation and RMSD are both normalized using the standard deviation of the observations.

Results

a. General Conditions in Umiujaq

Table 2 shows the conditions during the three summers of interest at the Umiujaq site. These conditions varied considerably in terms of air and soil temperatures, with summer 2019 being the warmest and summer 2018 being the coldest. Also, precipitation differed considerably, although there was no dry period in neither of the three summer, with a total of 237 mm in 2017 and 369 mm in 2019. Another marked difference, particularly in the beginning of summer, is the amount of fallen snow from the preceding winter. Summer 2018 was delayed by more than two weeks due to an exceptionally thick snow cover (≈ 1.4 m) from the previous winter. This was not the case for the other two summers, as the maximum snow depths from those winters were more typical (0.6-1.0 m).

Studying the one site over three summers allows us to investigate an interesting range of interannual variability, with total seasonal precipitation varying by up to 132 mm (factor of 1.56) and mean air temperatures varying by up to 2 • C. In general, we can classify the summer of 2017 as dry, the summer of 2018 as cold, and the summer of 2019 as warm and humid. the applied EBC criterion. The average Bowen ratio BR ( / ) fluctuated between 1.4 and 1.54 for and EBC ranges between 0.2 to 1.5 and 0.8 to 1.2. The energy budget closure was best in 2019, with a mean closure of 76.7%. Indeed, in 2019, there were more frequent winds from areas less affected by topography, and therefore more favorable to the use of the eddy covariance technique.

b. Evolution of the Surface Energy Budget

Although conditions varied considerably between the three summers, the general seasonal trends for and were comparable (see Figure A2). exhibited a strong increase directly after snowmelt and peaked in late June, with maximum values of 300 W m -2 . It then decreased throughout the summer before settling on a magnitude comparable to that of at the end of the summer. Unlike the clear peak and fall seen for , did not follow a clear seasonal trend (Figure A2). It stayed rather constant throughout the season and shows only a small peak towards the middle of the summer in the last two years, reaching maximum values of nearly 200 W m -2 . Consequently, BR is greater than 1 at the beginning of summer and subsequently decreases, dropping below 1 at the end of the summer.

By comparing the three years with one another, it is evident that a delayed start to the summer season in 2018 did not have a large impact on the overall patterns of turbulent heat fluxes. in 2019 showed a similar pattern to 2017, only differing in the second half of July due to an extended warm spell (Figure A1). In summary, it appears that evolves in a consistent way and follows the same trend after snowmelt, with some differences as a result of meteorological conditions.

c. Evapotranspiration and Soil Properties

In this section, we look at one of the limiting factors for evapotranspiration: water availability.

To further investigate whether a lack of water restricts the latent heat flux, Figure 6 displays the precipitation, evapotranspiration (ET), and soil water content in summer 2017 (summer 2018 and 2019 are available in the supplementary material).

Compared to other sites at similar latitudes, Umiujaq experiences high levels of precipitation (Figure 4), meaning there is no general lack of water supply at the surface. Rainfall in Umiujaq is also fairly equally distributed over the entire summer (Figure 2), therefore no distinct dry or wet periods were observed for any of the three summers (the corresponding figures for 2018 and 2019 are available in the supplementary material). However, as the red curve in Figure 6 indicates, the top soil volumetric water content under the lichen stayed rather low even during precipitation events, with peaks below 0.15 m 3 m -3 in 2017 and 0.2 m 3 m -3 for the three summers (the corresponding figures to Figure 6 are available in the supplementary material). Following rainfall, dropped rapidly to well below 0.1 m 3 m -3 , and at times dropped to 0.05 m 3 m -3 . This demonstrates that even though there is abundant water in the form of rain, the soil can only hold small amounts of the water. This is because the water drains rapidly due to the large particle size [START_REF] Gagnon | The carbon sink due to shrub growth on Arctic tundra: a case study in a carbon-poor soil in eastern Canada[END_REF] of the soil and the associated high hydraulic conductivity.

Two post-rain periods in Figure 6, one at the end of June and the other in the beginning of July, very clearly demonstrate the influence of soil water content on ET. Approximately two days after the drop in surface water content, ET declined from about 2 mm day -1 to less than 0.5 mm day -1 and then rose again in response to the next rainfall event. The presence of a ≈2-day lag time between both variables is particularly visible until mid-July as illustrated by the 3-day running mean of ET.

The running mean of ET also highlights the influence of soil moisture on ET. Consistently, the opposite behaviour was observed for BR, which increased after the rainy periods.

d. Comparing Umiujaq to Other High Latitude Sites

Figure 7 illustrates the mean net radiation and sensible and latent heat fluxes for the comparison sites (see Table 1), exclusively during daytime ( > 0). The observed net radiation peaked at the beginning of the study period or shortly after, and subsequently decreased. Maximum net radiation values tend to drop at higher latitudes (CA-NS recorded more than 300 W m -2 at 56 • N and the maximum mean value at RU-COK was 180 W m -2 at nearly 71 • N).

The general behavior of and can be roughly divided into three groups: that is typically higher than (BR<1; RU-Cok, US-ICH, RU-VRK, US-EML), and that are approximately equal (BR≈1;US-AN and CA-NS) and finally that is lower than (BR>1; CA-UM (Umiujaq)).

We argue that this general behavior of and is usually influenced by three main factors: (i) meteorological conditions, (ii) vegetation cover, and (iii) soil type. Factors (i) and (ii) are discussed in the remainder of this section, while factor (iii) is discussed in section 3e.

(i) There were considerable differences in the mean annual precipitation for all sites analyzed, ranging from 140 mm (US-AN) to more than 800 mm in Umiujaq. However, there is no correlation between the annual precipitation and the Bowen ratio or the magnitude of the latent heat flux. This can be seen at the drier sites in Alaska, which have BRs that are similar to or lower than the much more humid RU-VRK site. Moreover, the highest BR values were observed at Umiujaq, the most humid site. Thus, the annual mean precipitation does not limit , although in the aftermath of a rain event, levels will rise as seen for Umiujaq before.

(ii) Vegetation cover also alters the BR through transpiration. Through the shade they provide, shrubs cool ground temperature by reducing incoming radiation and consequently influence the relative importance of and . The relative importance was not analyzed as it was not available for all sites. Even though all the sites were classified as shrubland, the composition of the vegetation at these sites varied. Table 1 provides a brief overview of the vegetation at each site, from the sparse tussock tundra at US-AN with only a few shrubs, to the CA-NS site's dense shrub vegetation with small trees. However, in regards to precipitation, no correlation between the vegetation type and the Bowen ratio was found.

e. Soil Types

As discussed above, the influence of precipitation and vegetation on seasonal and its behaviour relative to is limited. This leaves soil as the dominant factor that controls the BR at Umiujaq.

Figure 8 shows the approximate position of each site in a soil texture diagram. The estimate for Umiujaq is based on [START_REF] Gagnon | The carbon sink due to shrub growth on Arctic tundra: a case study in a carbon-poor soil in eastern Canada[END_REF], while the others are based on the description of the soil from the references presented in Table 1.

Although the positioning of each site within the soil texture triangle is highly uncertain, there is a clear distinction between Umiujaq and the other sites: all the other sites exhibit more clay and silt-rich soils compared to Umiujaq. Sand-rich soils are known to hold less water, drain quicker, and are thus often drier compared to soils rich in clay and silt. For instance, the mean soil moisture for the sites US-ICH (0.39 -0.52 m 3 m -3 ) and US-AN (0.56 m 3 m -3 ) are considerably higher than for Umiujaq (0.1 -0.2 m 3 m -3 ). According to [START_REF] Clapp | Empirical equations for some soil hydraulic properties[END_REF], the saturated volumetric water content of sand is ≈0.1 m 3 m -3 less than for clay or silty clay. The largest difference is in the saturated hydraulic conductivity: 1.056 cm min -1 for sand compared to 0.0062 cm min -1 for silty clay [START_REF] Clapp | Empirical equations for some soil hydraulic properties[END_REF].

Permafrost, which forms an impermeable layer, is also an important factor at each site (see Table 1). Sites that have a shallow active layer, defined as the soil layer that thaws in summer, have a topsoil layer that is often wet, thus favoring . At the majority of the comparison sites, the active layer is thin, ranging from 15 cm to 70 cm. At Umiujaq however, there is discontinuous/intermittent permafrost and an active layer that is quite deep (100 -150 cm), which favors the percolation of surface water.

f. Model Comparison

1) S M

Results from the previous sections highlight a particularity of the Umiujaq site, that despite its abundant precipitation, the presence of sandy soil produced a low latent heat flux. In light of this, the ability of both models to simulate soil water content was first examined. As mentioned in section 2e, two runs were performed for each model. A default run was performed with no changes to the parameters of the model and another with adjusted hydraulic parameters was performed in order to improve soil simulations. Figure 9 compares the observed soil water content at two different depths with the outputs from ISBA and SVS with default parameters and Figure 10 is the corresponding figure for the model run with adjusted parameters.

For the default simulation, only the soil composition was input into the models. ISBA grossly overestimated the soil water content at 6 cm by a margin of 0.2 m 3 m -3 during peak runoff and produced results that were comparable to the observations about one month later in the simulation.

Even after the soil water content dropped to a range that was comparable to the observations, the model still slightly overestimated the soil water content (around 0.02 m 3 m -3 ). The observed soil water content at 6 cm and 50 cm differed considerably, while only minor differences between these two layers were observed in the estimates produced by the two models. The surface water content at 6 cm was overestimated more by SVS than by ISBA, but the default simulation revealed more pronounced differences between values at 6 cm and 50 cm.

After adjustments were made (see section 2e), the ISBA outputs more closely matched the observations. When applying the modified soil properties, the simulated SVS results were slightly less accurate than ISBA at the surface when compared to the observations. SVS generally underestimated the soil water content by about 0.03 m 3 m -3 , while it showed high peaks during rain periods.

At 50 cm in depth, SVS produced results that were comparable to ISBA and to the observed value.

To simulate the turbulent heat fluxes, the adjusted models were used. Neither model was able to simulate as accurately as , with the most notable problem for both models being the nighttime fluxes (Figure 12). ISBA rarely simulated very small negative fluxes, while SVS did not simulate them at all, despite the observational dataset presenting a regular occurrence of negative fluxes. A negative indicates the condensation of water vapor onto the soil surface, and this process seems completely absent from the models. Apart from this, both models yielded a data point distribution that exhibited a more positive trend compared to their simulations. For ISBA, the bias was 5.3 W m -2 and for SVS, it was 11.1 W m -2 . Particularly, patterns with very low simulated fluxes and high observed fluxes were more frequent for than for . The q-q plot for ISBA showed only minor deviations from the 1:1 line for fluxes with an amplitude of less than 100 W m -2 . In the region above 100 W m -2 , the number of points dramatically decreased, with only a few points responsible for the high deviations seen in the q-q plot. For SVS, the datapoints diverge from the 1:1 line at a lower observed flux (≈70 W m -2 ) and the slope following the divergence is higher than for ISBA. Thus, a significant portion of the data points were overestimated, with deviations as large as 150 W m -2 for observed fluxes in the 100 W m -2 region.

The left side of Figure 13 shows that ISBA performed slightly better in simulating than SVS.

The correlation (correlation coefficient of 0.91 for SVS versus 0.97 for ISBA) and the normalized RMSD (0.44 for SVS versus 0.27 for ISBA) show better results for ISBA than for SVS. Also, the standard deviation for ISBA was similar to that of the observations, and ISBA slightly overestimated the variability in , while the overestimation of SVS is a bit higher (1.08 for SVS versus 1.02 for ISBA). These results are consistent with the more dispersed distribution of the data points in the scatter plot that indicate the underestimation of by SVS.

The Taylor diagram on the right of Figure 13 illustrates the relatively poor simulation performance of the two models for compared to which affected both models to a similar extent. (ISBA: correlation coefficient of 0.84, standard deviation of 1.2 and normalized RMSD of 0.64; SVS: correlation coefficient of 0.74, standard deviation of 1.12 and normalized RMSD of 1.64)

Discussion

a. Energy Balance Closure and Measurement Error

A typical issue with eddy covariance measurements is the lack of energy budget closure [START_REF] Foken | The energy balance closure problem: An overview[END_REF]. Throughout most meteorological conditions, the energy budget remained unclosed at the Umiujaq site (on average 70.5%). However, certain wind directions can lead to better closure than others (see Appendix Figure A4). The most probable reasons for the SEB imbalance at this site are: (I) the complex topography and (II) the mismatch between the footprint areas of the different flux sensors.

(I): As seen in Figure 1b, the flux tower was located to the east of a steep ridge. The cliff acted as a backward-facing step, potentially leading to flow detachment in the presence of westerly winds.

The eddy covariance approach is based on the assumption that there is a surface located upstream that is flat and homogeneous, which is not the case here. However, in recent years, we have noticed an increase in studies conducted in non-ideal terrain (e.g. 2015)). However, considering that the closure of the energy balance exceeded 55% for westerly wind direction (compare Figure A4), most of these data were retained. [START_REF] Foken | The energy balance closure problem: An overview[END_REF] argues that in heterogeneous terrains, large eddies can substantially contribute to the energy exchange when generated at the boundaries between different land uses that are excluded from the flux measurements.

(II): The footprint of the eddy covariance system naturally varies depending on wind speed and direction, covering an area of up to a few thousand square meters whereof 20 -30% is dominated by lichen, and the rest by shrubs (Figure A3). Radiation measurements cover areas on the order of a few square meters, whereas the ground heat flux measurement is a point measurement. Thus, radiation and ground heat flux measurements that do not represent the flux footprint inevitably, which leads to an unbalanced budget. The measurements of and do not match those of and taken in the footprint area, which may therefore lead to a budget imbalance. The radiation and ground heat flux measurements were derived from lichen-covered surfaces, but the footprint of the eddy covariance system covered a mixture of lichen and shrubs of varying heights, with either lichen or moss understory [START_REF] Gagnon | The carbon sink due to shrub growth on Arctic tundra: a case study in a carbon-poor soil in eastern Canada[END_REF]). This was however, difficult to quantify as knowledge on the distribution of each vegetation type would have been required for every 30-min flux footprint area. However, we observed a worse energy budget closure (50% to 65%) for periods during which winds blew in from the south, where shrubs tended to be taller and the topography was more complex.

To explore the spatial variability of the surface albedo around the flux tower, point measurements of spectral albedo were made over the two dominant vegetation types [START_REF] Belke-Brea | Impact of shrubs on winter surface albedo and snow specific surface area at a low Arctic site: In situ measurements and simulations[END_REF], lichen and shrubs. The results showed slightly higher albedo levels over lichen, on the order of a few percentage points over the entire spectral range. Additionally, the ground heat flux was calculated under different vegetation types and only minor differences were observed, which alter the overall EBC by 1 or 2%. The 5TM probe did not cover the organic layer at the surface thus an additional error of the ground heat flux is probable. This is because the hydraulic and thermal properties of the organic layer differ from those of the mineral soil. Therefore, although the measurement footprint mismatch certainly contributes to the energy budget imbalance, it is not the sole cause.

Furthermore, measurement errors possibly contributed to the degradation of the energy budget.

b. Importance of the Soil

As rather sparse vegetation prevails within the eddy covariance footprint (20 -30% lichen and 60 -70% small shrubs with low LAI) we conclude that it is not the dominant factor which controls the 6. The delay of evapotranspiration relative to soil moisture can likely be explained by variations in air temperature which also show an increased importance once the changes in soil moisture are rather small as observed in the second half of the period depicted in Figure 6. [START_REF] Lemieux | Groundwater dynamics within a watershed in the discontinuous permafrost zone near Umiujaq (Nunavik, Canada)[END_REF] investigated the groundwater dynamics of the Tasiapik valley and used empirical formulas to estimate that 40% of the total precipitation returned to the atmosphere as evapotranspiration. In our study, 22% to 35% of the precipitated water evaporated in the three summers of interest. There are likely two probable reasons for this difference. One was the use of empirical formulas in the [START_REF] Lemieux | Groundwater dynamics within a watershed in the discontinuous permafrost zone near Umiujaq (Nunavik, Canada)[END_REF] study that were not specifically developed for the climatic conditions in Umiujaq. The second reason is the fact that the watershed contained areas with clay-rich soil, which can retain more water than sandy soil [START_REF] Clapp | Empirical equations for some soil hydraulic properties[END_REF].

Thus, the evapotranspiration rates for some parts of the Tasiapik valley are likely to be higher than those found around the flux tower used in our study.

ISBA performed slightly better in simulating soil conditions than SVS, revealing weaknesses of using a simplified soil scheme. However, both models exhibited problems when simulating the vertical profile of the soil water content and showed no or only a slight gradient of the volumetric water content in the default configurations (Figure 9). Thus, the models accounted for only a fraction of the actual changes in the soil properties that occurred with depth.

c. Turbulent Fluxes and Summer Energy Balance

Overall, the sensible heat flux is the dominant consumer of the energy from radiation due to the striking dependence of the latent heat flux on the soil water content. On average, the sensible heat flux is 61.3 W m -2 during summer, showing a decline from 81.5 W m -2 at the beginning of summer to 48.3 W m -2 at the end. As previously discussed, the latent heat flux is more constant throughout summer with an average value of 39 W m -2 . The net input of energy from the radiation on average is 187.5 W m -2 .

Thus, in Umiujaq, the sensible heat flux makes up 33% of the net radiation on average while 22.4% can be attributed to latent heat flux. At the comparison sites, the proportion of sensible heat flux to radiation is 22% to 40% and the proportion of latent heat flux varies between 29% and 55%. Thus, sensible heat fluxes in Arctic environments seems somewhat less variable, while the differences in the latent heat flux are higher and depend more on local peculiarities such as water availability. This makes latent heat fluxes more difficult to simulate. The poor simulation of the soil at the Umiujaq site by the default configurations in both models, which leads to wetter conditions, impacted simulations of the turbulent fluxes as the latent heat flux was larger due to more water availability. After adjusting the soil properties, the ISBA simulations for these fluxes improved, indicating the influence of soil moisture on the turbulent fluxes in both models. While was accurately simulated along the entire range of observed values, this was not the case for .

The simulated values for were overestimated, closer to the higher-end values of the observed . This difference may be attributed to the fact that the parametrization of is less adapted to the dry soils at that site.

Both SVS and ISBA apply the same parametrization for the above-canopy turbulent fluxes.

Therefore the reason for the more accurate simulations in comparison to most probably lies in the soil simulation. Both the calculations for the thermal and the hydraulic regime are considerably simpler in SVS. As and depend on the temperature and moisture content at the surface, respectively, an error in those variables directly affects the calculation of and .

Conclusions

We analysed the SEB of a Subarctic experimental site in eastern Canada using the eddy covariance technique. We then compared the SEB to other Arctic and Subarctic sites and assessed two current LSMs, ISBA and SVS, in their ability to simulate SEB at the experimental site.

Our results showed that the sensible heat flux was the most important way of dissipating the energy from net radiation. The sensible heat flux accounted for about 33% of net radiation while the latent heat flux made up only 22.4%. The experimental site near Umiujaq consisted of sandy soil that resulted in very low soil water content at the surface, despite abundant precipitation for this latitude. The dry soil conditions most likely constrained the latent heat flux by limiting the water available for evaporation. This condition favored the dissipation of energy through sensible heat flux rather than through latent heat flux. Thus, we conclude that soil has a major influence on the partitioning of the turbulent heat fluxes in contrast to other studies that favor vegetation as the dominant factor for this partitioning.

Other documented sites in the Arctic showed net radiation partitioning that was different. At these sites, the latent heat flux was larger than the sensible heat flux. We attribute this difference to soil properties that retain more water, allowing it to drain slower.

We have also shown that both LSMs either need to be tuned or are not sophisticated enough to account for this particular situation concerning soil content. For this reason, we adjusted the models to be able to account for this factor. However, once the soil water content was better simulated through the modification of the models, the sensible heat flux in particular was very accurately simulated. On the other hand, latent heat flux simulations were not as accurate. Nonetheless, the results confirm that both LSMs are generally suitable for operational hydrometeorological applications in the humid Subarctic tundra ecotone.

Due to the isostatic uplift that followed deglaciation and led to the emergence of beaches, we stress that the sandy soil that covers a significant portion of Arctic T 2: Hydrometeorological conditions for the three summers of interest. Soil temperature measurement depth is 5 cm. 1.

LIST OF FIGURES

  3 • C in eastern Russia (RU-COK), by far the coldest site, to -1 • C in Alaska (US-EML). Mean annual precipitation varies considerably, from the dry Alaskan site (122 mm; US-AN) to more humid sites in western Russia (538 mm; RU-VRK) and eastern Canada, with Umiujaq being by far the most humid site (800 -1000 mm).

Figure 5

 5 Figure 5 illustrates the cumulative values of net radiation , sensible heat flux , latent heat

  Figures 11 and 12 show observations compared with simulated sensible and latent heat fluxes

  [START_REF] Mauder | A strategy for quality and uncertainty assessment of long-term eddy-covariance measurements[END_REF] reported errors of 10%-15% for the eddy covariance data for flags 0 and 1, which were used in this study. Radiation and ground heat flux measurements yielded errors of up to 10% (compare https://www.campbellsci.com/cnr4). Although the energy balance closure was rather poor, taking climatic conditions into account, it was still in the range of other carefully designed experiments (e.g. Cui and Chui (2019)). The energy budget closure could most probably be improved by additional radiation measurements and by relocating the tower further away from the steep ridge to reduce perturbations of the mean flow for westerly winds.

Fig. 1 .Fig. 2 .Fig. 3 .Fig. 4 .Fig. 5 .Fig. 6 .Fig. 7 .Fig. 8 .Fig. 9 .Fig. 10 .Fig. 11 .Fig. 12 .Fig. 13 .Fig. A3 .Fig. A4 .

 12345678910111213A3A4 Fig. 1. a) Location of the study site and the experimental setup in the Tasiapik valley close to the village of Umiujaq on the eastern shore of Hudson Bay, Canada and b) photo of the upper portion of the valley where the flux tower is installed. The angle of view in panel b) is shown in panel a). Note the presence of a cuesta west of the tower. . . . . . . . . . . 46 Fig. 2. Mean, high and low 2-m air temperatures as well as the mean amount of snow and rain per month. The shaded area indicates the period of interest for this study. . . . . . . . . 47 Fig. 3. Soil profile under lichen cover in the vicinity of the flux tower. . . . . . . . . . 48 Fig. 4. FLUXNET site locations and spatial pattern of circumpolar yearly precipitation. Source: Harris et al. (2014) . . . . . . . . . . . . . . . . . . . . . 49 Fig. 5. Partitioning of net radiation into , , and for the three summers. The residual after the partitioning is also given. Only data with EBC between 0.2 and 1.5 were used. . . . 50 Fig. 6. Daily sum of evapotranspiration (ET) and precipitation and daily mean of soil water content (SWC) at 6-cm depth, for summer 2017 (June to August). . . . . . . . . . . . 51 Fig. 7. Mean net radiation (black lines), sensible and latent heat fluxes (red and blue lines; uncorrected for energy imbalance) of all available summers and selected study sites. The shaded areas illustrates the mean standard deviation. Only data during daytime ( > 0) were used. Sites are presented in ascending order of latitude. The bottom right plot shows the mean of Umiujaq compared to the mean for all sites. For RU-COK only one year of data was available. . 52 Fig. 8. Ternary diagram of soil textures according to the USDA classification. Circles indicate the estimates for each site. For site US-An, no data about the soil was available. The color of the circles informs on the BR: purple: BR<1; orange: BR≈1 and green: BR>1. Source: see references in Table 1. . . . . . . . . . . . . . . . . . . . . 53 Fig. 9. Observed and modeled soil water content at 2 different depths for 2017 for models with default parameters. . . . . . . . . . . . . . . . . . . . . . 54 Fig. 10. Observed and modeled soil water content at 2 different depths for 2017 for models with adjusted parameters. . . . . . . . . . . . . . . . . . . . . . 55 Fig. 11. Left column: Scatter plots comparing observed and modeled sensible heat fluxes at the Umiujaq site for all summers with the colors indicating the hour of day for (a) ISBA and (c) SVS. Right column: q-q plots, for which the quantiles of the distributions are plotted against each other, for (b) ISBA and (d) SVS. . . . . . . . . . . . . . . . . 56 Fig. 12. Same as Figure 11 but for latent heat flux . . . . . . . . . . . . . . . 57 Fig. 13. Taylor diagrams indicating how closely modeled sensible heat flux (left) and latent heat flux (right) match observations. . . . . . . . . . . . . . . . . . . . 58 Fig. A1. Weekly running mean of the air temperature for summer 2017, 2018 and 2019. . . . . . 59 Fig. A2. Weekly running mean of summer sensible and latent heat fluxes for all three study years.Summer 2017 starts later as data were available only from mid June.. . . . . . . . 60

  Ternary diagram of soil textures according to the USDA classification. Circles indicate the estimates for each site. For site US-An, no data about the soil was available. The color of the circles informs on the BR: purple: BR<1; orange: BR≈1 and green: BR>1. Source: see references in Table

Fig. A1 :

 A1 Fig. A1: Weekly running mean of the air temperature for summer 2017, 2018 and 2019.

Fig. A2 :

 A2 Fig. A2: Weekly running mean of summer sensible and latent heat fluxes for all three study years. Summer 2017 starts later as data were available only from mid June.

Fig. A3 :

 A3 Fig. A3: Drone image of the study site with the flux footprint climatology of summer 2017, indicating 80% of the source area.

Fig. A4 :

 A4 Fig. A4: Dependence of the energy budget closure on the wind direction for all three summers.

  .1% of the hourly values of the sensible heat flux (2017: 50.8%; 2018: 63.5%:

	2019: 54.2%) and 47.8% of the latent heat flux (2017: 46.7%; 2018: 54.2%: 2019: 42.4%) passed
	all quality control steps.

  Geissbühler et al. (2000); Hiller et al. (2008); Etzold et al. (2010); Nadeau et al. (2013); Stiperski and Rotach (

  Canada causes energy partitioning that is very different from what is presented in other available studies. A priori soil conditions found at the Tasiapik valley site seem uncommon. However, consulting this map (https://ftp.maps.canada.ca/pub/nrcan_rncan/publications/STPublications_ PublicationsST/295/295462/gsccgm_195_b_2014_mn01p1.pdf), it can be seen that many regions in the Arctic have similar soil conditions. Indeed, large parts of the coasts of Hudson Bay and Ungava Bay, as well as other parts of the coast of numerous other regions in the eastern and central Canadian Arctic are covered by coarse sediments such as sand and gravel. Moreover, important parts of the Canadian Arctic comprise rocky outcrops that likely show even higher Bowen ratios. The soil water content and latent heat flux might therefore be overestimated by surface models across the Canadian Arctic. Furthermore, the disappearance of permafrost due to climate change may lead to drier soil and subsequently, altered latent heat fluxes. Description of the study sites. MAT and MAP are the mean annual temperature and the mean annual precipitation. C: Continuous, DC: Discontinuous, S: Sporadic.

	T	Site Location Coordinates Elevation (m ASL) MAT ( • C) MAP (mm) Vegetation Soil Type Permafrost Reference 2017 1 July -5 September CA-UM Quebec, Canada 56.33 • N 76.28 • E 140 -4.0 800-1000 shrub tundra sand DC/S This study 2018 2019 1: T CA-NS US-ICH Manitoba, Canada Alaska, US 55.92 • N 98.96 • E 68.61 • N 149.30 • E 244 940 -3.1 -7.4 495 318 heath tundra tussok tundra clay porous peat DC/S C Goulden (2016) Bret-Harte et al. (2016) ( • C) 9.6 11.5 8.1 9.3 10.1 11.7 Mean 9.3 10.8	US-EML Alaska, US 63.88 • N 149.25 • E 700 -1.0 378 tussok tundra till/ loess C Belshe et al. (2012) Precipitation US-AN Alaska, US 68.93 • N 150.27 • E 600 -11.0 122 shrub tundra n/a C Rocha et al. (2016) (mm) Maximum snow height RU-VRK RU-COK Vorkuta, Russia Chokurdakh, Russia 67.05 • N 62.94 • W 70.83 • N 147.49 • W 100 48 -6.1 -14.3 538 212 shrub tundra shrub tundra loamy and peaty fluvial silt DC/S C Biasi and Friborg (2016) van der Molen et al. (2007) (cm) in previous winter 237 65 290 140 369 100 298 101

Data availability statement. Data available on request from the authors.

Acknowledgments. This study was funded by Sentinel North, a Canada First Research Excellence Fund, under the theme 1 project titled, "Complex systems: structure, function and interrelationships in the North". We are grateful to Denis Sarrazin for technical support of the CEN weather stations and the community of Umiujaq for permission to conduct our research. We acknowledge helpful and very constructive reviews by Werner Eugster and 2 anonymous reviewers.

APPENDIX

Supplementary Figures

LIST OF TABLES

, and for the three summers. The residual after the partitioning is also given. Only data with EBC between 0.2 and 1.5 were used.