Climate and energy balance of the ground in University Valley, Antarctica

Abstract

We report 3 years of data from one meteorological and three smaller stations in University Valley, a high-elevation (1677 m) site in the Dry Valleys of Antarctica with extensive dry permafrost. Mean air temperature was -23.4°C. Summer air temperatures were virtually always < 0°C and were consistent with the altitude lapse rate and empirical relationships between summer temperature, distance from the coast and elevation. The measured frost point (-22.5°C) at the 42 cm deep ice table is equal to the surface frost point and above the atmospheric frost point (-29.6°C), providing direct evidence that surface conditions control ground ice depth. Observed peak surface soil temperatures reach 6°C for ice-cemented ground > 15 cm deep but stay < 0°C when it is shallower. We develop an energy balance model tuned to this rocky and dry environment. We find that differences in peak soil surface temperatures are primarily due to the higher thermal diffusivity of ice-cemented ground compared to dry soil. Sensitivity studies show that expected natural variability is insufficient for melt to form and significant excursions from current conditions are required. The site's ice table meets the criteria for a Special Region on Mars, with 30% of the year > -18°C and water activity > 0.6.

Introduction

Dry permafrost, defined as ground that never warms > 0°C and has a negligible ice content, is rare on Earth. In the high elevations of the Antarctic Dry Valleys, there is dry permafrost with no ice below (Campbell & Claridge 2006, Bockheim et al. 2007), as well as dry permafrost overlying ice-cemented ground (McKay et al. 1998, Campbell & Claridge 2006, Bockheim et al. 2007). Year-round atmospheric and subsurface temperature measurements by McKay et al. (1998) at Linnaeus Terrace in Upper Wright Valley at an elevation of 1600–1650 m indicated dry permafrost extending from 12.5 cm below the surface to the top of the ice-cemented ground at 45 cm. Dry permafrost over ice-cemented ground has also been reported near Mount Dolence in Ellsworth Land (Schaefer et al. 2017, McKay et al. 2019).

Dry permafrost results when the average frost point at the surface is lower than the average frost point of ice-cemented ground below the surface while the temperature remains < 0°C year-round. This causes any ice that is already present to sublime into the atmosphere and retreat deeper, leaving a dry permafrost soil layer; or conversely ice cannot deposit into these soils as it is unstable. The frost point decreases with depth (McKay 2009), and the ice table will stabilize at a deeper level. Initial assessments of the stability of the ice table in the upper Dry Valleys assumed that the atmospheric boundary condition determined the stability of the ice. However, recent work has shown that the effective frost point of the top of the soil column can be different from the atmosphere due to the presence of frost and snow or due to increased moisture content. The temperature and moisture conditions at the surface determine the depth and stability of the ice; the surface conditions are controlled by but not equal to the atmospheric values (Hagedorn et al. 2007, McKay 2009, Liu et al. 2015, Fisher et al. 2016, McKay et al. 2019).

University Valley (77°52'S, 160°45'E; ~1700 m above sea level), one of the upper valleys in the Quatermain Range in the Dry Valleys, is of particular interest in the study of dry permafrost because observations there indicate that the average depth to ice-cemented ground increases linearly across the length of the valley floor (~1.5 km) from near zero to over 1 m (McKay 2009, Marinova et al. 2013). This provides a natural experiment regarding the environmental conditions that determine the ice content and temperature regime of dry permafrost. The conditions that allow for ice-cemented ground to reach the surface in these arid and dry locations are not well understood.

In this paper, we report on the climate conditions at University Valley as determined by 3 continuous years of observations of the atmosphere, surface and subsurface dry permafrost and underlying ice-cemented ground. The data were collected from one full meteorological station and three smaller permafrost stations, with the four data collection sites distributed in the Valley (Fig. 1) and having different depths to ice-cemented ground. We also compare this site to other locations in the Dry Valleys for which extensive environmental data exist.

Fig. 1. a. Location of University Valley, Antarctica. b. Aerial photograph of University Valley and the locations of the meteorological weather station (W) and smaller stations over locations of deep (D), medium-depth (M) and shallow (S) ice-cemented ground; yellow bars show elevations in metres. North is up. c. The meteorological station setup.

Lacelle et al. (2016) incorporated this University Valley dataset in a comparison of the relations between solar radiation and air and ground temperatures in University Valley and nearby Beacon Valley with those in ice-free Victoria Land and Arctic Canada, focusing their analysis on thermal offset factors.

In this analysis, we develop a detailed energy balance model and use the observational data to constrain the model parameters. Key energy balance model parameters such as albedo and emissivity, but especially thermal diffusivity and roughness length scales, are not well known for the high-elevation Dry Valleys, and this modelling approach develops a way to predict and study surface and subsurface conditions based on limited (past or future) datasets. The model uses atmospheric (temperature, humidity, wind) and solar insulation inputs to determine the surface and subsurface temperature profiles in University Valley. The model can be validated by comparing its results to the measured ground temperatures. The model is used to explore the sensitivity of surface and subsurface conditions to varying environmental properties, corresponding with the natural variability observed within a single valley and at nearby locations. In addition, we examine the changes to the environmental conditions and parameters that would be required to reach key inflection points where the environmental regime changes, such as temperatures above freezing occurring and allowing the presence of liquid water. The presence of liquid water specifically has key implications for habitability and weathering processes.

The model is analogous to that developed by Hunt et al. (2010) for the lower elevations in the Dry Valleys. Hunt et al. (2010) developed an energy balance model for soil temperature and water percolation on the valley floor of Taylor Valley patterned after standard approaches for bare soils without significant vegetation (Matthias 1990). The water balance included snowmelt, freezing/thawing of soil water, soil capillary flow and vapour flows. By computing the surface fluxes of sensible heat exchange between the atmosphere and ground surface, subsurface heat conduction and shortwave and longwave radiation, the model could account for 96–99% of the variation in soil temperature. A significant simplification compared to the model of Hunt et al. (2010) is that, unlike the floor of Taylor Valley, liquid water movement or freezing in the soil column is not important at University Valley. Ice and water vapour are present and exchange by sublimation and condensation (e.g. Lacelle et al. 2013). A complexity of our problem compared to that of Hunt et al. (2010) is determining the properties of the dry permafrost layer, in particular the thermal diffusivity as the temperature and moisture content vary.

University Valley is of special interest as an analogue for Mars environmental and habitability studies. While dry permafrost is rare on Earth, it is widespread on Mars, where dry permafrost begins at the surface and is underlain by ice-cemented ground (Mellon & Jakosky 1993, Mellon et al. 2009, Smith et al. 2009, Mellon & Sizemore 2021). The Antarctic Dry Valleys represent one of the driest and coldest places on Earth, and as such they provide an interesting comparison point to conditions in the Martian subsurface. We make this comparison based on the temperature and water activity limits defined for Special Regions on Mars (Rummel et al. 2014).

Methods

Meteorological measurements

A full meteorological station and three smaller permafrost stations were deployed in University Valley on 11 December 2009, including air, surface and subsurface instruments. Data were collected through 31 January 2013. The overall site and the specific locations of the stations are shown in Fig. 1. As can be seen in Fig. 1, all four stations are between 1650 and 1700 m elevation on rocky ground on the valley floor.

The meteorological station is at 77°51.729'S, 160°42.606'E, elevation 1677 m, and it is based on a Campbell CR1000 data logger operated by solar-recharged batteries. The air temperature and relative humidity (RH) were measured 1.2 m above the surface with a Campbell 207 probe in a ventilated radiation shield, as can be seen in Fig. 1. Wind was recorded using an RM Young 05103 wind anemometer (wind speed and direction, 1.08 m above ground) and solar insolation was recorded using a LI-COR sky radiation pyranometer (LI200X, 1.27 m above ground). In the (sub)surface, Campbell 107 temperature probes were placed at the surface (covered with a thin layer of soil), at 10 cm depth, at the top boundary of the ice-cemented ground (42 cm depth) and at 7 cm into the ice-cemented ground (49 cm depth from soil surface). Three Onset U23 Pro v2 External Temperature/Relative Humidity loggers were also placed at the meteorological station site to measure temperature and humidity at the surface (covered with a thin layer of soil), at 20 cm depth and at the ice-cemented ground interface at 42 cm depth. All Campbell instruments were sampled every 30 min, while all Onset sensors were sampled every 60 min (due to limitations in memory capacity). The deep, medium-depth and shallow ice-cemented ground stations were placed in locations with ice depths of 36, 22 and 8 cm, respectively. Their respective locations are 77°51.946'S, 160°43.720'E (elevation 1697 m), 77°51.773'S, 160°43.527′E (elevation 1661 m) and 77°51.584′S, 160°43.019′E (elevation 1684 m). Each site has a combination of an Onset 4 channel smart logger (U12-008) and U23 Pro v2 External Temperature/Relative Humidity loggers with sensors at/close to the surface, halfway to the ice-cemented ground and at the ice-cemented ground interface.

The complete processed datasets are archived at the National Snow and Ice Data Center.

The placement of the probes aimed to minimally disturb the soil and ensure good thermal connectivity with the soil. A narrow hole was manually dug at each site. The surface sensors were placed ~30 cm from the hole, slightly depressed into the ground and covered with a thin layer of soil, aiming to ensure that the surface albedo was maintained. As will be discussed later, the removal of this overlying material by wind and exposing the probe surface directly to the sun may contribute to some non-representative temperature measurements. Sensors in the dry permafrost were inserted 5–10 cm into the side of the dug hole, aiming to minimize the soil disturbance. For the probe in the ice-cemented ground, a hole was drilled using a drill bit with the same diameter as the sensor; the sensor was then snuggly inserted into the hole, ensuring good thermal contact.

The operating temperature range for the Campbell system is -50°C to 50°C, which encompasses the experienced temperature range. The error in the temperature measurement of the 207 and the 107 probes is ±0.2°C. The error in the humidity measurement is < 10%. An important caveat is that the Campbell 207 RH sensor has high errors for RH values < 15%, usually tending to systematically overestimate values; the measured humidity, however, rarely drops to < 30%. The LI-COR LI200X is a solar pyranometer that is sensitive between 400 and 1100 nm with an absolute error in natural daylight of ±5% maximum (±3% typical). The RM Young 05013 Wind Monitor has a stated measurement range of 0–100 m s^-1 with a threshold of 1 m s^-1 and an error in speed of ±0.3 m s^-1 or 1% of the value. The range of directions is 0–360°, with an error of ±3°.

For the Campbell data, the weather station battery value stayed within working limits, including through the winter, reaching a low of 12.2 V and increasing to 14.5 V in summer, supporting the integrity of the data. Full datasets, spanning from 5 December 2009 to 31 January 2013, are available for air temperature and humidity, solar insolation, wind speed and direction and temperatures at 10, 42 and 49 cm depths. The surface temperature is only available for 11 December 2009 to 8 December 2010 and from 26 April 2012 to 11 December 2012; this is the temperature sensor used in all cases in this paper when surface temperature is reported.

For the Onset U23 Pro v2, which is a combined humidity and temperature probe, the operating range is -40°C to 70°C. The error estimates for the temperature sensor are listed as ±0.25°C from -40°C to 0°C and ±0.2°C from 0°C to 70°C. The resolution of the temperature data is 0.04°C. Drift is < 0.01°C per year. Because the temperature-sensing element is a 10 kΩ thermistor, we assume that the main source of error is due to systematic offsets and not random noise. Hence, the error in average values is not significantly decreased with increased sample size. The accuracy of the RH sensor is ±2.5% in the 10–90% RH range and ±5% at < 10% or > 90% RH. The resolution of the sensor is 0.05% and its expected drift is < 1% per year. At temperatures < -20°C or > 95% RH, an additional error of 1% may be present.

For the Onset data, battery information is available for the first year of operation and stays within an appropriate operational range of 3.21–3.56 V. Full datasets, spanning from 5 December 2009 to 31 January 2013, are available for 20 cm depth temperature and humidity and 42 cm depth temperature and humidity. For the surface temperature and humidity, data are available from 10 December 2009 to 4 December 2010 plus some additional short periods of data near the end of the recorded period. The surface humidity data measured by this Onset sensor are used in all cases reported here.

The temperature data at the surface and 42 cm depth were recorded by both the Campbell and Onset sensors. For the surface data, during sunny summer times, the Onset sensor recorded a temperature on average 0.19°C lower than the Campbell sensor, while in winter, the Onset sensor recorded a temperature on average 0.14°C higher than the Campbell sensor. The standard deviations were 0.10°C and 0.13°C, respectively. The maximum differences in temperature were notably higher during the sunny summer times. This is attributed to either or both sensors becoming partially or fully uncovered at the surface by wind moving the thin layer of soil that was placed over the sensors and the Campbell probe having a shiny metallic sensor tip while the Onset sensor has a white tape tip. At 42 cm depth, at the ice-cemented ground interface, on average the Onset sensor recorded a temperature 0.24°C higher than the Campbell sensor, with a standard deviation of 0.13°C. The minimum and maximum deviations were -0.17°C and 0.61°C, respectively.

At the deep (36 cm), medium-depth (22 cm) and shallow (8 cm) ice-cemented ground stations, we also used the Onset U12-008 four-channel outdoor/industrial data logger with the Onset TMCx-HD temperature sensors. The U12 logger has a listed operating range of -20°C to 70°C, but previous experience had demonstrated its operation at lower temperatures. The temperature probes themselves have a measurement range of -40°C to 50°C in soil. The TMCx-HD probe used with the U12 logger has an accuracy of ±0.25°C at 20°C, a drift of < 0.1°C per year and a resolution of 0.03°C at 20°C.

The collected data were checked for integrity. The most notable deviation is periods of missing data, the cause of which is not well understood. Additionally, the Campbell data were corrected for ~5 or fewer points out of a total dataset of > 55 000 data points per sensor. These errant values were probably due to cosmic ray strikes. An exception was the surface temperature, which had ~5% of data points removed due to unrealistically high values (e.g. > 20°C, where this cannot be explained by surface property variation such as rocks) and also had extensive stretches of missing data.

For minimum and maximum temperatures and thaw depths, the instantaneous readings are reported. For the modelling, a time step of 1–2 min was used.

Relative humidity correction and drift

The sensor in contact with the ice table provides a useful test of the response of the humidity sensor, as it is expected to show a constant 100% humidity. McKay et al. (2019) used similar sensors to define the correction needed to determine the humidity values at temperatures below freezing. Anderson (1994, 1995) and Koop (2002) assumed that RH sensors record the RH with respect to liquid water even for temperatures < 0°C, and previous work (e.g. Doran et al. 2002, Hagedorn et al. 2007, Andersen et al. 2015, Liu et al. 2015) used this assumption in correcting Antarctic data. The RH data from the University Valley Campbell 207 probe are corrected this way. However, McKay et al. (2019) found that the results from the Onset sensors at the ice table indicated a correction of

(1)$$RH_i = RH_w- 2 - 0.65T$$

where RH_i is the RH over ice, RH_w is the sensor reading (approximately the RH over water) and T is the temperature in degrees Celsius. When corrected using this method, our data for the ice table are consistent with a RH of 100% within approximately ±1% (see Appendix).

The 3 year duration of our dataset with the RH sensors undisturbed allows for an assessment of the drift of the sensors over time. We find that the sensor drift is an increase in the indicated reading of < 1% per year (see Appendix). Nair et al. (2015) tested the response of humidity sensors stored for 7 months at room temperature and tested at RH values from 30% to 70%. They reported errors of 0.4% ± 0.2% RH (n = 5 sensors). If the sensing film is stable, then the primary cause of sensor ageing when exposed to air is the accumulation of contaminant particles on the sensor material, changing its properties and response time (Nair et al. 2015 and references therein). The constant low-temperature environment of the ice table, with little or no airborne particulates, may minimize sensor drift.

Ground energy model

The energy balance model uses as inputs the local environmental conditions and calculates the surface and subsurface temperatures, which can then be compared to the measured values. The model can be forced using flux boundary conditions at the top and bottom or with enforced temperature boundary conditions to calculate the temperature profile in the in-between layers. In the case of a flux forcing at the top, we nominally use the values measured at the site; however, the flux can also be calculated using representative theoretical values (such as expected solar insolation based on the location of the site and air temperature based on an average temperature and idealized diurnal and yearly fluctuations). The values used for all constants are described here and in more detail later.

The flux of heat at the surface is given by:

(2)$$F_{surf} = Solar + Longwave + Sensible + Latent-Blackbody$$

where F_surf is the net flux at the surface (W m^-2), Solar is the net incident solar radiation, Longwave is the net longwave sky radiation, Sensible heating is due to turbulent heat transfer between the surface and air, Latent heating is due to turbulent heat flux from the movement of water vapour between the surface and air and Blackbody is the blackbody radiation from the surface. We set downwelling radiation (i.e. energy that is added to the surface) as positive, except in the case of Blackbody, which is always positive and denotes energy radiated from the surface. In solving for the subsurface temperatures, F_surf sets the upper boundary condition, and the lower boundary condition is set to no net flux.

Solar is given by the solar flux (S_surf) measured by the weather station (W m^-2) or as calculated by the sun model and adjusted for the surface albedo, or reflectivity, A, at the site of interest:

(3)$$Solar = ( {1-A} ) S_{surf}$$

Longwave is given by (as used by Mölg & Hardy 2004):

(4)$$Longwave = \sigma T_{air}^4 ( {C_1 + C_2E_{air}} ) ( {1 -shadow} ) $$

where σ is the Stefan-Boltzmann constant (5.67 × 10^-8 W m^-2 K^-4), T_air is the air temperature (K), C₁ and C₂ are fitting parameters, which are set to 0.585 and 6.2 × 10^-4 Pa^-1, respectively, following Mölg & Hardy (2004), and E_air is the water vapour partial pressure (Pa). A shadowing term, shadow, takes into account that the surrounding valley walls limit the surface from seeing the full sky.

The sensible heat flux, Sensible (W m^-2), is given by:

(6)$$Sensible = C_{\,p, air}\rho _0\displaystyle{P \over {P_0}}\displaystyle{{K^2v( {T_{air}\;-\;T_{surf}} ) } \over {ln\left({\displaystyle{{z_m} \over {z_{0m}}}} \right)ln\left({\displaystyle{{z_h} \over {z_{0h}}}} \right)}}$$

where C_p,air is the heat capacity of air (J kg^-1 K^-1), ρ₀ is the density of air at standard temperature and pressure (kg m^-3), P is the local air pressure (Pa), P₀ is the air pressure at STP (Pa), K is the von Kármán constant, v is the wind velocity (m s^-1), which is taken from the measured wind velocity at the weather station, T_air is the air temperature (K), T_surf is the surface temperature (K), z_m is the height at which the wind measurement is taken (m), z_0m is the momentum roughness length scale (m), z_h is the height of the air temperature measurement (m) and z_0h is the roughness length scale of temperature (m). The atmospheric pressure is assumed (not measured) to be constant and proportional to the altitude: P = 82.8 kPa for University Valley.

The latent heat flux, Latent (W m^-2), is calculated using:

(7)$$Latent = 0.623L_{sublime}\rho _0\displaystyle{1 \over {P_0}}\displaystyle{{K^2v( {E_{air}\;-\;E_{surf}} ) } \over {ln\left({\displaystyle{{z_m} \over {z_{0m}}}} \right)ln\left({\displaystyle{{z_v} \over {z_{0v}}}} \right)}}$$

where L_sublime is the latent heat of sublimation from the ice in the subsurface to vapour in the atmosphere (J kg^-1), E_air is the vapour pressure in the air (Pa), E_surf is the vapour pressure at the surface (Pa), z_v is the height at which the humidity measurement was taken (m) and z_0v is the roughness length scale of water vapour (m). The vapour pressure in the air is calculated using the measured RH. For the surface, the water vapour pressure is calculated using the measured surface humidity when depth of ice-cemented ground (d_ice) > 20 cm. However, in the case where the ice depth is 3 cm or less, the surface humidity is set to 100%, as the surface is effectively in contact with the ice. For ice depths between 3 and 20 cm, the surface humidity is linearly interpolated so that it equals 100% when d_ice = 3 cm and it equals the measured surface RH when d_ice = 20 cm.

The blackbody radiation is given by:

(8)$$Blackbody = \varepsilon \sigma T_{surf}^4 ( {1 -shadow} ) $$

where ɛ is the emissivity of the surface and T_surf is the surface temperature (K). The shadowing by the surrounding cliffs, shadow, as also used in the equation for Longwave radiation, accounts for the surface not radiating to the cold sky in all 2π steradians due to the surrounding valley walls. The temperature of the surface of the cliffs is assumed to be close to that of the measurement site. Only in the case of Blackbody does a positive value denote energy being lost by the surface; note that Solar and Blackbody always have positive or null values. For all other terms in the energy balance equation, positive values denote energy being deposited into the ground.

When fitting parameters or determining how well a certain modelling run fits the data, we use as input the measured solar flux, calculate the longwave radiation from the measured temperature and air humidity, calculate the sensible and latent heat fluxes using the measured air properties and the calculated surface properties and calculate the blackbody radiation from the calculated surface properties.

Numerical method

Below the surface, we solve the thermal diffusion equation (Eq. (13) in the Appendix) using the Crank-Nicolson method (Acton 1970). This numerical method is a combination of the forward Euler and the backward Euler methods, is implicit in both the time and spatial variables and is generally stable for diffusion equations. The radiative, convective and conductive fluxes together specify the upper boundary condition. The lower boundary condition is a zero-flux condition (∂T/∂x = 0). We can also set the upper or lower boundary condition to specified temperature values over time, which can be set to the measured values. The key parameter that combines the physical properties of the soil and the numerical spacing of the time and spatial variables is α:

(9)$$\alpha = \displaystyle{1 \over 2}\displaystyle{k \over {\rho C_p}}\displaystyle{{\Delta t} \over {\;\Delta x^2}}$$

where k is the thermal conductivity (W m^-1 K^-1), ρ is the density (kg m^-3), C_p is the heat capacity (J kg^-1 K^-1), Δt is the time step that is used (s) and Δx is the thickness of each of the modelled subsurface layers (m). The values of k, ρ and C_p are specific to each layer and can change with depth in the simulated subsurface. We thus simulate the dry soil using the designations of k_soil, ρ_soil and C_p,soil and the ice-cemented ground using the designations k_icy-soil, ρ_icy-soil and C_p,icy-soil. The smaller the value of α, the more stable the result. The time step used is nominally 120 s, although in some cases it is decreased to 60 s to reduce the possibility of the model developing numerical instabilities. The layer thickness was set to 1 cm.

The model can also be run by enforcing temperature boundary conditions; that is, forcing the top and bottom boundary conditions with specified temperatures rather than setting a flux condition. The model then solves for the in-between temperatures. This approached was used in trying to solve for the soil thermal diffusivity, as will be described below.

Parameters used in the model

One of the benefits of collecting in situ data is the ability to extract or at least confirm that the parameters being used in associated models reasonably represent the environment. The ability to do this is determined by limitations in the collected dataset (number of sensors available, disturbing the site during placement and sensor capabilities and failures) and the challenge of using one or a few measurement locations to meaningfully describe and characterize the complex and varying natural environment. In choosing the environmental parameters to represent the University Valley site, we use a multi-pronged approach: attempting to extract these parameters from the available data, reviewing the literature and using other analytical and empirical fits that have been developed. It should be noted that a first attempt was made to do a large parameter space search in order to determine what combination of values would provide the best fit to the data. This approach did not yield useful results for a number of reasons. Primarily, the limitations of the collected data affect the fitting statistics that are calculated. Additionally, the idea of ‘best fit’ is ambiguous, with there being several statistical measures that could be minimized, including mean absolute error, mean error, median error or ensuring no significant overshoots in summer warm periods (the period arguably of most interest). Each of these criteria leads to a different set of fitting parameters that produce a ‘best fit’. Because of ambiguity in which statistical measure is most important to minimize, we chose to use the combination of targeted fitting and environmental parameter values and methods from the literature. The parameters used in the model are summarized in Table I.

Table I. Nominal model parameters and values used.

Surface albedo

The best fit surface albedo is determined by fitting the surface warming rate during sunny, cloudless and low-wind days. The high solar insolation results in the Solar component of the energy balance equation dominate, while the low wind reduces the sensible and latent heat fluxes. Additionally, fitting the shape of the rise in surface temperature specifically focuses on the solar warming rather than requiring an overall accurate fit. The snow albedo is fit using a similar methodology; we infer that the surface is snow-covered when the surface humidity is high (> 90%) and concurrently the atmospheric humidity is low (< 80%). There are few periods appropriate for snow albedo fitting because snow is not common in summer when solar insolation is high.

Using this methodology and fitting three separate summer periods, the surface albedo at the weather station site is determined to be ~0.33. This albedo is overall consistent with published values (e.g. An et al. 2017), although it is towards the higher end. Overall, significant variability in albedo values is apparent in the published data. For snow, we determine the albedo to be 0.86, which is a bit low but close to the range of reported values (e.g. Wiscombe & Warren 1980). However, this value has less of an overall impact as snow is not common during high-insolation and warmer summer periods, which are of most interest for this work.

While we do not have a direct measurement of snow cover at University Valley, we infer that snow-like albedo and emissivity should be used when the measured air humidity is < 80% but the surface humidity is > 90% at the weather station site. It is unclear how well this criterion applies to the dark winter months and how often surface frost rather than snow may produce a similar response in humidities; our modelling and the data cannot provide further insight on this point.

Surface emissivity

The surface emissivity was found by fitting winter surface temperatures as the surface is warmed by katabatic winds but then cools during a period of low wind. The low wind as the surface cools minimizes the latent and sensible heat flux contributions and the winter darkness eliminates the solar insolation term. The energy balance equation reduces to a relationship between the blackbody outgoing radiation, which is a function of surface emissivity and dominates, and the incident longwave radiation. The properties of the subsurface also play a role, as they set the rate at which energy diffuses from or into the surface; however, this term did not appear to appreciably affect the best fit emissivity. We find best fit values of ɛ_soil = 0.92 (from three fit periods), and we chose ɛ_snow = 0.97 to be consistent with the literature (Hunt et al. 2010), as the overall low blackbody signal during the cold winter months made fitting snow emissivity challenging.

Shadowing

The shading at the weather station was measured in 10° radial increments using an inclinometer (valley wall profile shown in Fig. 2). The sky fraction (i.e. the fraction of the hemisphere that is sky) is calculated using the equation for the normalized area of a spherical cap: A_cap = 2π(1 - sinα), where α is the angle from the horizontal to the horizon. For a hemisphere, α = 0°, and this reduces to A_cap = 2π steradians, as expected. To calculate the sky fraction with a valley wall profile that is not at a constant angle, we sum over the measured profile:

(10)$$Sky = \mathop \sum \limits_{0^\circ \;< \;\theta \;< \;360^\circ } 2{\rm \pi }( {1\;-\;\sin \alpha_\theta } ) \left[{\displaystyle{{\Delta \theta } \over {2{\rm \pi }}}} \right]$$

where θ is the radial angle and Δθ is the radial increment over which the angle α was measured. We then get shadow = 1 - Sky. Using the valley wall profiles shown in Fig. 2, we find the shadowing at the weather station site to be 35.5%.

Fig. 2. Profile of the valley walls from the weather station site in University Valley: 35% of the sky is obscured by the valley walls.

Thermal diffusivity

The parameter that appears directly in the diffusion equation (Eq. (13)) is the thermal diffusivity, κ, which depends on fundamental soil properties: ρ (soil density), C (heat capacity) and k (thermal conductivity) as κ = k/ρC. We derive values of the thermal diffusivity considering several different numerical approaches, as discussed by Pringle et al. (2003), and using direct fitting with the model. These approaches are described in detail in the Appendix. Importantly, we find that the thermal diffusivity is a function of temperature, where κ_soil = 6.9 × 10^-7 m² s^-1 for -5°C ≤ T ≤ -25°C and rises from κ_soil = 6.9 × 10^-7 to 1.0 × 10^-6 m² s^-1 for T = -25°C to -45°C. The increase in κ_soil at low temperatures is probably due to frost cementation of the dry grains, as described in the Appendix, with a smaller effect due to the change in thermal conductivity and heat capacity of mineral grains and ice with temperature. This is directly seen by the strong correlation of thermal diffusivity with humidity increasing to ~100% for temperatures < -25°C and concurrently the thermal diffusivity increasing below this temperature. Using a representative soil density of 1630 kg m^-3 and a heat capacity of 800 J kg^-1, as per McKay et al. (1998), this gives a soil thermal conductivity (k_soil) of 0.9 W m^-1 K^-1 for -5°C ≤ T ≤ -25°C, and k_soil rises from 0.9 to 1.3 W m^-1 K^-1 for T = -25°C to -45°C.

For the ice-cemented ground, we use the values published in McKay et al. (1998) of κ_{icy_soil} = 1.0 × 10^-6 m² s^-1, using their representative density (2022 kg m^-3), heat capacity (1200 J kg^-1) and thermal conductivity (2.5 W m^-1 K^-1).

Latent and sensible heat fluxes

The latent and sensible heat fluxes require as input the roughness length scales for momentum, temperature and water vapour. These length scales in effect characterize the turbulent flow that is generated by the surface features at the site: rocks, boulders and vegetation. These values are accordingly specific to the site and can differ by orders of magnitude between sites (e.g. between a smooth glacier surface and a forested area).

The momentum roughness length scale is a measure of the aerodynamic roughness of the surface, while the temperature and water vapour roughness length scales represent the surface height below which the representative surface value is reached for the respective quantity. No measurements of these values, either by profile measurements or sublimation measurements (e.g. Denby & Snellen 2002, Wagnon et al. 2003), have been made for University Valley. We also did not find measurements for surfaces that are similar to the rough and rocky surface of the study site. Mölg & Hardy (2004), in their modelling of the glaciers on Kilimanjaro, used a value of 0.001 m for all three roughness length scales based on measurements at Canada Glacier, Antarctica (Lewis et al. 1999).

Both theoretical and field evidence suggests that the momentum, temperature and water vapour length scales should not be expected to be the same (e.g. Andreas 1987, Greuell & Smeets 2001), and generally the temperature and water vapour length scales are similar and one to three orders of magnitude smaller than the momentum roughness length scale (Mölg & Hardy 2004). In the case of University Valley, we expect the length scales to be larger than those measured at the glaciers due to the greater roughness of the surface.

The surface of University Valley is covered by material ranging from sand grains to boulders approaching 1 m in height (Fig. 1; also see fig. 3 in Heldmann et al. 2013) and patterned ground with polygons that range in size from 10 to 25 m in diameter (Mellon et al. 2014). Most of the published studies of turbulent heat and momentum transfer in the Antarctic are over ice or snow surfaces that are relatively smooth compared to University Valley.

In the context of aeolian transport of sand, Lancaster (2004) conducted a systematic field study of roughness length in a variety of Dry Valley surfaces ranging from flat sand to rocky moraines. For the smoothest surface he derived values of z_0m = 0.0009 m, while for the roughest surface he obtained z_0m = 0.0360 m. Liu et al. (2019) computed z_0m at Zhongshan Station at a rocky site (see fig. 1 in Liu et al. 2019), inferring a momentum (aerodynamic) roughness length of z_0m = 0.0036 m. They computed the thermal roughness length from the momentum roughness length using the approach of Andreas (1987), obtaining a value of z_0h = 0.00012 m and a ratio of z_0m/z_0h = 30.

For our modelling, we use the rockiest site in Lancaster (2004) to set z_0m = 0.0360 m. We follow Liu et al. (2019) in using a ratio of 30 between z_0m and z_0h, and we follow Andreas (1987) in setting the temperature and water vapour roughness length scales equal to each other: z_0h = z_0v = 0.0012 m.

We also tried to independently determine the roughness length scales at University Valley by fitting the data during relevant times. Specifically, after determining the surface albedo, emissivity and subsurface thermal diffusivity, we looked for low-insolation periods with strong winds and sudden changes in wind speeds. This increased the relative importance of the latent and sensible heat fluxes in the energy balance equation. Searching the plausible parameter space, this fitting supported the assumption that the temperature and water vapour roughness length scales are approximately equal and that the values of the roughness length scales selected above are approximately in the middle of the range of well-fitting values.

von Kármán constant

For the von Kármán constant, we use the traditionally accepted value of 0.4 (e.g. Högström 1985), but more recently Andreas et al. (2006) have suggested that a value of 0.387 may be more accurate. We proceed with using a value of K = 0.4, but note that the lower value will decrease both the latent and sensible heat fluxes by 6.8%. The resulting effects are explored in the sensitivity studies described below.

Longwave radiation

The longwave radiation parameterization that we used is only one of the possible parameterizations. In particular, significant work has been done on the change in longwave radiation parameterization depending on whether the sky conditions are clear or cloudy. In the case of our dataset, we can use apparent atmospheric optical depth to determine cloudy conditions during the sunny months; however, this is not possible for the long, dark winter, when longwave radiation is a more significant component of the radiation budget.

The nominal values used in the modelling are shown in Table I. These values are primarily determined using the methods described above. Where the more general fit suggests a range of values, the ones that preferentially show a good fit during the warm summer days are chosen, as it is during the warm summer days that the question of the availability of water is most important and the surface temperatures peak.

Results

Collected weather data

We recorded 3 complete years of data (2010–2013) at the weather station site (marked with a ‘W’ in Fig. 1). The averages and extrema of temperature and RH are shown in Table II. For each location, the frost point is calculated based on averaging the water vapour density computed from individual values of the temperature and RH (e.g. McKay 2009). Table III lists the average sunlight and wind speed, which are important parameters for comparison with other meteorological stations throughout the Dry Valleys. Figure 3 provides an overview of the entire dataset, showing the daily average air temperature, sunlight, RH and wind speed.

Fig. 3. Daily averages for a. air temperature, b. sunlight, c. air relative humidity and d. wind speed at 1.2 m above the surface for the University Valley weather station site. Numerical values are in the Supplemental Material.

Table II. Summary of temperature and humidity data (2010–2013).

^a The surface and frost point data could not be calculated for the 2010–2013 dataset directly because of missing data, as described earlier. The reported values are calculated by taking the continuous 358 days of surface temperature and humidity data that was collected using the Onset U23 Pro v2 temperature/relative humidity sensor from December 2009 to December 2010 and repeating the last 7 days again to make a complete 1 year (365 day) dataset. Additional surface temperature data from the Campbell T107 temperature probe are available; however, as these data do not provide full-year datasets, they are not reported in this table.

^b A maximum surface temperature of 12.3°C was recorded using the Campbell T107 temperature probe in the summer of 2012.

Table III. Sunlight and wind speed for each year (2010, 2011, 2012).

From the peak temperatures at 10 and 20 cm depths in Table II, we estimate that the depth of the active layer (the maximum depth of the 0°C isotherm) is ~15 cm. Given the variability in conditions across University Valley, this active layer depth is only applicable locally to the weather station site. This compares well to the value of 12.5 cm determined for a similar site (Linnaeus Terrace in Wright Valley; McKay et al. 1998). Considering all four summers of the dataset, the depth of the active zone varied by 4 cm and the mean was 13.2 ± 2.2 cm.

As can be seen in Table II, the frost point of the surface and the frost point of the ice table 42 cm below the surface are identical, providing direct confirmation that the surface conditions are the relevant boundary conditions for ground-ice stability and not atmospheric conditions measured at 1–2 m above ground height.

Table IV lists the summer monthly averages and extrema of temperature for those summer months for which the data are complete. Figure 4 shows the average wind speed as a function of wind direction and also shows the correlation between air temperature and RH with wind direction and the frequency of wind directions.

Fig. 4. Distributions of a. wind velocity and c. relative humidity (RH) with wind direction, as well as b. frequency of wind directions for the weather station site.

Table IV. Monthly summer air temperature values.

Model validation

We tested the model for fit against the continuous 358 days of surface temperature and humidity data from the Onset temperature/RH sensor from December 2009 to December 2010. Comparing the model and measured surface temperatures, the mean error is 1.3°C, the median error is 1.7°C and the mean absolute error is 2.1°C. By tuning the parameters noted in Table I within their realistic range, this error could be considerably reduced. For example, by reducing the momentum roughness length scale to 0.018 m (from 0.036 m) and setting the temperature and water vapour length scales to be 100 times smaller, the mean error reduces to 0.9°C, the median error reduces to 1.2°C and the mean absolute error reduces to 1.6°C. These possible corrections notwithstanding, we proceed with the values as given in Table I, as they still provide a reasonable fit for surface and subsurface temperatures for the whole year and a very good fit for the summer time period used for the further sensitivity studies, as described next.

We use a reference sunny and warm summer period to assess the fit of the model to the surface and subsurface temperature data (Fig. 5). As described later, this same time period is used to look for the effects of changes of parameters on the surface and subsurface temperatures. In this as in the other cases described, the model uses the atmospheric conditions (air temperature and humidity and wind velocity) and solar insolation as inputs to the energy balance, and the surface and subsurface temperatures are the calculated outputs. The reference period is 7–17 November 2010 (days 310–320). This period is unusual and particularly interesting as it includes multiple days of above-freezing surface temperatures, with a maximum surface temperature of 5°C. The model fits the data for this period quite well, as can be seen in Fig. 5, where the surface temperature has a mean error of -0.4°C, a median error of -0.12°C and a mean absolute error of 1.6°C. The peak temperatures are fit within 1°C. The largest errors are due to small offsets in the rise and fall rate during each day. As the temperature changes sharply when the sun rises or sets behind the valley walls, an offset of just half an hour in the response can produce errors of 5°C in the modelled vs measured temperature. The response time of the temperature sensors may also play a role in these cases. The overall close fit to the data allows for meaningful conclusions to be drawn about temperature changes as a result of varying environmental parameters.

Fig. 5. Reference warm period used for assessing the sensitivity of the environment to changing parameters; data collected at the weather station site. The model fits the data to within a few degrees Celsius, which allows for meaningful assessment of temperature changes when the environmental parameters are varied.

It is important to note, however, that this period does not have the warmest surface temperatures recorded during the 3 years of data collected from the site; the maximum recorded surface temperature was 12°C compared to 5°C in the reference warm period. From the available data, this is also seen in the 10 cm depth temperature measurement; for the period being analysed (up to day 341, as described later), the maximum temperature is -5.5°C, while the overall highest temperatures recorded at this depth during the summers of 2010–2011, 2011–2012 and 2012–2013 are +2.7°C, +2.8°C and +2.5°C, respectively. Thus, the results below should be interpreted specifically regarding the warm period being analysed, but they do provide a useful test of the ability of the model to match warm surface temperatures.

The model is initialized with a subsurface profile approximating this time of the year (confirmed by the subsurface temperature measurements), starting 50 days preceding the start of the reference warm period shown in Fig. 5 (day 260). This ensures that the model has equalized and the subsurface profile is expected to accurately represent the subsurface temperatures. The 50 day damping depth is 97 cm, which extends significantly past the 42 cm dry soil to the ice-cemented ground interface. The model is also run for 21.6 days after the end of the period of interest to allow for the thermal wave to propagate to the ice and for a more accurate maximum ice temperature to be obtained. For 21.6 days, the damping depth is 64 cm; however, the peak of the thermal wave is not expected to reach 42 cm in this timeframe. A longer timeframe is not simulated as this is the last data point that is available for the surface humidity dataset, and running the model further would require extrapolating the humidity data.

Sensitivity of surface and subsurface temperatures to environmental variations

The warm, sunny 10 day period described above (Fig. 5) is used to assess the sensitivity of the surface and subsurface conditions to changing environmental factors (Tables V & VI), with the goal of understanding both the effect of natural variability seen in the local University Valley environment (e.g. surface optical properties and textures) and the effect of variations from year to year and between neighbouring local valleys (e.g. wind profiles and air humidity). In the former case, the fitting parameters, such as albedo, emissivity and roughness length scales, are varied within the reasonable range and assessed uncertainty of the values, capturing the natural variability in the area (as described earlier for each parameter in the ‘Methods’ section). In the latter case, we vary the parameters within a plausible range of year-to-year and local inter-valley variability. The sensitivity studies also point to what parameters are especially important to determine with more certainty in future studies in order to more accurately model local conditions.

Table V. Modelling of the sunny, warm period in November 2010 and how changing local conditions and environmental parameters modifies the surface and subsurface temperatures. Unless otherwise noted, the environmental and subsurface parameters are those described in Table I. In all cases, z_0h = z_0v. Thaw depth is the depth to which the subsurface is warmed above freezing, regardless of water content.

Table VI. Effect of changing the depth of ice-cemented ground. Unless otherwise noted, the environmental and subsurface parameters are those described in Table I. In all cases, z_0h = z_0v.

Four quantitative metrics were used to compare the system response to parameter variations: maximum surface temperature, maximum ice temperature, the depth to which the subsurface is warmed above freezing regardless of water content (thaw depth) and the degree days above freezing at the surface (integral of the temperature-time curve for when the surface is above freezing). These quantitative metrics are evaluated as follows: maximum surface temperature and degree days above freezing for days 310–320 and maximum ice temperature and thaw depth during days 310–341.6 (Table V).

For this time period, the baseline model, using the best fit and as-measured environmental properties summarized in Table I, gives a maximum surface temperature of 5.8°C, a thaw depth of 4 cm and a maximum ice temperature (at 42 cm) of -18.9°C. These values are in close agreement with the measured values. The very shallow thaw depth demonstrates how the low thermal diffusivity of the warm, dry soil quickly dampens out the thermal wave. However, as noted earlier, this modelling run is more specific to the timeframe being analysed, and the thaw depth computed here is not to be confused with the depth of the active layer. The depth of the active layer is the maximum thaw depth achieved over the entire summer and, as discussed above, it averages 13.2 cm.

The variation in parameters that may be expected at the site has a moderate effect on the surface temperatures. Specifically, the variations in albedo, emissivity, thermal diffusivity of the soil, wind variability, roughness length scales and solar variability result in up to ~3.5°C warming in the surface temperature and ~1°C warming of the ice-cemented ground at 42 cm, keeping it well below the freezing point. The thaw depth varies by only a few centimetres. The warmest surface results from increasing the solar insolation by ~20%, which is at the upper limit of what may be plausible (i.e. effectively a significant decrease in cloud cover). The highest ice temperature occurs as a result of increasing the thermal diffusivity of the soil, which then allows more of the surface heat to be carried to depth.

The sensitivity to surface albedo is approximately ±1.3°C of the surface temperature for a corresponding change of approximately ±0.05 in albedo. The albedo value is expected to be most strongly affected by the surface material properties, grain size and presence of rocks, with a probably insignificant effect of whether pore-filling ice reaches the surface.

Emissivity changes of ±0.02 result in corresponding surface temperature changes of ±0.2°C.

Wind changes of ±25% result in corresponding surface temperature changes of approximately ±1.8°C, with higher winds making the ground colder. Higher winds result in a stronger coupling between the air and surface temperature, cooling the ground to the air temperature in summer and warming the ground in winter.

Solar flux changes of ±10% result in corresponding surface temperature changes of approximately ±1.9°C.

Changes in the roughness length scales have an effect that is comparable to the values discussed above. Thus, roughness length scales may be important in comparing results between sites with different surface characteristics.

There is an observed change in depth to ice-cemented ground down the length of University Valley. We noted that when the ice-cemented ground is exposed at the surface there was no melting observed during the three summer field visits. In using the model to simulate different ice depths, with all other parameters being as described in Table I, we find a similar result (Table VI): ice-cemented ground at the surface does not melt in the summer. There is some ambiguity in these results as the observation period during a field visit is brief, the modelling does not use the warmest period recorded in the 3 years of data and the warmest temperatures recorded at 10 cm at the weather station site were ~2.5°C, although the data from an ice-cemented ground depth of 42 cm cannot be directly related to temperature profiles when the ice reaches the surface.

Changing the depth of the ice in the model changes the computed subsurface temperature profile dramatically, as expected. Modelling ground ice all the way to the surface gives a maximum surface/ice temperature of approximately -2.5°C. Our measured temperature at the shallowest ice site (depth to ice of 8 cm) gave a maximum ice temperature of -16°C during this same timeframe and a maximum ice temperature for the year of -4°C. Modelling an ice depth of 8 cm at the weather station site gives a maximum ice temperature of -7°C during the modelled timeframe. These results are overall consistent, and some discrepancy is expected as the narrowness of University Valley results in rapidly varying shadowing and solar insolation throughout the valley floor; other environmental factors are expected to vary as well, changing the thermal environment.

It should be noted that to model intermediate depths of ice we must specify a surface humidity value. This is done as described earlier after Eq. (7). As there is uncertainty in how to treat the surface humidity for shallow ice depths, the latent heat component of the energy balance is more uncertain.

Shallower depths to ice-cemented ground increase the thermal diffusivity, which enables heat from the surface to be more easily conducted into the subsurface. In addition, the higher thermal heat capacity of ice-cemented ground results in the same heat input resulting in a smaller bulk temperature increase. Thus, even when ice is present close to what is normally expected to be a warm surface from solar heating, the presence of the ice modifies the subsurface properties in such a way that the ice is still kept below freezing and no melting occurs. This general observation is expected to hold during average years, but exceptional weather excursions can still result in conditions that may melt subsurface ice, particularly for shallower ice. This is discussed further in the ‘Discussion’ section.

Discussion

Climate comparison

To place the climate of University Valley and the conclusions we draw from our data and modelling into the broader context of climate in the Dry Valleys, we compare the overall climate of University Valley to the extensive datasets that exist for other locations in the Dry Valleys. These include early datasets (Keys 1980, Clow et al. 1988), the Long-Term Ecological Research (LTER) network (Doran et al. 2002), the stations of the Soil Climate Monitoring Project (Seybold et al. 2009) and the stations established at high elevations to study the cryptoendolithic colonization in the Beacon Sandstone (Friedmann et al. 1994, McKay 2015).

The average air temperature at University Valley is -23.4°C, comparable to the values for similar elevations in Wright Valley (Friedmann et al. 1994), and this is only a few degrees colder than the average at some of the low-altitude lake sites near the coast. For example, Doran et al. (2002) report an average temperature at Lake Fryxell of -20.2°C. Comparing University Valley to Lake Fryxell, the long, cold winters at University Valley are not much colder than at Lake Fryxell, but the summer months are much colder. This is due mainly to the effect of elevation on summer air temperature. Lake Fryxell has a peak summer air temperature of 9.2°C and an average of 25.5 degree days above freezing (integrated time during which the air temperature is above freezing). At University Valley, the peak summer air temperature effectively never gets above zero; there are no summer degree days above freezing. Over the 3 year data period, one temperature point (an average over 30 min) was +0.4°C with a measurement uncertainty, as discussed in the ‘Methods’ section, of ±0.25°C. The value of degree days above freezing has been shown to be an important indicator of meltwater production from glaciers in the Dry Valleys (Wharton et al. 1992, Doran et al. 2008). Wharton et al. (1992) showed that increases in the lake level of Lake Hoare tracked the degree days above freezing over the summer. At University Valley, the lack of degree days above freezing also correlates with the observed lack of melt.

Many features of the climate in the Dry Valleys are determined by the wind patterns. The winds vary between mild coastal winds that bring cold, moist air up the valleys and strong downslope winds that create warm and dry conditions (Keys 1980, Clow et al. 1988, Doran et al. 2002, Nylen et al. 2004, Fountain et al. 2010). Figure 4a shows the average wind speed as a function of direction at University Valley. The distribution includes strong winds coming down the valley (90°) and milder winds moving up the valley (270°). Figure 4c complements the data in Fig. 4a, showing the variation of RH with wind direction. In the lower valleys, the winds show a strong correlation between humidity and direction (e.g. Clow et al. 1988). This is less evident in University Valley. For example, the strong downslope winds at University Valley have relatively high humidity compared to the corresponding downslope winds at Lake Hoare and Lake Fryxell. This is because as they descend from the Polar Plateau the winds follow an approximately adiabatic trajectory through the valleys with constant absolute water content. Thus, the further down the valleys, the warmer these winds become and the lower the RH. At the high elevation of University Valley, the effect is smaller.

Doran et al. (2002) suggested a relationship between summer temperatures at stations in the Dry Valleys of the form:

(11)$$\Delta T = 0.09\Delta d-10\Delta z$$

where ΔT is the difference in mean annual air temperature between any two stations, Δd is the difference in distance from the coast (km) and Δz is the difference in elevation (km). Doran et al. (2002) showed that the temperature difference between sites correlated with the cold coastal winds that pick up more heat as they move inland. Thus, the inland stations are cooled less than the coastal stations and thus the temperature was correlated with the shortest direct distance from the coast. The second term on the right-hand side of Eq. (11) is the reduction of air temperature moving upward with elevation by the dry adiabatic lapse rate (rounded to 10°C km^-1). McKay (2015) showed that Eq. (11) predicts the temperature at high-altitude stations in Wright Valley, including stations as high as 2178 m, and 80 km distant from the coast.

Table VII compares the summer monthly average air temperatures at University Valley to two LTER stations: Beacon Valley and Lake Hoare. The Beacon Valley station was located at 77°49.681′S, 160°38.422′E, elevation 1176 m. This station is 500 m lower and 4 km distant from our University Valley station, which we take as the values of Δz and Δd, respectively. The Lake Hoare station was the first of what is now the LTER stations (Clow et al. 1988) and is located at 77°37.523′S, 162°53.999′E, elevation 73 m, and at a distance to the coast of 15 km (Doran et al. 2002). The summer monthly averages generally follow the Doran et al. (2002) rule given the spread in the data.

Table VII. Comparing monthly average air temperatures to the predictions of Doran et al. (2002).

LTER = Long-Term Ecological Research.

From these comparisons, we can conclude that the meteorological conditions at University Valley follow the expected trend of elevation vs temperature change as seen in high-elevation stations in other valleys and given the published correlations between stations at low and high elevations.

The correlation of temperatures between University Valley and other stations can be used to estimate maximum air temperatures over an extended period. The measured maximum air temperature at Lake Hoare from 1985 to 2018 was 10°C and the maximum air temperature at the LTER Beacon Valley station from 2000 to 2012 was 2.8°C (Obryk et al. 2020). Using the observed and predicted temperature differences between these two stations and University Valley (listed in Table VII), the maximum air temperature at University Valley is < 0°C based on the Beacon Valley data and +0.1°C (observed) and +1.4°C (predicted) based on the Lake Hoare data. This analysis suggests that over the instrumental record from 1985 to 2018 University Valley did not have maximum summer air temperatures > 0°C, which is within the uncertainty of the measurements and extrapolations.

Ice stability

Ice-cemented ground under dry permafrost occurs when the depth to ground ice exceeds the depth of the active zone: the depth at which the peak temperatures reach 0°C. As discussed before (e.g. McKay 2009, Fisher et al. 2016, McKay et al. 2019), the depth to ground ice is determined by the frost point temperature of the surface. As illustrated in fig. 1 of McKay (2009), the frost point temperature of ice-cemented ground is maximal at the surface and decreases with depth until it becomes equal to the mean annual temperature at depths below ~2–3 damping depths. If the frost point temperature of the surface is lower than the mean annual temperature of the surface, then ice is unstable at all depths. If the frost point of the surface is equal to or above the frost point corresponding to the depth of the active zone, then the top of the ground ice is in the active region and melts in the summer. For our weather station site, the frost point of the surface is -22.5°C and the mean annual temperature of the subsurface is -24.1°C. A reduction in the absolute water content of the surface by a factor of 0.86, as measured by the surface humidity sensor, would lower the frost point of the surface from -22.5°C to -24.1°C and ice would not be stable at any depth. Similarly, an increase in the absolute water content of the surface by a factor of 1.12 would raise the frost point temperature of the ice table from -22.5°C to -21.3°C (the frost point temperature at the bottom of the active zone at ~15 cm if the RH at that level is set to 100% all year; see the last column in Table II) and the top of the ice table would just reach the melting point in the summer.

In University Valley, depth to ground ice varies significantly along the valley (Marinova et al. 2013), which we suggest is primarily due to reduced surface moisture. At the upper end of the valley, ice-cemented ground is present at the surface due to lower surface temperatures resulting from reduced sunlight in the narrow semicircle of the valley headwalls; this area may also receive additional snow cover from snow blown in from the Polar Plateau. Throughout the valley floor, the depth to ground ice in and across polygons varies considerably, which we suggest is due to surface moisture variations associated with snow that preferentially persists in the cracks between polygons, as is evident in fig. 2 of Mellon et al. (2014).

Ice melting and habitability insights from the modelling

The previously discussed sunny days are also used to explore a related question of interest: what conditions are required for ice melting to occur or, conversely, what is the sensitivity of water availability to changing environmental conditions? The same period of sunny, warm days is used (days 310–320 in 2010), as was already discussed in the ‘Results’ section and shown in Fig. 5. The effects of varying the environmental conditions and whether they resulted in ice melt are summarized in Table VIII. It should be noted again that the 10 days used for the modelling here represent one of the warmest periods for which we have a complete dataset from the 3 years of collected data; however, it is not the single warmest period. This period was used due to its data completeness and suitability for the intended modelling.

Table VIII. Climate sensitivity study. Unless otherwise noted, the environmental and subsurface parameters are those described in Table I. In all cases, z_0h = z_0v. Elevation changes use a dry adiabatic lapse rate of 9.8°C km^-1 (McKay 2015) with a reference elevation of 1677 m (University Valley). Parameter changes that are not realistic for the current environment are intended to be illustrative of the change required for ice-cemented ground melting to occur.

The modelling showed that the expected magnitude of local variation in environmental conditions (e.g. ice depth, wind speed, etc.) does not seem to be sufficient to allow shallow (surface) ice to melt in University Valley (Tables V & VI). This is consistent with the presence of soluble salts and clay-sized particles preferentially concentrated near the surface (Tamppari et al. 2012) and isotope data that indicate that the uppermost ice layers formed by condensation-diffusion of water vapour (Lacelle et al. 2013). The consistent modelling and observations speak to the extremely dry nature of this environment. For liquid water to be present in the soil, large temperature excursions are required. Such warm events are expected to be relatively rare, although we do not have the data to quantify their frequency. The conditions required to melt the ice are discussed below.

It should be noted that the model does not take into account the heat of fusion required for the ice to melt, and thus the temperatures of the ice-cemented ground that are slightly above freezing probably represent thin water films or an ice-water mixture rather than full melting and further temperature increase.

The air temperature has a strong effect on the subsurface temperature through its sensible heat flux contribution: the air and surface temperatures are coupled by the wind and mixing length scales. The air temperature, which is cooler than the sun-warmed surface, effectively cools the ground when there are strong winds and/or short mixing length scales. Conversely, slow winds and larger mixing length scales mean that the air and surface temperature are effectively decoupled and the surface can warm up significantly to above the air temperature. Specifically, modelling an increase in air temperature by 5°C and 10°C give increases in thaw depth of 4 and 9 cm, respectively, for ice depths of 42 cm. In the case where the ice-cemented ground reaches to the surface, an air temperature increase of 5°C allows the ice to reach melting. This suggests that during the warmest recorded periods, when surface temperatures at the weather station site were up to 7°C warmer, in areas where the ice reaches the surface a melting temperature may have been reached. We do not have comparison data for this time period at the shallow ice site to make a further assessment. Interestingly, the model results show that if the ice-cemented ground is even 5 cm below the surface, an air temperature increase of even 10°C means a thaw depth of 3 cm and a maximum ice temperature of -1.1°C, suggesting that the ice is not melting. For the warmest times in the dataset, when above-freezing temperatures are reached at the 10 cm depth in dry permafrost (d_ice = 42 cm), it is unclear whether above-freezing temperatures would reach the 10 cm depth if the ice was at 10 cm.

Wind plays an important role in the maximum temperature reached by the surface, as both the sensible and latent heat fluxes scale with the wind velocity. Through the sensible heat flux, in summer, stronger winds cool the surface to the air temperature, while weaker winds allow the surface to warm significantly in sunny conditions. In winter, stronger winds warm the surface. During the modelled reference warm period (Fig. 5), the wind speed varies from ~0 to 7.7 m s^-1. Dividing all wind speeds by a factor of 2 changes the maximum surface temperature from 5.8°C to 10.8°C for the nominal ice depth. This effect on the surface temperature is similar to increasing the air temperature by 5°C. With an ice depth of 5 cm, the maximum surface temperature using half the wind speed is 3.5°C and the maximum ice temperature is -4.6°C. With ice-cemented ground to the surface and half the wind speed, the surface temperature reaches just above freezing, at ~0.5°C.

Changing the mixing length scales within what is considered to be the reasonable range for the location does not result in the surface reaching freezing when the ice-cemented ground reaches the surface.

To evaluate the variability in conditions in the broader Dry Valleys area, we examine the expected surface and subsurface ice temperatures at a valley at 500 m elevation and its summer air temperature that correspondingly increases by 11.5°C (using a dry adiabatic lapse rate of 9.8°C km^-1); this elevation is equivalent to Pearse Valley, where we also have observational data that can be compared to these results. The modelling results, which use the University Valley climate conditions except for the increase in air temperature by 11.5°C (corresponding to the elevation change), show that the thaw depth reaches a depth of 14 cm when the ice-cemented ground depth is 42 cm. Bringing the ice-cemented ground depth to the surface clearly results in melting (notional, but overestimated, ground/ice temperature of 4.9°C). The maximum depth at which melting occurs is ~5 cm, although this depth is affected by at least a few centimetres based on surface humidity and soil thermal diffusivity, as well as the other surface properties, which may be different from those in University Valley. These results are consistent with observations at Pearse Valley, which show yearly melting of near-surface ice (Heldmann et al. 2012).

Looking at the importance of the subsurface thermal diffusivity on the thermal wave, we note that the model results for the reference warm period give a maximum surface temperature of ~5.8°C, while when the ice-cemented ground depth is set to 10, 5 and 0 cm, which overall increase the effective thermal diffusivity of the subsurface column, the maximum surface temperatures decrease to 1.9°C, -0.4°C and -2.5°C, respectively. The modelling shows that increasing the ice table depth beyond 10 cm has little effect on the peak soil temperature. This is not surprising, as the daily damping depth in dry permafrost is ~14 cm. Essentially, dry permafrost dampens the thermal wave quickly, inhibiting heat from the warm surface from propagating to any significant depth, while ice-cemented ground conducts the surface heat very efficiently, not allowing the surface and near-surface ice to warm as much. These two effects make it difficult for ice melting to occur in these extremely cold environmental conditions.

These sensitivity studies show that significant excursions in temperature and/or other environmental conditions are required for melting to occur at University Valley. This has important implications for both geological processes, such as weather and salt transport in the soil, and for habitability. In the latter case, the work by Goordial et al. (2016) has shown that the lack of available liquid water due to melt is severely constraining microbial life, as the culturable and total microbial biomass in University Valley soils is extremely low and microbial activity was undetectable in laboratory assays. This is in stark contrast with reports from the nearby lower-elevation Dry Valleys, where more abundant and metabolically active soil microbial populations are found (Cary et al. 2010).

Dry permafrost on Mars

Ice-cemented ground under dry permafrost is widespread on Mars, making University Valley an analogue to Martian conditions. The distribution of subsurface ice on Mars was predicted by models (e.g. Mellon & Jakosky 1993). These predictions were largely confirmed by the orbital data from the Odyssey mission, which indicated widespread subsurface ground ice in the polar regions of both hemispheres on Mars (Feldmann et al. 2004). Direct observations by the Phoenix mission, which landed at 68.2°N on Mars, found an ice table below the dry surface with depths ranging from 1.3 to 11 cm and with an average depth of 4.6 cm (Mellon et al. 2009).

The similarity and comparison between University Valley and Mars extend to the availability of water and hence habitability. The potential for cold, ice-rich regions on Mars to support biological growth has prompted the delineation of Special Regions on Mars. Special Regions are currently defined (Rummel et al. 2014) as follows: ‘Special Regions on Mars continue to be best determined by locations where both of the parameters (without margins added) of temperature (above 255 K) and water activity (a_w > 0.6) are attained’. Here, we compare the temperature and water activity of the ice-cemented ground in University Valley to the limits of Special Regions (temperature > -18°C, a_w > 0.6). Figure 6 shows the cumulative time spent above a given temperature in the permafrost at 42 cm depth at the University Valley weather station site. Also shown is the water activity of ice computed as the ratio of the vapour pressure of ice divided by the vapour pressure of water (Murphy & Koop 2005). This plot shows that the water activity at the ice table is always above the limit set for Special Regions on Mars and the temperature exceeds the limit for considerable periods. For example, 129 h were recorded with temperatures > -10°C and 2560 h (~30% of the year) were recorded with temperatures > -18°C. In summary, the temperature and water activity values at the ice-cemented ground are within the limits set for Special Regions on Mars. By contrast, numerical models based on Martian environmental data predict that near-surface ice-rich regions on Mars never cross the limit set for Special Regions (Rummel et al. 2014). This makes University Valley a particularly relevant analogue for the study of Martian habitable regions. While the water activity and temperature conditions in University Valley soils meet the definition of Mars Special Regions for a significant fraction of the year, the work by Goordial et al. (2016) shows that the soils do not appear to support observable microbial activity and growth, and microbial activity was undetectable in laboratory assays under ambient conditions using techniques that have worked in every other permafrost environment assessed to date.

Fig. 6. Temperature and water activity at the ice table (42 cm deep) at the University Valley weather station site. The black curve is the cumulative time spent above a specified temperature in hours per year. The blue curve is the water activity set by the temperature of ice computed using the formulae in Murphy & Koop (2005).

Conclusions

We have recorded the meteorological parameters at University Valley, a high-elevation valley in the Dry Valleys of Antarctica, over 3 years. We have developed a detailed energy balance model of the surface overlying dry permafrost and ice-cemented ground in this valley. Based on the comparison of the observations to other locations in the Dry Valleys and from the comparison of the model to the observational data, we draw the following conclusions.

University Valley's climate follows the trends reported from other Dry Valley stations that are mostly at lower elevations. In particular, we find that the summer temperatures at University Valley scale with elevation using the dry adiabatic lapse rate and distance from the coast. The mean annual air temperature is -23.4°C with virtually no summer degree days above freezing.

Values were obtained for surface and subsurface environmental parameters (e.g. albedo, emissivity, thermal diffusivity, roughness length scales) by fitting data periods where the effect of a parameter was important or dominated the energy flux. Using this method, we can fit a sunny and warm 10 day period to within ~1°C and the full dataset to a mean error of 1.3°C, suggesting that these fitting methods can also be used elsewhere to better characterize surface and subsurface properties.

Thermal diffusivity is fitted using multiple methods for comparison and is found to increase as the temperature decreases, making it time-variable. This correlates with the increase in humidity in the subsurface, and we interpret this as ice condensing on the soil grains and thereby significantly increasing the conductivity of the soil.

The depth to ice-cemented ground plays an important role in modulating the temperature of the surface. With a thick dry soil layer (> 10 cm), the surface can warm up in summer and reach temperatures above freezing (the maximum surface temperature was 7.5°C as measured by the Onset sensor), as the dry soil effectively acts as an insulator. When ice-cemented ground is closer to the surface, including reaching the surface, our modelling shows that its higher thermal diffusivity effectively wicks heat away and limits the warming of the surface. This effect means that when ice reaches the surface it does not experience above-freezing temperatures and does not melt under these conditions. This is supported by observations. Extreme weather events (not seen in our dataset) are required for the ice to melt.

The surface and subsurface temperatures are notably increased by slower winds, lower thermal diffusivity of the subsurface, smoother (less rough) surfaces, larger ratios between the momentum and temperature/water vapour length scales and higher air temperatures (equivalent to lower elevations in the Dry Valleys). However, varying all parameters independently within the plausible range did not result in melting occurring for the simulated time period and the observed ice conditions in University Valley.

The co-occurrence of dry permafrost and ice-cemented ground in University Valley depends on the difference between the frost point of the surface and the frost point at depth in the ground. We confirm previous suggestions that the frost point of the ice table is numerically equal to the frost point of the surface and considerably warmer than the frost point of the atmosphere. We find that a reduction in the absolute water content of the surface by a factor of 0.86 would lower the frost point temperature of the surface and ice would not be stable at any depth. Conversely, an increase in the absolute water content of the surface by a factor of 1.12 would raise the frost point temperature enough that the ice table would be at the depth of the active zone and there would be no dry permafrost. We conclude that this dependence on surface moisture is the main factor in the variation in the average depth to ground ice along the valley floor, and this is also the main factor impacting the variation in the depth to ground ice seen in and around polygons.

The temperature and water activity at the surface of the ground ice are often within the limits that define Special Regions on Mars. University Valley and similar sites in Antarctica thus provide operational analogues on Earth of Mars Special Regions.

Future work that collects extensive environmental information for similar environments will be very beneficial in further tuning these models and allowing for a more thorough investigation of these unusual environments using a combined in situ and modelling approach. In such studies, care should be taken to install subsurface temperature sensors in such a way as to be able to obtain the subsurface diffusivity properties based on the methods described in the Appendix. Additional work to measure the mixing length scales in these environments will significantly reduce the uncertainty in modelling results. Observing the site during the warmest time of the year, January, would provide additional insights on localized melt processes.