Vegetation type is an important predictor of the arctic summer land surface energy budget

Surface energy fluxes and components

In our study, we focused on the circumpolar land north of 60° latitude, and specifically on the extent of the circumpolar Arctic vegetation map (CAVM²⁰, Supplementary Fig. 1–3). We obtained half-hourly and hourly in situ observations of energy fluxes and meteorological variables from the monitoring networks FLUXNET²⁸ (fluxnet.org; FLUXNET2015 dataset), AmeriFlux²⁹ (ameriflux.lbl.gov), AON^31,32 (aon.iab.uaf.edu), ICOS (icos-cp.eu), GEM^35,36 (g-e-m.dk), GC-Net^33,34 (cires1.colorado.edu/steffen/gcnet) and PROMICE³⁰; (promice.dk; Supplementary Table 3). We did not include observations from the Baseline Surface Radiation Network (BSRN; bsrn.awi.de) and Global Energy Balance Archive (GEBA; geba.ethz.ch) because they typically lack information on non-radiative energy fluxes. Finally, we did not include observations from the European Flux Database Cluster (EFDC, europe-fluxdata.eu) because these data are largely located outside the domain of the CAVM²⁰.

We aggregated surface energy fluxes and components (Supplementary Table 1) to daily resolution as follows: (i) we extracted only directly measured data and excluded gap-filled data by filtering according to quality information; (ii) we performed a basic outlier filtering (excluding shortwave and longwave radiation flux values >1400 Wm⁻² and in case of incoming/outgoing radiation <0 Wm⁻², excluding albedo values <0 and >1, excluding air and surface temperatures < −100 °C; (iii) we converted all times of measurements to local standard time (i.e. without daylight saving time); (iv) we calculated daily (i.e. 24 h) mean, minimum and maximum values for all days and surface energy fluxes and components where a minimum of 65% of data was available with a maximum temporal gap of 4.8 h; (v) we extracted metadata (units, measurement heights, original variable names, instrumentation) for all sites and variables; (vi) when data for the same time and location were available from several networks (partial overlapping of FLUXNET and AON with Ameriflux), we averaged respective daily values across networks; (vii) in case data for the same time and location were measured by several sensors at one site (replicated measurements), we averaged across replicates; (viii) we harmonized units for all variables; (ix) we harmonized the flux direction convention for non-radiative energy fluxes (i.e. H, LE, and G; Supplementary Table 1) as „positive away from surface“; (x) we derived net radiation (R_net), net shortwave radiation (SW_net), net longwave radiation (LW_net) and albedo from corresponding daily aggregated incoming and outgoing shortwave and longwave radiation, respectively, if not otherwise available; (xi) we derived normalized fluxes for R_net (n.R_net), SW_net (n.SW_net), LW_net (n.LW_net), H (n.H), LE (n.LE) and G (n.G) as the percentage (%) of daily maximum potential incoming shortwave radiation based on location and topographical conditions⁵⁹. The final surface energy budget dataset (named SEB-data hereafter) consists of daily mean, minimum and maximum values for surface energy fluxes and components for 64 tundra and glacier sites across the years 1994–2021 (Supplementary Figs. 1–3; Supplementary Table 2; nr. sites = 64, nr. site years = 652).

Drivers of the surface energy budget (SEB-drivers)

We characterized all study sites available in the SEB-data according to environmental covariates (named SEB-drivers hereafter; Supplementary Table 1) that potentially have important effects on the magnitude and seasonality of surface energy fluxes from ancillary geographic data (Supplementary Fig. 4, Supplementary Table 3). We identified these relevant SEB-drivers according to our literature synthesis and consequent correlation analysis (Supplementary Discussion, Supplementary Figs. 7–9).

As the main SEB-driver of interest, we focused on the local-scale, in situ vegetation type (Vegetation type variable), which reflects the main classes described in the CAVM²⁰ map (barren complexes, graminoid tundra, prostrate dwarf-shrub tundra, erect-shrub tundra, wetland complexes, glacier), plus boreal peat bog.

To derive this in situ vegetation type variable, we extracted for each study site in the SEB-data the vegetation descriptions from adequate literature references (Supplementary Table 2). We categorized each site’s vegetation according to the most adequate CAVM-class by using the following decision chain: (i) assign the CAVM-class where most described species from in situ vegetation descriptions and Table S1 in ref. 20 are in accordance. Take into account species‘ ecological niche sizes (e.g. stenotopic plants are better indicators than eurytopic plants); (ii) if there is no species list and/or in case there are several CAVM candidate classes, take into account the general habitat and ecosystem descriptions (including community vertical and horizontal structure, soil moisture and acidity, topography); (iii) if there are still several CAVM candidate classes, take into account the bioclimatic subzone and evaluate the exact location of the site with available satellite data; (iv) if there are still several CAVM candidate classes, assign the most dominant CAVM-class; (v) in case of uncertainties, consult with the study author(s); (vi) If the vegetation is not similar to any CAVM-class, describe it according to the description of the study author(s) and/or adequate references (Supplementary Table 2).

This categorization revealed 1 barren tundra site (1 B3), 12 graminoid tundra sites (10 G4, 1 G3 and 1 G1), 3 prostrate-shrub tundra sites (2 P2 and 1 P1), 4 erect-shrub tundra sites (2 S2 and 2 S1), 8 wetland sites (3 W3, 3 W2 and 2 W1), 33 glacier sites (GL), and three sites that were not similar to any CAVM-class but clearly identified as boreal peat bogs (Supplementary Figs. 1–3; Supplementary Table 2).

For each study site in the SEB-data, we extracted the landscape-scale, dominant vegetation type in the site-surrounding area (radius = 500 m, CAVM type variable) from the raster version of the CAVM^20,60 (Supplementary Fig. 4a; Supplementary Table 3). This map is a refined version of the widely used original version^20,60 and provides information at 1 km² spatial resolution and pan-Arctic scale north of the Arctic treeline²⁰. Compared to the Vegetation type variable, the CAVM type is based on the analysis of a combination of large-scale satellite and environmental data, but both variables refer to the same class definition. The CAVM-classes are based on the plant physiognomy of the zonal vegetation in a given area, analogous to the widely applied Braun-Blanquet approach for plant communities on the ground²⁰.

The CAVM types used in this study reflect these main plant physiognomies (plus one glacier class) that are named according to the following dominant plant growth forms: B: barren and barren complexes (4 sites), G: graminoid tundra (10 sites), P: prostrate dwarf-shrub tundra (1 site), S: erect-shrub tundra (10 sites), W: wetland complexes (3 sites), GL: glacier (30 sites; there are 6 non-Arctic sites in our SEB-data). A confusion matrix showing the classification of study sites into Vegetation type and CAVM type categories is shown in Supplementary Table 6.

For each study site, we averaged the mean annual air temperature (°C), the annual sum of precipitation (mm) and the annual sum of snow water equivalent (Snow amount variable, mm; defined as daily precipitation when daily mean air temperature < =0 °C) across the years 1979-2018. Therefore, we used high resolution downscaled model output estimates of temperature and precipitation (CHELSA V2.1; „tas“ and „pr“ variables^61,62,63; Supplementary Fig. 4d, e; Supplementary Table 3).

Using the same air temperature data as above^61,63, for each study site we calculated the average Conrad Continentality Index (CCI⁶⁴) across the years 1979–2018. Specifically, we used the formula:

whereby T_max and T_min (°C) refer to the mean air temperature of the warmest and coldest month, respectively, and φ (radians) refers to latitude.

Using the same air temperature data as above^61,63, for each study site we calculated the average Summer Warmth Index (i.e. the annual sum of monthly mean air temperatures above 0 °C; SWI⁴²) across the years 1979–2018.

For each study site, we extracted the bioclimatic subzone (CAVM subzone variable) as described in the circumpolar Arctic vegetation map (CAVM^20,60). These five bioclimatic subzones (A-E) and additional classes for glacier and non-Arctic zones are largely based on a combination of Arctic phytogeographic zones⁶⁵, dominant growth forms of plants and summer temperatures⁶⁶. Hence, these bioclimatic subzones are generally well aligned with summer warmth index (SWI) classes⁴². A confusion matrix showing the classification of study sites into Vegetation type and CAVM subzone categories is shown in Supplementary Table 7.

For each study site, we extracted the median snow cover duration (Snow duration variable) from satellite-sensed daily snow cover (MODIS MOD10C1⁴¹; Supplementary Table 3). Snow duration was calculated for each year as the number of days between the spring snow-free date (i.e. mean of days of year for last ‚snow‘ day and first “no snow” day MODIS categories) and the autumn snow-onset date (i.e. mean of days of year for last “no snow” day and first “snow” day MODIS categories). Yearly snow duration, snow-free date and snow-onset date for each study site were then aggregated by calculating the median across the years 2000–2020.

For each study site, we extracted mean annual cloud cover (%) and cloud-top temperature (°C) from monthly satellite imaging radiometer data (“cldamt” and “tc” products; ISCCP-Basic-H series^67,68; Supplementary Fig. 4f, g; Supplementary Table 3) and averaged these variables for each study site across the years 1984–2016.

For each study site we extracted the corresponding permafrost extent as described in the Circum-Arctic permafrost and ground ice map⁶⁹ (NSIDC gdd318_map_circumArctic version 2; Supplementary Fig. 4b; Supplementary Table 3). This map describes 4 categories of permafrost extent, based on the percentage of the ground that is underlain by permafrost. There are additionally 5 separate categories for glaciers, relict permafrost, inland lakes, oceans and land with no permafrost. For the Permafrost extent variable used in this study, we aggregated these categories into the following five classes: continuous permafrost (C; 90–100% extent; 21 sites); discontinuous permafrost: (D; 50–90% extent; 5 sites); sporadic or isolated patches of permafrost (Si; <50% extent; 4 sites), ocean/inland seas (o; 2 sites) and glaciers (g; 32 sites). A confusion matrix showing the classification of study sites into Vegetation type and Permafrost extent categories is shown in Supplementary Table 8.

Using the same permafrost data⁶⁹ (Supplementary Fig. 4b; Supplementary Table 3), for each study site we extracted the corresponding permafrost ground-ice content class: high (h; >20% ice content; 8 sites), medium (m; 10–20% ice content; 7 sites), low (l; 0–10% ice content; 15 sites), ocean/inland seas (o; 2 sites), and glaciers (g; 32 sites). A confusion matrix showing the classification of study sites into Vegetation type and Permafrost ice content categories is shown in Supplementary Table 9.

For each study site we extracted the average altitude (m a.s.l.), slope (°) and northness of the aspect (1 if north-exposed, −1 if south-exposed; derived from the cosine of aspect in radians) in the surrounding area (radius = 500 m) from the satellite-sensed digital elevation model raster mosaic at 100 m spatial resolution (ArcticDEM: Arcticdem_mosaic_100 m_v3.0; ref. 70; Supplementary Fig. 4c; Supplementary Table 3).

Data analysis

For our analyses, we focused on the surface energy fluxes and components net radiation (R_net, Wm⁻²), sensible heat flux (H, Wm⁻²), latent heat flux (LE,Wm⁻²), ground heat flux (G, Wm⁻²), net shortwave radiation (SW_net, Wm⁻²), net longwave radiation (LW_net, Wm⁻²), albedo (unitless), and surface temperature (T_surf, °C), air temperature (T_air, °C) and the difference between surface and air temperature (T_surf-T_air; Supplementary Table 1). We included albedo and T_surf because they are directly related to the radiative SEB as follows¹⁴:

where SW_in and LW_in are the incoming shortwave and longwave radiation, respectively (Wm⁻²), ε is the surface emissivity (≅1) and 훔 is the Stefan-Boltzmann constant¹⁴. We repeated our analyses for normalized fluxes of R_net (n.R_net), SW_net (n.SW_net), LW_net (n.LW_net), H (n.H), LE (n.LE) and G (n.G), which are all expressed in percent of maximum potential incoming shortwave radiation⁵⁹ (%; Supplementary Table 1).

We conducted all data processing and analyses using the R-software version 3.6.0⁷¹. We analyzed the full SEB-dataset (nr. sites = 64, nr. site years = 652) and separately a “vegetation” data subset excluding glacier sites (nr. of sites: 31, nr. of site years: 234). Furthermore, we conducted all analyses at the yearly timescale and for the summer season (JJA; June, July, August), which is when in situ observations for vegetation sites have more complete data coverage (Supplementary Discussion). Summer is also the time when the SEB is dominated by absorbed solar radiation, and local controls by the land surface are expected to be more important than in the winter season^27,38. In the winter season, absorbed solar radiation is negligible and the SEB is largely influenced by synoptic processes, such as advection from lower latitudes^27,38. We chose to use the standard meteorological summer season (as opposed to “snow-covered” vs. “snow-free” season) because (1) the JJA season is largely snow free for relevant cases (i.e. except glacier sites), and (2) standard seasons are consistently defined in time; reducing confounding of results due to seasonal changes in solar irradiance.

We compared the relative importance of SEB-drivers for explaining variance in the surface energy fluxes and components R_net, H, LE, G, SW_net, LW_net, albedo and T_surf. We repeated this analysis for normalized fluxes n.R_net, n.SW_net, n.LW_net, n.H, n.LE and n.G (Fig. 1 and Supplementary Fig. 5). Specifically, we focused on the importance of the SEB-driver “Vegetation type” compared with the importance of other drivers related to landscape-scale dominant vegetation type (CAVM type), climate (Temperature, Precipitation, Snow amount, Continentality, Summer warmth, CAVM subzone), clouds (Cloud cover, Cloud temperature), snow (Snow duration), permafrost (Permafrost extent, Permafrost ice content), topography (Altitude) and geographic location (Latitude). These 15 SEB-drivers were selected based on our previous literature synthesis and consequent correlation analysis (Supplementary Discussion, Supplementary Figs. 7–9).

Using the “vegetation” data subset excluding glacier sites (nr. of sites: 31, nr. of site years: 234), we averaged surface energy fluxes for each summer season (JJA), year and study site, where at least 80% of daily measurements were available. To calculate the relative importance for each SEB-driver, we applied a variance partitioning method^72,73 to predict surface energy flux values averaged for each study site. Specifically, we used the set of 15 selected SEB-drivers to build all possible models with 2 predictors (a number high enough to allow for the pairwise assessment of statistical confounding among predictors and low enough to avoid model overfitting). For each model, we quantified the variance explained by each predictor in a predictor pair when fitted first, when fitted last and when averaged over all possible orderings in the models. Finally, we averaged the “first”, “last” and “average” explained variance (%) for each SEB-driver across all models for each surface energy flux variable. This led to the testing of 105 unique SEB-driver pairs × 2 predictor orders × 14 surface energy flux variables = 2940 models for the vegetation data subset.

We estimated the magnitude of the surface energy fluxes R_net, SW_net, LW_net, H, LE and G for each Vegetation type at seasonal (JJA) and yearly timescale (Y) using the full SEB-dataset (nr. sites = 64, nr. site years = 652; Table 1; Supplementary Table 4). We repeated this analysis for the monthly aggregated, normalized (i.e. potential incoming radiation-adjusted⁵⁹) summer fluxes n.R_net, n.SW_net, n.LW_net, n.H, n.LE, n.G, albedo and T_surf using the vegetation data subset (excluding glacier sites; nr. of sites: 31, nr. of site years: 234; Supplementary Table 5).

In both analyses, we averaged the daily mean values of the surface energy fluxes for each timescale, year and study site, where at least 80% of days with measurements were available. We then averaged surface energy flux values across years for each study site and estimated the mean ± 95% confidence interval (CI) as a function of Vegetation type by using a linear mixed-model analysis. Specifically, we modeled the study-site aggregated means of each surface energy flux as a function of Vegetation type (fixed effect) and the corresponding data distribution network (i.e. Ameriflux, FLUXNET etc.; random effect; Supplementary Table 2). To compare differences of summer surface energy flux estimates among Vegetation types, we applied a consequent post-hoc pairwise comparison with bonferroni correction of significance estimates⁷⁴ (Supplementary Table 4).

We derived the typical seasonal change of the surface energy fluxes R_net, SW_net, LW_net, H, LE and G for each Vegetation type (Fig. 2). We repeated this analysis for the normalized (i.e. maximum potential incoming radiation-adjusted⁵⁹) fluxes n.R_net, n.SW_net, n.LW_net, n.H, n.LE, n.G, albedo, T_surf, T_air and T_surf-T_air (Supplementary Fig. 6).

Specifically, using data constrained to the years 2000–2021 and excluding barren vegetation type (B; because of missing Rnet data; nr. of sites = 61, nr. site years = 617), we averaged the daily mean values of the surface energy fluxes across all available years for each day of year (DOY) and study site. In a second step, we grouped study sites by Vegetation type and derived mean ± s.e. for each DOY and surface energy flux variable. Finally, we smoothed the resulting mean ± s.e. values for each Vegetation type by calculating the centered 15 day moving average for each DOY (Fig. 2).

We used these smoothed seasonalities for the five Vegetation types Boreal peat bog, Wetland complex, Graminoid tundra, Erect-shrub tundra, Prostrate-shrub tundra (excluding Glacier; nr. of sites = 28, nr. site years = 217), to assess the start and the end of the “summer-regime” period for R_net, H, G, albedo and T_surf. In the case of R_net, H, G and T_surf, we defined the summer-regime as the time period where values exceed 0 Wm⁻² and 0 °C, respectively. In the case of albedo, we defined the summer-regime as the period when values fall below the mean of the yearly minimum and maximum value^75,76. We excluded LE since fluxes are >0 Wm⁻² all year in most cases. Using this data, we compared the timing of the summer regime of R_net, H, G, T_surf and albedo to the timing of the snow-free period (Snow duration⁴¹) for each Vegetation type. Specifically, we aggregated the spring snow-free and autumn snow-onset dates for each Vegetation type and subtracted these from the start and end of the summer-regime period of the selected surface energy fluxes, respectively. Finally, we aggregated the start and end of summer-regime timings across Vegetation types for each selected surface energy flux variable, to test the differences in the timing of the surface energy flux summer-regime and snow-free/snow-onset dates, using corresponding two-sample Welch’s t-tests (Fig. 3).

Source: Ecology - nature.com