2.1 Constructing a conceptual surface energy budget

The SEB is defined as the equilibrium between net radiation ${\mathit{Q}}_{\mathrm{S}}^{*}$ at the surface and the sum of sensible heat flux Q_H, latent heat flux Q_LE, and ground heat flux Q_G:

The change of heat storage at the surface was neglected, as no vegetation or built infrastructure is present in this environment. In the summer season, we expect evaporation to be the dominating phase change of ice in valley floors generating atmospheric latent heat flux, which is why we defined Q_LE as the evaporative heat flux. The magnitude of the ground heat flux in a cold-desert energy budget can be substantial because thawing of ice beneath the ice table – the upper boundary of currently frozen soil – requires a large amount of energy (Lloyd et al., 2001). We divided the ground heat flux Q_G into the temporal change in heat storage in the thawed layer ΔS_TL and into the soil heat flux Q_IT at ice table depth, which consists of the latent and sensible heat fluxes into the frozen layer beneath the ice table (Fig. 1):

where ΔS_TL equals the temporal heat storage change between surface and ice table depth z_IT(t) and was determined using a calorimetric approach following Liebethal and Foken (2007):

where C_G(z) ($\mathrm{J}{\mathrm{m}}^{-\mathrm{3}}{\mathrm{K}}^{-\mathrm{1}}$) is the volumetric heat capacity of the soil and z denominates depth in the soil.

Figure 1Surface energy balance concept used in this study, with net radiation ${\mathit{Q}}_{\mathrm{S}}^{*}$, sensible heat flux Q_H, latent evaporative heat flux Q_LE, temporal storage change in the thawed layer ΔS_TL, and soil heat flux at ice table depth Q_IT. Ground heat flux is given by adding together ΔS_TL and Q_IT. Black arrows depict observed directions of energy fluxes. The soil profile can be divided by the active layer depth z_AL into the seasonally thawed active layer and the perennially frozen permafrost, or by the ice table depth z_IT into the currently thawed and frozen layers, where the ice table equals the thawing front advancing into the soil, denoted by grey arrows.

Download

The turbulent sensible heat flux and latent evaporative heat flux were computed using the eddy-covariance method in combination with flux footprint modeling. The micrometeorological eddy-covariance method is a tool for quantifying turbulent energy and mass net exchange between the land surface and the atmosphere with high precision, and is used at more than 700 sites globally across all climate zones and biomes (e.g., Baldocchi, 2003; Aubinet et al., 1999). It is based upon turbulence observations and the application of the Reynolds decomposition ($X=\stackrel{\mathrm{‾}}{X}+X^{\prime }$), which divides a scalar or vector quantity X into its temporal mean $\stackrel{\mathrm{‾}}{X}$ and temporal perturbation X^′ to compute the net flux. According to Reynolds' second postulate, the total vertical net flux $\stackrel{\mathrm{‾}}{\mathit{w}X}$, which is computed from the measurements of the vertical wind speed w and a quantity X, can be calculated as

where the covariance $\stackrel{\mathrm{‾}}{{\mathit{w}}^{\prime }{X}^{\prime }}$ is the turbulent flux. Averaged over a sufficiently long period, the mean vertical wind speed $\stackrel{\mathrm{‾}}{\mathit{w}}=\mathrm{0}$ m s⁻¹ in the surface layer of the atmosphere. Hence, Eq. (4) is simplified to

Thus, mean vertical fluxes of scalars in the surface layer can be expressed as turbulent fluxes via covariances of w and a scalar specifying the flux. For the computation of Q_H and Q_LE, the covariance $\stackrel{\mathrm{‾}}{{\mathit{w}}^{\prime }{T}^{\prime }}$ (m K s⁻¹) of w and air temperature T and the covariance $\stackrel{\mathrm{‾}}{{\mathit{w}}^{\prime }{q}^{\prime }}$ ($\mathrm{kg}{\mathrm{kg}}^{-\mathrm{1}}\mathrm{m}{\mathrm{s}}^{-\mathrm{1}}$) of w and specific humidity q were used, respectively, to calculate heat fluxes in energetic units (W m⁻²) with Eqs. (6) and (7):

where c_p ($\mathrm{J}{\mathrm{K}}^{-\mathrm{1}}{\mathrm{kg}}^{-\mathrm{1}}$) is the specific heat capacity of air, ρ (kg m⁻³) is the air density, and λ (J kg⁻¹) is the latent heat of vaporization. Note that we used the latent heat of vaporization and not that of sublimation assuming that the water evaporating from the wetted water track soil is liquid already. The complete SEB equation applied in this study then equates to

Flux footprint modeling was used to distinguish between turbulent fluxes originating from water tracks and those from non-water track surfaces. Our intention was to isolate the effect of these small-scale, linear features, which may occupy only several tens to hundreds of square meters, onto the SEB. Flux footprint modeling allows for connecting the recorded turbulence signals at sensor location to their source area (see, e.g., Leclerc and Foken, 2014, and references therein). Based upon well-known laws for isotropic and homogeneous turbulent airflows, flux footprint models compute a spatially explicit probability density function, i.e., the flux footprint, quantifying the contribution of each grid cell of the land cover matrix to the total observed flux. Since these contributions and thus the source area vary with sensor height, surface properties, and airflow properties including turbulence statistics, wind speed, and wind direction, one may select intervals for which the turbulent flux predominantly originates from a certain land cover type of interest.

2.2 Field observations

The study was conducted in Taylor Valley, a polar desert with 3 to 50 mm annual precipitation (Fountain et al., 2009) and −18 ^∘C mean annual temperature (Doran, 2002). Continuous permafrost soils in lower Taylor Valley show active layer depths between 45 and 75 cm (Bockheim et al., 2007) and are inhabited by a simple, nematode-dominated soil ecosystem (Priscu, 1998).

Measurements were taken over 26 d under summer conditions from 26 December 2012 to 21 January 2013 at three closely collocated sites near the Ross Sea shore in the valley floor of lower Taylor Valley (Fig. 2a, Table 1), with slopes <7^∘ (Fig. S1 in the Supplement). At any time during the experiment, two eddy-covariance and SEB stations were operated: one station was installed throughout the whole period at the investigated water track (Fig. 2b), referred to as WT. The other station was operated as a reference representing the dominant non-water track bare soil surfaces in lower Taylor Valley; it was successively installed at two sites with different soil textures: PLD was located on a paleolake delta dominated by fine surficial sediments (Fig. 2c), while GT represented coarse glacial till (Table 1, Fig. 3). PLD and GT combined together are referred to as non-water track reference NWT.

Table 1Installation notes for the three eddy-covariance stations, with height of anemometer z_an, height of infrared gas analyzer z_irga, horizontal distance between both devices d_an,irga, and estimated roughness length z₀; all lengths in meters. PLD and GT are non-water track reference surfaces, and were summarized as NWT.

Download Print Version | Download XLSX

Figure 2Location of the three study sites near Ross Ice Shelf in lower Taylor Valley, with inset map showing location of the sites on the Ross Sea, gridded by 10^∘ latitude–longitude lines (a). WT site (water track) (b). PLD site (paleolake delta) with eddy-covariance and surface energy balance station (c).

Each station consisted of a net radiometer (NR01, Hukseflux Thermal Sensors B.V., Delft, NL) for measuring incoming and outgoing longwave and shortwave radiation components, and an ultrasonic anemometer (81000 VRE, R.M. Young Company, Traverse City, MI, USA) in combination with an infrared gas analyzer (LI-7500, Li-cor Inc., Lincoln, NE, USA) for eddy-covariance measurements (Fig. 2c). Sonic anemometer measurements providing wind and acoustic temperature data and infrared gas analyzer measurements of water vapor were sampled and recorded at 20 Hz (Table 1).

Soil temperatures were recorded in several depths in the thawed layer with thermistors and thermocouples (Table C1). Ice table depths were determined at WT and PLD by measuring depth to refusal (DTR). We measured volumetric heat capacity C_G of the soil in situ with a thermal properties analyzer (KD2 Pro, Decagon Devices, Pullman, WA, USA), averaged from 4 and 31 soil samples from the surface at WT and PLD, respectively, along with thermal conductivity K_G and thermal diffusivity D_G.

2.3 Calculating surface energy balance components

Eddy-covariance fluxes were computed using a fixed perturbation timescale of 30 min using the bmmflux tool software of the Micrometeorology Group of the University of Bayreuth (see Appendix in Thomas et al., 2009). First the raw data were filtered by instrument flags and plausibility limits. A despiking routine was applied to exclude unphysical turbulence data (Vickers and Mahrt, 1997). Time lags between gas analyzer and anemometer were corrected by maximizing the covariances of the measured quantities. A three-dimensional rotation routine was used to rotate the flow into the mean streamlines and eliminate the mean vertical wind potentially caused by a tilt in the sonic anemometer, surface conditions, or semi-stationary eddies of timescales exceeding the perturbation timescale (Wilczak et al., 2001). Computed fluxes were corrected for low- and high-pass filtering following Moore (1986). The buoyancy flux was converted into sensible heat flux by a post hoc buoyancy correction (Liu et al., 2001). A post hoc density correction was applied to the latent evaporative heat flux (Webb et al., 1980). Eddy-covariance quality flags for turbulent fluxes following Foken et al. (2004) were used to filter out intervals in which the assumption of stationarity and well-developed turbulence were not satisfied. The scheme runs from 1 (best quality) to 9 (worst quality) and we discarded data with flags ≥7.

ΔS_TL was determined as temporal heat storage change via calorimetry with Eq. (3), from soil temperature profiles which were logarithmically interpolated between the measurements with increments of 0.01 m. The thickness of the thawed layer was assumed to be constant in SEB calculations throughout the measuring period, although in reality it likely changed over time due to lowering of the ice table as a result of continued soil thawing in the active layer, which is typically not fully thawed before mid-January to early February in Taylor Valley (Adlam et al., 2010; Conovitz et al., 2006).

For WT, temperature at ice table depth was assumed to be 0 ^∘C. At PLD ground temperature was recorded at the initially measured ice table depth of 0.3 m with a mean of 1.9 ^∘C. We decided to fill gaps in the ground temperatures at ice table depth for NWT with this mean value, which was the case for the entire record at GT. Thereby we avoided unphysical spikes in the calculated ΔS_TL for NWT, which would have resulted from setting the temperature for these cases to 0 ^∘C at ice table depth.

Despite the simplified assumption of constant ice table depth in SEB calculations and a deviation of temperature at ice table depth from 0 ^∘C for PLD, we argue that the possible error in ΔS_TL is small compared to the flux magnitudes (see Appendix A) since most of the energy storage change occurred close to the ground surface. The computed ΔS_TL was discarded when calculated from temperature measurements at one depth only, which applied to 34 % of the data for WT.

The second component of the ground heat flux, Q_IT, was approximated as the SEB residual after quantifying ${\mathit{Q}}_{\mathrm{S}}^{*}$, Q_H, Q_LE, and ΔS_TL according to the methods laid out above. Estimating Q_IT as a residual assumes the SEB to be closed and therefore poses limitations to its interpretation as the conductive energy transported away from the surface into the ground, available to deepen the thawed layer by moving the ice table depth further into the ground toward the depth of the permafrost layer as the summer season progresses.

Almost all experimental SEB studies in non-permafrost environments using direct flux measurements from eddy covariance have reported a significant residual on the order of 10 % to 30 % of the net radiation ${\mathit{Q}}_{\mathrm{S}}^{*}$ (see, e.g., list in Foken, 2008) with a few rare exceptions (e.g., Mauder and Foken, 2006). The proposed explanations for the observed systematic residual are manifold and include, after compensating for measurement artifacts and impacts of post-field data processing, surface heterogeneity creating semi-stationary convective structures which are systematically neglected in near-surface observations. We intentionally did not consider applying any of the common methods to eliminate the SEB residual as the results from non-permafrost surfaces may not be applicable to permafrost-dominated ecosystems including the MDV.

While the discussion of the potential causes for the SEB residual is outside of the scope of our study, its implications for estimating Q_IT need to be considered. Any systematic observational residual ϵ≥0 W m⁻² would systematically enhance Q_IT as $\mathit{ϵ}+{\mathit{Q}}_{\mathrm{IT}}=-{\mathit{Q}}_{\mathrm{S}}^{*}-{\mathit{Q}}_{\mathrm{H}}-{\mathit{Q}}_{\mathrm{LE}}-\mathrm{\Delta }{S}_{\mathrm{TL}}$, potentially even beyond physically reasonable limits, and lead to a much deeper thawed layer. We recall that for our MDV sites the energy flux directed toward deepening the melt front may be the dominant term or of similar magnitude compared to the sensible and evaporative heat fluxes. Hence, significant errors should not go undetected and some certainty if estimating Q_IT as the SEB residual should be gleaned from comparing the magnitude of all SEB components, their signs and their relative timing, and by comparing the observed DTR to that computed from Q_IT assuming a certain ice concentration in the soil. This evaluation is included in the discussion below.

To date, there is no experimental technique based on direct flux measurements which can eliminate the observational residual with certainty. The residual was found to be greatly reduced when using indirect techniques such as large-aperture scintillometry based upon similarity theory or large-eddy simulation modeling to estimate SEB components (Foken et al., 2010). For comparison, any SEB model would by definition assume the SEB to be closed, and models are routinely constrained by observations potentially suffering from the observational residual. We therefore do not consider our approach to estimate Q_IT a flaw, but rather an assumption which needs to be evaluated against the background of the results and our understanding of the physical heat transfer and melting processes.

The expected lag of the diurnal temperature signal t_lag penetrating from the surface ground to the ice table depth z_IT was calculated after Woo (2012) as

with period L=1 d and thermal diffusivity D_G.

All SEB components were filtered for physically implausible values. The resulting data gaps corresponded to between 2 % and 34 % of available data for SEB components. The gaps were filled by linear interpolation whenever possible, leading to a data coverage between 66 % and 100 % for WT and 93 % to 96 % for NWT (see also Supplement Fig. S2).

2.4 Identifying turbulent heat fluxes from the water track via footprint modeling

The time-variant footprint of the turbulent energy and mass fluxes was modeled for WT with the Lagrangian-statistics backward TERRAFEX model of the University of Bayreuth (Göckede, 2001). The land cover matrix needed for footprint modeling was generated by mapping the main land cover classes on a QuickBird satellite image. The used land cover classes include modern stream channels, water tracks, paleolake delta sediments, glacial tills, and exposed ice. Pixel classification was done on the basis of albedo, texture, and field descriptions. Water tracks were mapped as multi-segment lines and were assigned a constant width of 10 m derived from the in situ visual observations. Stream channels were mapped as multi-segment lines and were assigned a fixed width of 20 m. Landscape regions were given a single classification code avoiding any overlap. Features were mapped to provide continuous plan-view coverage with no gaps between features and no unassigned cells. Vector landscape features were rasterized at 10 m px⁻¹ and were exported into a local Lambert conformal conic projection to produce gridded land cover values. Each class was assigned an adequate momentum surface roughness length from literature and used in footprint calculations.

Since flux footprint modeling can provide an accurate estimate of the shape, size, and orientation of the source area, we selected the turbulent sensible and latent evaporative heat fluxes of the water track at WT if a cumulative percentage of at least ≈80 % was contributed by water track pixels as indicated by the land cover map, which corresponds to ≈50 % of the footprint spatial extent (Fig. 3). As there was a good match between the size of flux source area and the surface area of the water track, this spatiotemporal filtering is expected to isolate the contribution of the water track at WT to the SEB in a meaningful fashion. Lending to the significant insolation throughout the diurnal cycle in the Antarctic summer in combination with the low albedo of the MDV surface, net radiation was consistently negative, showing a net energy gain at the surface. Hence, no dynamically stable atmospheric conditions were observed which would have resulted in increased flux footprint sizes and complicated selection of meaningful heat fluxes dominated by water tracks. Therefore, no further filtering was necessary to avoid errors in turbulent heat fluxes caused by weak turbulence in the case of dynamic stability.

Figure 3Excerpt of the land cover matrix for lower Taylor Valley with eddy-covariance station locations on water track (WT), paleolake delta (PLD), and glacial till (GT) surfaces. Extent of the flux footprint of WT is indicated by three probability density lines.

Download

2.5 Statistics

Average diurnal variations were calculated for all SEB components with all representative data for comparing individual energy fluxes between WT and NWT (see Supplement Fig. S2). We used linear regression models to compare magnitudes of these average diurnal variations between WT and NWT.

For SEB evaluation average diurnal variations were only calculated from data where all SEB components were available, which was fulfilled for 36 % of the data set (see Supplement Fig. S2). Average energy fluxes for the complete measuring period were computed by averaging mean diurnal variations in the SEB since the data coverage for SEB was unevenly distributed across the diurnal course.

We estimated average energy fluxes for each SEB component Q_x spatially integrated for bare soil in the entire lower Taylor Valley with

where r_WT, r_PLD, and r_GT denote the ratio of bare soil cells in the land cover matrix covered by water tracks, paleolake deltas, and glacial till, respectively, and Q_WT and Q_NWT signify energy fluxes at WT and NWT, respectively.

3.1 Net radiation

Magnitudes and average fluxes of ${\mathit{Q}}_{\mathrm{S}}^{*}$ were increased by a factor of 1.1 at WT relative to NWT (Fig. 4) determined from a linear correlation (R²=0.99, p<0.001). A lower albedo was observed at WT (0.13±0.01) compared to NWT (0.16±0.01), averaged over radiation measurements between 11:00 and 13:00 UTC+12. We interpret the increased ${\mathit{Q}}_{\mathrm{S}}^{*}$ at WT as being partly caused by the reduced albedo (Levy et al., 2011). Levy et al. (2014) found even greater differences in albedo of 0.22 for on-track and 0.14 for off-track soils, suggesting that the difference in energy uptake between water tracks and dry soils may even be greater at other locations.

Lower radiative surface temperatures were observed for WT (5.0±3.2 ^∘C) than for NWT (6.3±5.1 ^∘C) when averaging over the entire record. Differences in surface temperatures between the stations were lowest at low solar angles and greatest at high solar angles, where temperatures reached up to 12.7 ^∘C for WT and 19.5 ^∘C for NWT. These observations suggest that increased ${\mathit{Q}}_{\mathrm{S}}^{*}$ at WT relative to NWT can be partly explained by the reduced energy loss through upwelling longwave radiation due to lower surface temperatures, which are caused by both higher evaporative cooling and higher soil heat capacity (Table 2). This reduced surface heating at WT explains why the differences in ${\mathit{Q}}_{\mathrm{S}}^{*}$ between WT and NWT are greatest around solar noon when solar zenith angles are smallest (Fig. 5).

Table 2Measurements of ice table depth as depth to refusal DTR (m) and soil thermal properties close to the surface: soil volumetric heat capacity C_G ($\mathrm{MJ}{\mathrm{m}}^{-\mathrm{3}}{\mathrm{K}}^{-\mathrm{1}}$), thermal conductivity K_G ($\mathrm{W}{\mathrm{m}}^{-\mathrm{1}}{\mathrm{K}}^{-\mathrm{1}}$), thermal diffusivity D_G (mm² s⁻¹); mean ± standard deviation. Measurements were taken for WT (water track) and PLD (paleolake delta).

Download Print Version | Download XLSX

Figure 5Comparison of diurnal variations in inverted net radiation $-{\mathit{Q}}_{\mathrm{S}}^{*}$ (a), latent evaporative heat flux Q_LE (b), sensible heat flux Q_H (c), ground heat flux Q_G (d), energy storage change in the thawed layer ΔS_TL (e), and soil heat flux at ice table depth Q_IT (f) between WT (water track) and NWT (non-water track), averaged over the entire measuring period of 26 d. Energy fluxes for (a) and negative energy fluxes for (d) are directed towards the ground surface; energy fluxes for (b)–(d) and (f) and positive energy fluxes for (d) point away from the ground surface. Positive energy fluxes for (e) are directed from the ground surface or the ice table into the thawed layer and negative fluxes move out of the thawed layer towards the ground surface or the ice table. Error bars depict standard errors of the mean.

Download

3.2 Turbulent heat fluxes

Averages of the combined turbulent heat fluxes Q_H+Q_LE reached ≈130 W m⁻² at both WT and NWT (Fig. 4b). The partitioning between Q_H and Q_LE, however, strongly differed: at WT magnitudes of Q_LE computed from linear correlation (R²=0.84, p<0.001) were increased by a factor of 3.0 relative to NWT, while magnitudes of Q_H were reduced to 0.7 (R²=0.97, p<0.001). Average heat fluxes at WT were increased by a factor of 3.2 for Q_LE and reduced to 0.8 for Q_H compared to NWT (Fig. 4b). Mean Bowen ratio $Bo={\mathit{Q}}_{\mathrm{H}}{\mathit{Q}}_{\mathrm{LE}}^{-\mathrm{1}}$ over the measuring period was 3.3 for WT and 15.5 for NWT. Q_LE peaked about 1 h earlier during the day and Q_H 1 h later at WT compared to NWT. Turbulent heat fluxes reached their daily maximum in the afternoon irrespective of the surface (Fig. 5).

Increased water content of the water track soil (Langford et al., 2015; Levy et al., 2011, 2014) was likely the cause for increased Q_LE and reduced Q_H, confirming the findings in high-arctic environments (e.g., Westermann et al., 2009). The high soil moisture also led to the earlier, longer, and more pronounced daily peak of Q_LE at WT due to the greater heat diffusivity. The resulting strong evaporative cooling constrained Q_H by weakening surface warming at WT. This finding likely explains why Q_H showed a later and less distinct daily maximum at WT compared to NWT, where energy is preferentially exchanged through Q_H. The maximum turbulence intensity was observed in the late afternoon hours, which led to the concurrent peak in turbulent fluxes.

3.3 Conductive heat fluxes

Mean soil temperatures in the topmost 0.3 m of the soil were lower for WT (2.9±2.0 ^∘C) than for NWT (3.7±2.6 ^∘C) (see also Supplement Fig. S4), contradicting previous observations that soil temperatures of water tracks are higher than their surroundings in summer (Levy et al., 2014). The DTR measurements yielded greater ice table depths for WT than for PLD (Table 2), which confirmed contrasting active layer depths between water tracks (0.45 m) and dry adjacent soils (0.19 m) observed in Lake Hoare and Lake Bonney basins (Levy et al., 2011; Gooseff et al., 2013). Modeled active layer depths for lower Taylor Valley of 0.45 to 0.75 m (Bockheim et al., 2007) indicate that in our experiment the active layer was not fully thawed and the ice table was still being lowered.

Linear correlation analysis of the conductive heat fluxes between WT and NWT yielded an increase by a factor of 1.4 for Q_G (R²=0.93, p<0.001) and by a factor of 1.8 for ΔS_TL (R²=0.92, p<0.001) at WT. Average diurnal variations in Q_IT showed a substantial scatter preventing a meaningful correlation. Average heat fluxes at WT were increased by a factor of 1.3 for Q_G and Q_IT relative to NWT while they were zero for ΔS_TL (Fig. 4). The distinctly higher ΔS_TL at WT compared to NWT at the time of the daily maximum can be explained by the higher thermal diffusivity of the water track soil compared to dry paleolake delta sediments (Table 2) due to increased soil moisture mentioned earlier (Ikard et al., 2009; Levy et al., 2011).

On average, Q_G was directed into the soil except for the period of low solar elevation angles between 21:30 and 03:30 UTC+12 when the net energy transported was directed out of the ground to the surface. ΔS_TL changed sign approximately every 12 h, turning the thawed layer from an energy sink into a source and back. Q_IT remained positive throughout the day when averaged over the entire observational period. The simultaneous occurrence of negative Q_G and positive Q_IT at low solar elevation angles reveals a flux divergence in ΔS_TL, which was directed towards both the ground surface and the ice table at these times, leading to significant net energy loss from the buffering thawed layer.

Since the active layer is typically not fully thawed before mid-January in the MDV (Adlam et al., 2010; Conovitz et al., 2006), we argue that the ice table was being lowered throughout the recording period and most of Q_IT was consumed by latent heat thawing the frozen layer, which explains the large energy fluxes matching soil heat flux observations in other permafrost regions (Lloyd et al., 2001; Lund et al., 2014; Westermann et al., 2009). The increased Q_IT at WT relative to NWT was likely due to a larger amount of ice in the frozen layer and stronger transport of energy through the thermally more conductive thawed layer.

The mean diurnal peak in Q_IT occurred 2.5 h earlier for WT compared to NWT (Fig. 5). Calculation of the diurnal temperature wave penetration lag t_lag from the surface to ice table depth with Eq. (9) yielding a 15 and 13 h lag for WT and NWT, respectively. This finding suggests that the observed daily peak in Q_IT at one specific day corresponded to the peak in Q_G from the previous day delayed by 21 and 24 h for WT and NWT, respectively. The relatively earlier peak of Q_IT at WT lends strong support to our claim of a faster transport of energy through the thawed layer due to its increased thermal diffusivity in water tracks (Table 2) despite the greater depth of the thawed layer.

As mentioned above, our findings of Q_IT need to be carefully discussed in the context of potential artifacts arising from a nonzero observational residual in the SEB, which could obscure true physical behavior of this significant SEB term in permafrost-dominated ecosystems. On the one hand, the increased scatter in the correlation analysis between WT and NWT points to some random error in Q_IT estimated as the residual of the SEB. Random error is largest for the residual term as the random uncertainty of all other directly measured or observation-based computed terms is projected into variations in Q_IT, while the individual random errors may not compensate for each other for an individual averaging period and across the observational period.

On the other hand, many findings support the claim that our estimates of Q_IT represent the true physical heat transport and ice melting process: (i) Q_IT was increased by 30 % at WT compared to NWT, which over the course of the summer melting season should lead to a deeper thawed layer and lower ice table depth, both confirmed by our DTR measurements. (ii) The sign of Q_IT was always positive and its magnitude was similar to that of $({\mathit{Q}}_{\mathrm{H}}+{\mathit{Q}}_{\mathrm{LE}})/\mathrm{2}$ confirming that a significant fraction of net radiation is consumed for ice table lowering by melting in the peak summer season. The magnitudes of our Q_IT estimates were similar to those found in other permafrost landscapes. If Q_IT were dominated by ϵ, then the residual would correspond to roughly 30 % of ${\mathit{Q}}_{\mathrm{S}}^{*}$, which is the maximum observed value for non-permafrost surfaces. Given the relative homogeneity in surface conditions in MDV compared to that in vegetated non-permafrost surfaces, it is unsure whether such a large residual would be realistic. (iii) The diurnal course of Q_IT showing an earlier peak for WT matched our expectation from the time lag analysis in combination with the soil thermal property measurements. Those would lead to a small dominance of the effect of enhanced thermal diffusivity over the greater transport time required for the daily heat pulse from the ground surface to penetrate to the ice table depth due to the greater distance to be covered. One must recall that the thermal diffusivity is the ratio of the energy conducted to that stored in a given medium. (iv) The oscillatory nature of the measured ΔS_TL over the course of the day equating to zero when averaged over the ensemble averages at any location lends strong support to our interpretation as the temporary heat storage reservoir and thus heat communicational layer. We claim that when it is combined with the other measured SEB terms in the form of $-{\mathit{Q}}_{\mathrm{S}}^{*}-{\mathit{Q}}_{\mathrm{H}}-{\mathit{Q}}_{\mathrm{LE}}-\mathrm{\Delta }{S}_{\mathrm{TL}}$ it forces physically meaningful diurnal dynamics in Q_G and Q_IT.

However, the existence of some influence of a true SEB observational residual ϵ was confirmed, but judging from our observations it induces a random component in the diurnal course of the ensemble-averaged components and is not expected to have a systematic diurnal signal leading to an unphysical timing or magnitude in Q_IT. Our method is unable to rule out the existence of a constant, time-invariant systematic ϵ, which would effect changes in the magnitude of Q_IT. However, we would expect its magnitude and thus its impact to be small. An independent verification of our method to estimate Q_IT as the SEB residual could have been gleaned from observations extending into the after-summer season when the ice table depth is decreasing due to the refreezing of the soil moisture. Under those conditions, a change in sign of Q_IT would be expected, which could lend support to a physically meaningful quantification of all SEB terms. Unfortunately, we were unable to extend the observational period for logistical reasons.

3.4 Climate change response of soils

How will expected 21st century warming in the MDV (Arblaster and Meehl, 2006; Chapman and Walsh, 2007) impact mass and energy fluxes of soils? Though water tracks and other wetted soils show high spatial consistency in the MDV, their spatial extent is largely dependent on snow and ground ice melt (Gooseff et al., 2013; Langford et al., 2015; Levy, 2015).

Snow abundance, distribution, and duration will be an important factor for predicting changes to wetted area since snow is a major water source for many water tracks in the MDV (Langford et al., 2015). Snowfall is projected to rise during the next century (Christensen et al., 2013), and snowdrift accumulation, providing an amount of snow equal to snowfall in the valley floors, will mainly depend on dynamics in katabatic winds delivering snow from the valley walls and ice sheet (Fountain et al., 2009). Melt of existing snow patches in the summer is likely to increase in response to a warming climate. As water tracks also feed on ground ice melt (Harris et al., 2007), deepening active layer thaw is all but guaranteed to produce additional melt. This, in turn, may drive subsidence, which could trap new snow banks, encouraging further water track activity. Therefore, we expect water tracks to expand in surface fraction cover and activity in response to climate change: (i) as active layers thaw more deeply and melt previously stable ground ice due to positive feedbacks of increased thermal diffusivity and energy uptake of wet soils (Gooseff et al., 2013; Ikard et al., 2009; Levy and Schmidt, 2016), (ii) as perennial snow patches melt and thin during summer warming, and (iii) as increased evaporation from the Ross Sea and decreased sea ice raises humidity, allowing enhanced soil salt hydration. All these mechanisms collectively or individually point to the potential for expanded surface fractions and connectivity of wetted soil in the MDV (Ball et al., 2011; Wall, 2007). This anticipated development of water track abundance is contrary to projections for the high arctic where primarily snow-fed water tracks are likely to be reduced in area because precipitation will increasingly occur in liquid form (Comte et al., 2018).

To provide some quantification of the anticipated response of energy and matter fluxes in the MDV in the summer season based upon our evaluation of the SEB of water tracks, we calculated the increase in average energy fluxes integrated for bare soil in the entire lower Taylor Valley resulting from arbitrary but realistic increases in the relative water track fraction of bare soil using Eq. (10) (Table 3). A doubling scenario is a realistic increase in water track coverage considering that the spatial extent of wetted soils in the MDV varies by at least a factor of 2 between cold and warm years (Langford et al., 2015). This doubling would lead to an increase in evaporation from bare soil in the entire lower Taylor Valley by 6 % to 0.36 mm d⁻¹, which implies a significant increase in evaporation in the MDV given the expected response mechanisms in soils mentioned above.

Table 3Integrated energy fluxes from bare soil in lower Taylor Valley calculated with Eq. (10). “Factor” denotes different increments of water track coverage (WTC) relative to the observed value of 2.9 %. Net radiation (${\mathit{Q}}_{\mathrm{S}}^{*}$), sensible heat flux (Q_H), latent evaporative heat flux (Q_LE), and ground heat flux (Q_G) in watts per square meter.

Download Print Version | Download XLSX

One likely implication of our results demonstrating the intensified heat cycling for water tracks is that internal physical processes of water tracks, i.e., energy exchange and hydrology, could respond faster and more vigorously to climate change than dry off-track soils in addition to an increase in coverage in the MDV. Since increased thermal conductivity and energy uptake enhance positive soil thawing feedbacks for water tracks (Gooseff et al., 2013; Ikard et al., 2009; Levy and Schmidt, 2016), soil moisture may rise and active layers deepen, leading to increases in discharge and evaporation compared to the current conditions. These effects can be assessed in future monitoring efforts for water tracks, where observed trends in internal processes of water tracks would serve as indicators for landscape-scale climate change response in the MDV.

3.5 Climate change scenarios

Since we have presented many indications of increasing spatial extent of water tracks with climate change in the MDV, a discussion of possible pathways of climate change and their respective roles in determining the change in spatiotemporal dynamics of water tracks will benefit the framing of energy and mass flux scenarios for the MDV, and could serve as a basis for advanced microscale modeling to provide constraints on the timing and spatial feedbacks of these predicted mechanisms. Here we briefly consider two contrasting scenarios informed by meteorological first principles: one where insolation increases but regional air temperatures remain constant (Gooseff et al., 2017), and one where temperatures increase but insolation remains constant, which corresponds to a more typical Arctic-type climatology.

While we lack sufficient information to make quantitative estimates for the two alternative scenarios, we can offer predictions based upon mechanistic reasoning: in the first scenario of constant insolation and rising temperatures in the MDV, an increase in vapor pressure deficit would lead to delayed cloud formation while net radiation would increase because of reduced cloud coverage, causing a positive feedback on temperature. The increased temperature would act to intensify snowmelt and soil thawing. This effect would even be enhanced by the occurrence of strong down-valley winds which can increase temperatures and melting rates in the summer dramatically (Doran et al., 2008). In the second alternate scenario temperatures stay constant and insolation increases as observed by Gooseff et al. (2017). Higher net radiation will cause an evaporation increase which will reduce the vapor pressure deficit since air temperatures remain constant. The increased cloud formation would then cause a negative feedback on insolation and lead to a weaker increase in snow and ground ice melting than in the first scenario. Both scenarios represent pathways of climate change which would increase spatial extent of water tracks and other wetted soils, though the magnitude of this climate change response may vary between them.

3.6 Implications for soil ecosystems

Here, we speculate on the implications of the differences in SEB found between water tracks and dry soils for the biologic communities and geochemistry: since water tracks are wetter than adjacent off-track soils at the surface, it is possible that these conditions promote enhanced chemical weathering of the water track soils (Campbell et al., 1998). One might expect that with expanded water tracks, solute fluxes to water tracks would increase as chemical weathering initiates on what were cold, dry soils. Increased weathering could also result in increased phosphorous fluxes to soils, enhancing microbial and invertebrate habitat suitability (Heindel et al., 2017). Our SEB findings suggest that water tracks go through fewer freeze–thaw cycles than adjacent off-track soils, due to the release of latent heat during freezing. This reduces biotic stresses, which in turn might improve habitat suitability for both microbes and invertebrates (Yergeau and Kowalchuk, 2008; Convey, 1996). Active layer depths are typically several decimeters deeper in water tracks than they are in adjacent dry soils (Levy et al., 2011). That suggests that occupation of currently dry soils by water tracks may result in an expansion of the seasonally thawed portion of the soil column, producing expanded habitat for soil invertebrates.

However, these effects may only lead to increase in biomass, productivity, and respiration if episodic overland flow flushes solutes out of the wetted soil (Ball and Levy, 2015; Zeglin et al., 2009). Overland flow can only occur when meltwater availability exceeds infiltration capacity and hydrologic conductivity of water track soils (Levy, 2015). In combination, the discussed processes could lead to more hydrologically connected and less saline soil habitats at increased and more uniform soil moisture in the MDV and, by that, increased biogeochemical cycling (Ball et al., 2011; Gooseff et al., 2013). This habitat change could suppress dry-tolerant taxa like the dominant nematode Scottnema lindsayae and favor wet-tolerant generalist taxa at lower biodiversity levels (Buelow et al., 2016; Simmons et al., 2009; Wall, 2007). However, it is likely that overall improved habitat suitability would increase biological productivity and respiration in soils in the MDV.