2.1 Constructing a fine-resolution regional snow data set

In a previous study (Yi et al., 2018) we used coarse-resolution (0.5^∘) snow depth data from the MERRA-2 global reanalysis as input for the permafrost model over Alaska. In the prior study, we first interpolated the MERRA-2 data over a finer 1 km spatial grid using an inverse-distance weighting scheme and then used the MODIS 500 m SCE data to identify snow-free pixels within each coarser MERRA-2 grid and adjust the 1 km snow depth estimates accordingly. However, a more sophisticated downscaling scheme was needed to better account for the influence of local topography on the 1 km snow distribution pattern. This information can be derived from the MODIS SCE data; however, persistent cloud cover and patchy snow conditions constrain the ability of the MODIS SCE data to capture snow cover variability, especially during the transitional season. To overcome this constraint, in this study we developed an elevation-based spatial filtering algorithm to predict the snow occurrence for MODIS cloud-contaminated pixels; we then used the cloud-free MODIS SCE data to downscale the MERRA-2 snow depth data (Fig. 1). For each snow-covered 1 km pixel indicated by the MODIS data, we estimated the snow depth based on the snow depth of neighboring MERRA-2 0.5^∘ grid cells, with weights predicted using a similar spatial filter.

2.2 The remote-sensing-driven permafrost soil process model

The newly developed snow depth data were used with other satellite remote-sensing data sets as primary input to an established permafrost soil model (Yi et al., 2018) for simulating soil freeze onset and zero-curtain period across Arctic Alaska. The model was developed based on a detailed permafrost hydrology model (Rawlins et al., 2013; Yi et al., 2015) but has a flexible structure designed to exploit remote-sensing observations as key model drivers and for model parameterization. The remote-sensing-based permafrost model, as described in Yi et al. (2018), uses a numerical approach to simulate soil F–T processes and the temperature profile down to 60 m below the surface using 23 soil layers, with increasing layer thickness at depth (soil nodes from 0 to 1 m: 0.01, 0.03, 0.08, 0.13, 0.23 ,0.33, 0.45, 0.55, 0.70, 1.05 m). Up to five snow layers are used to account for the effects of seasonal snow cover evolution on snow density and thermal properties. Both snow heat capacity and thermal conductivity vary with snow density and are estimated using empirical methods (Calonne et al., 2011). The model also accounts for the effects of organic soils and soil-water phase change on the soil F–T process as described below.

The soil model simulates snow and ground thermal dynamics by solving a 1-D heat transfer equation with phase change (Nicolsky et al., 2007; Rawlins et al., 2013):

where T(zt) is the temperature (^∘C) at a specific soil depth (z) and time step (t), L is the volumetric latent heat of fusion of water (J m⁻³), ζ is the volumetric water content (m³ m⁻³), and θ is the unfrozen liquid water fraction (range: 0–1). C and λ are the volumetric heat capacity ($\mathrm{J}{\mathrm{m}}^{-\mathrm{3}}{\mathrm{K}}^{-\mathrm{1}}$) and thermal conductivity ($\mathrm{W}{\mathrm{m}}^{-\mathrm{1}}{\mathrm{K}}^{-\mathrm{1}}$) of soil respectively, varying with soil moisture, F–T state and depth. The upper boundary condition is set as the surface temperature (i.e., LST) at the snow–ground surface (z_s), while a heat flux characterizing the geothermal gradient is applied at the lower boundary (z_b). The soil heat capacity is a function of heat capacities for soil solid and liquid water, and ice components (Farouki, 1981):

where C_m, C_o, C_w and C_i are the heat capacity of mineral and organic soil solid, liquid water and ice respectively, weighted by their volumetric function. f is the soil organic fraction; θ_sat, θ_w and θ_i are the respective saturated water content, volumetric liquid water and ice content, where $({\mathit{\theta }}_{\mathrm{w}}+{\mathit{\theta }}_{\mathrm{i}})\le {\mathit{\theta }}_{\mathrm{sat}}$. The thermal conductivity λ ($\mathrm{W}{\mathrm{m}}^{-\mathrm{1}}{\mathrm{K}}^{-\mathrm{1}}$) was estimated as a normalized thermal conductivity of the dry (λ_dry) and saturated (λ_sat) soil thermal conductivity weighted by soil saturation:

where the Kersten number (Ke) is a function of the soil saturation degree, generally using a logarithm form for unfrozen soils or linear form for frozen soils (Lawrence and Slater, 2008). λ_dry is estimated from the soil bulk density; λ_sat is estimated as a geometric mean of the thermal conductivity of different soil components (Farouki, 1981), which can vary several-fold from pure organic soil (∼0.5 $\mathrm{W}{\mathrm{m}}^{-\mathrm{1}}{\mathrm{K}}^{-\mathrm{1}}$) to mineral soils (1.5–3 $\mathrm{W}{\mathrm{m}}^{-\mathrm{1}}{\mathrm{K}}^{-\mathrm{1}}$):

where λ_m, λ_o, λ_w and λ_i are the thermal conductivity of mineral and organic soil solid, liquid water and ice, respectively; and θ_sat,w and θ_sat,i are the respective unfrozen liquid water and ice fractions under saturated conditions.

The unfrozen liquid water fraction (θ) is estimated empirically as follows:

Soil water usually freezes at sub-zero temperatures depending on solute concentration and other factors, and the constant T_∗ is used to represent this freezing-point depression, with values generally above −1 ^∘C (Banin and Anderson, 1974; Woo, 2012). b is a dimensionless parameter determined by fitting the unfrozen water curve (Romanovsky and Osterkamp, 2000; Schaefer and Jafarov, 2016). A significant amount of liquid water can exist even when the soil temperature is considerably lower than T_∗, characterized by different values of b. Fine-grained soils that can have a larger amount of liquid water below freezing are generally associated with smaller b values (Woo, 2012).

Figure 2(a) The study area (Arctic Alaska, >66.55^∘ N) used for the soil model simulations and the locations of in situ sites used for model calibration and validation; (b) an Alaska permafrost probability map (Pastick et al., 2015) is also shown, indicating higher permafrost occurrence in Arctic Alaska. The in situ sites include the Prudhoe Meadow SoilSCAPE site, GTN-P soil temperature sites and CALM ALT sites. Airborne P-band radar data (denoted by green lines) were obtained in late August and early October in 2014 and 2015. Another radar flight line was collected over the Barrow area (not shown) but was not included in this analysis due to a large percent of surface water in this area.

Download

2.3 Model driver data sets and in situ data

Model simulations were conducted in the Arctic Alaskan domain (>66.55^∘ N, Fig. 2), encompassing an area of ∼ 400 000 km² and spanning a 16-year period (2001–2016). Primary model drivers include MODIS 8-day composite 1 km LST (MOD11A2; Wan and Hulley, 2015) and 500 m SCE records (MOD10A2; Hall and Riggs, 2016), SMAP 9 km NatureRun (Version 4) and level 4 daily surface (≤5 cm depth) and root zone (0–1 m depth) soil moisture (L4SM, Reichle et al., 2017), and daily snow depth and snow density from MERRA-2 global reanalysis data (Gelaro et al., 2017). The MODIS LST and SMAP L4SM products were used to define model boundary conditions and soil thermal properties. The soil process model was run at 1 km resolution and 8-day time step consistent with the MODIS LST and SCE inputs. All model input data sets were reprojected to a consistent 1 km Albers projection. The model snow depth inputs were derived by combining the MODIS SCE and MERRA-2 snow depth records as discussed in Sect. 2.1. Compared with snow depth, snow density shows much smaller spatial and temporal variability (Sturm et al., 2010); therefore, we used the 1 km snow density data generated using a simple spatial interpolation scheme as described in Yi et al. (2018). Other ancillary inputs to the soil model included the 30 m National Land Cover Database (NCLD) 2011 (Jin et al., 2013), 2 arcsec (∼60 m) DEM for Alaska (U.S. Geological Survey, 2017), 50 m SOC estimates for Alaska (to 1 m depth; Mishra et al., 2016), and the global 9 km mineral soil texture data developed for the SMAP L4SM algorithm (De Lannoy et al., 2014). The dominant NLCD land cover type within each 1 km pixel was used to define the modeling domain, with open water and perennial ice and snow areas excluded from the model simulations. The soil texture and SOC data were used to define model soil properties including thermal conductivities and heat capacities. The SOC inventory data was distributed through the top 10 model soil layers (≤1.05 m depth) following an exponentially decreasing curve (Hossain et al., 2015) to calculate the SOC fraction and adjust the soil physical properties of each soil layer based on the weighted mineral and organic soil components. More details on the data processing can be found in Yi et al. (2018).

Three in situ data sets were used for model calibration and validation (Fig. 2), including half-hourly soil dielectric constant (ε) and temperature profile measurements from a Soil moisture Sensing Controller and oPtimal Estimator (SoilSCAPE) site (Moghaddam et al., 2010); daily soil temperature profile measurements from Global Terrestrial Network for Permafrost (GTN-P) sites (Biskaborn et al., 2015), and ALT measurements from regional Circumpolar Active Layer Monitoring (CALM) network sites (Brown et al., 2000). The SoilSCAPE soil temperature and ε measurements were obtained from four different depths (0.05, 0.15, 0.35, 0.56 m) and four different nodes of a wireless sensor network deployed near Prudhoe Bay, Alaska (70^∘13^′47${}^{\prime \prime }$ N, 148^∘25^′19${}^{\prime \prime }$ W) in the summer of 2016. ε was measured using a METER TEROS 12 soil moisture sensor operating at 70 MHz. For the GTN-P in situ measurements, we only selected the sites where shallow ground-temperature measurements (generally down to 1 m depth) were available for at least 2 consecutive years. Most GTN-P sites meeting these criteria in Arctic Alaska are located along the Dalton Highway (Table S1 in the Supplement). In addition, daily snow depth measurements using an ultrasonic sensor were available at SNOTEL (SNOwpack TELemetry) sites across Alaska (Schaefer and Paetzold, 2000; http://www.wcc.nrcs.usda.gov, last access: 1 September 2018, Fig. S1 in the Supplement) and used to validate the 1 km snow depth product (Sect. 2.1).

Figure 3In situ measurements of soil temperature (a) and dielectric constant (b) for a single SoilSCAPE sensor node (S6) at the Prudhoe Meadow site in northern Alaska (http://soilscape.usc.edu, last access: 1 May 2018).

Download

2.4 Model parameterization

Soil dielectric properties are strongly correlated with soil moisture, texture and F–T state (Mironov et al., 2010); ε can capture the soil-freezing process well due to large ε differences between liquid water and ice (Dobson et al., 1985), especially during the zero-curtain period, when soil temperatures hover around 0 ^∘C and are a relatively poor indicator of F–T conditions (Fig. 3). The SoilSCAPE measurements were used to calibrate the model unfrozen water content curve (Eq. 11) assuming a linear relationship between ε and liquid water content (Mironov et al., 2010; Park et al., 2017). However, the slope of this linear relationship may change during the freezing period due to different dielectric properties of free and bound water, and ice (Mironov et al., 2010). The ε measurements were also used to determine the timing of complete soil freeze-up at each soil depth. The timing of soil freeze-up was defined when the observed ε drops below a critical level:

where ε_max and ε_min are the maximum and minimum dielectric constants for each soil layer, and the threshold δ is an empirical parameter. The zero-curtain period at different depths was then calculated as the difference between soil freeze-up and land-surface freeze onset. The land-surface freeze onset was determined from MODIS 1 km LST records extracted at the SoilSCAPE site and defined as the date at which the mean LST during three consecutive 8-day periods dropped below 0 ^∘C.

The model-simulated soil temperature profile is very sensitive to the soil thermal conductivity, largely determined by soil texture (organic or mineral soils) and soil saturation (Eqs. 9–10), which are also two major factors affecting ε. Therefore, we used the in situ ε data at the SoilSCAPE site to guide the parameterization of model soil thermal properties. Since the soil is mostly saturated at this site, much larger ε values in the top two layers during the thaw season (Fig. 3) should be related to organic-rich soils with large soil porosity (thus high volumetric soil moisture). We defined the top five model soil layers (0–0.23 m) as organic soils and adjusted the model soil thermal properties accordingly. The unfrozen soil thermal conductivity within the organic layer was assumed to gradually increase with depth from ∼0.5 $\mathrm{W}{\mathrm{m}}^{-\mathrm{1}}{\mathrm{K}}^{-\mathrm{1}}$ in the surface organic soil layer to ∼1.2 $\mathrm{W}{\mathrm{m}}^{-\mathrm{1}}{\mathrm{K}}^{-\mathrm{1}}$ at 0.33 m depth for mineral soils, accounting for increases in the soil bulk density (Letts et al., 2000). Similarly, soil porosity was assumed to gradually decease from 0.8 at the surface to ∼0.4 in the deeper mineral soil layers. The soil thermal conductivity for frozen conditions can then be determined from Eq. (10). Using this soil thermal conductivity profile, the model-simulated temperatures agree well with the in situ observations (R>0.97, RMSE<2.24 ^∘C for all measured soil depths). We then tested different soil dielectric thresholds (δ) ranging from 0 % to 50 % and selected the threshold that produced minimum bias and RMSE between the zero-curtain period determined using in situ ε measurements and model-simulated unfrozen water content. Using this trial and error method, an optimal threshold of 15 % was selected, which produced a mean RMSE of 10.3 days in the simulated zero curtain from 0.15–0.56 m soil depth in 2016 and 2017 (Fig. 4a).

Figure 4Soil-freezing characteristics at the Prudhoe Meadow SoilSCAPE site: (a) comparisons of modeled and observed zero-curtain period at two sensor nodes (S5 and S6) using a cutoff threshold of 0.15 (Eq. 12) for both in situ soil dielectric constant and model-simulated unfrozen water content to determine soil freeze onset; (b) changes in in situ soil dielectric constant during soil freezing at the S6 node. The soil dielectric constant was normalized using the maximum and minimum values of soil dielectric constant during the observation period.

Download

The resulting threshold δ was then used to determine the critical threshold of soil temperature at soil freeze-up using both the GTN-P measurements and soil model simulations. The observed changes in the normalized ε with soil freezing is presented for a selected SoilSCAPE sensor node (S6) in Fig. 4b. ε below the freezing point and above −10 ^∘C ranges from ∼5 % to 20 % of the ε value for unfrozen soils. Assuming a linear relationship between ε and liquid water content (Mironov et al., 2010), Fig. 4b can approximate changes in unfrozen water content during soil freezing. Assuming soil freeze-up starts when ε drops below 15 %–20 % of the annual amplitude, the corresponding soil temperatures range from −0.01 to −1 ^∘C at depths between 0.05 and 0.56 m. We selected a temperature threshold of −0.35 ^∘C for soil freeze-up, which is at the higher end of the range indicated from the S6 node, but closer to the other SoilSCAPE nodes showing more rapid ε changes below 0 ^∘C. This temperature threshold is also consistent with our model simulations, which show a −0.3 to −0.5 ^∘C temperature range and 15 %–20 % liquid water content during freeze-up. This temperature threshold was used to determine the soil freeze onset and zero-curtain period at the GTN-P sites with only soil temperature measurements available. However, there is large variability in the relationship between ε and liquid water content at freezing temperatures due to changes in free and bound water and ice components (Mironov et al., 2010; Naeimi et al., 2012), which can result in large uncertainties in the above estimated thresholds.

2.5 Soil dielectric constant retrievals from airborne P-band radar

The soil model unfrozen water content curve used to define the soil freeze-up and zero curtain was only calibrated using limited SoilSCAPE soil dielectric measurements. We therefore evaluated the sensitivity of surface soil dielectric properties derived from local-scale (∼50 m) airborne low-frequency (P-band) radar acquisitions along regional transects in northern Alaska and associated F–T patterns to snow cover variations. Multiple flight lines were acquired in late August (fully thawed) and early October (partially frozen) of 2014 and 2015 in northern Alaska, including a regional transect along the Dalton Highway (DHN, 148.39–149.05^∘ W, 68.78–70.40^∘ N). Multiple soil parameters, including active-layer thickness (ALT) and soil moisture (converted from soil dielectric constant), were derived from NASA Airborne P-band (430 MHz) PolSAR radar backscatter measurements in August and October (Chen et al., 2019). The radar soil retrievals examined in this study differ from Yi et al. (2018), which used active-layer retrievals derived from multi-frequency (L + P-band) radar backscatter acquired in October 2015. The alternative single-channel (P-band) time series algorithm used in the current study avoids potential artifacts introduced from variable P and L-band radar acquisition times. The soil parameters were derived from the radar backscatter using a three-layer soil dielectric model. In August, the three layers represent the surface thawed layer, middle and bottom active layer, and the top of the upper permafrost layer. In October, the two surface layers represent a partially frozen active layer with a frozen surface layer overlying a deeper unfrozen active layer. An iterative optimization scheme was used to estimate the soil parameters by minimizing differences between the observed radar backscatter and radar scattering model simulations using the above three-layer soil dielectric model. The retrieved parameters include the thickness and soil dielectric constant of the surface layer in August and October, the depth and dielectric constant of deeper active layer. Initial validation indicated favorable radar ALT retrieval accuracy in relation to in situ CALM measurements along the DHN transect (Chen et al., 2018, 2019). The differences in surface soil dielectric constant between August and early October represent the changes in soil liquid water content that occur during F–T conditions. Therefore, we examined differences in soil dielectric constant under variable MODIS snow cover fraction to clarify the influence of snow cover on active-layer freezing.

2.6 Regional soil model simulation and analysis

We conducted an integrated analysis of in situ ground observations, soil model simulations and airborne radar retrievals to investigate the sensitivity of soil freezing to variable snow cover conditions across Arctic Alaska. For the regional simulation, the soil process model was spun up for 50 years to bring the top 10 m soil temperature profile into dynamic equilibrium using model inputs in the year 2000 (Yi et al., 2018), followed by a model transit run from 2001 to 2017. Unfrozen conditions in the deeper active layer may persist well into the winter season and into the subsequent calendar year. In order to accurately estimate the active-layer freeze onset and zero-curtain period for the current year, model simulations of the next calendar year were also needed. Therefore, the soil freeze onset and zero-curtain period in year 2017 were not calculated. The soil freeze onset for each soil layer was determined when model-simulated soil temperature dropped below −0.35 ^∘C as discussed in Sect. 2.4. The zero-curtain period at each soil depth was defined as the duration between land-surface freeze onset and freeze onset of the given soil layer. The regional correlation between snow onset calculated from the MODIS SCE data and the zero-curtain period for each soil layer was used to examine relations between the timing of early snow accumulation and soil freeze-up. Snow onset was determined as the center of the 8-day period with more than three adjacent snow-covered periods within a 40-day moving window; the relatively long temporal window was used to account for more variable snow cover conditions during fall.

We selected the DHN flight transect as the focus area for the integrated analysis due to the relatively dense network of GTN-P soil temperature and CALM ALT sites in this area relative to other transects (i.e., ATQ and IVO, Fig. 2). Because the MODIS 8-day SCE product (i.e., MOD10A2) only provides binary snow data (i.e., snow vs. nonsnow), the data set was binned for each 0.1^∘ latitudinal region along the radar flight transect to calculate the snow-covered area fraction (SCF). We also averaged the airborne radar-retrieved surface soil dielectric constant for each 0.1^∘ latitudinal bin and analyzed the correlations between SCF and selected variables representing the soil-freezing process, including the radar-retrieved surface dielectric constant changes between October and August, and the zero-curtain period derived from both the soil model simulations and in situ data. For the site analysis, the SCF from the 0.1^∘ latitudinal bin including the site was used. Arctic Alaska is generally fully snow covered by the end of October or early November. Therefore, we used the SCF averaged from September to October for the zero-curtain analysis, while only the September SCF was used for the airborne radar data analysis, since the radar data were obtained in early October.

3.1 Model validation

Accurate simulation of early cold-season soil freezing requires accurate characterization of landscape-scale snow cover conditions, which was addressed in this study by gap filling the MODIS SCE record to mitigate data loss from pervasive cloud cover and other factors. The gap-filled MODIS SCE products were then combined with other ancillary data to downscale the MERRA-2 reanalysis snow depth data, as one of the main driver data sets for the permafrost soil model. The accuracy of the gap-filled MODIS SCE product was cross-checked using the two MODIS sensors (Terra and Aqua); the downscaled snow depth data were evaluated using in situ SNOTEL observations across Alaska. The model-simulated soil-freezing process and ALT dynamics were conducted over a smaller Arctic Alaska domain and evaluated using a diverse set of regional observations.

3.1.1 Regional 1 km snow cover product

The cloud-free Aqua MODIS SCE data were used to evaluate the accuracy of filled pixels identified as cloud covered in the Terra MODIS SCE data and vice versa, assuming relatively consistent snow conditions between morning (Terra) and afternoon (Aqua) SCE acquisitions during the 8-day composite period. Our results indicate an accuracy of more than 80 % in the cloud-filtering algorithm with no obvious differences observed between the two sensor SCE records (Table 2 and Fig. S2). The cloud cover fraction for the 8-day temporal composite Terra MODIS SCE data represents 0.5 % to 10.1 % of the entire state throughout the year, and the percentage of cloud-free Aqua MODIS pixels that overlapped with cloud-covered Terra MODIS pixels ranges from 0.4 % to 4.6 % (Fig. S2). There are significantly more cloud-covered pixels in the Aqua MODIS record (1.1 %–15.2 %) and thus more cloud-free Terra MODIS pixels (0.9 %–9.8 %) overlapping with cloud-covered Aqua MODIS pixels. Cloud cover mostly occurs in the spring and fall shoulder seasons, resulting in larger SCE uncertainties during those periods. There is no obvious bias in the misclassification of cloud-contaminated pixels (Table 2), which indicates that using a cut-off threshold of 50 % for the snow occurrence probability (P_cutoff, Eq. 4 and Table 1) to classify snow or nonsnow conditions works well. Using a higher threshold (e.g., 60 %) generally results in more snow pixels misclassified as land pixels and vice versa.

Table 1Key model parameters used in this study.

Download Print Version | Download XLSX

Table 2The accuracy of the spatial filter algorithm applied to the Aqua MODIS SCE data during spring (April to June) and fall (September to November) transitions averaged over Alaska from 2003 to 2015. Pixels that identified cloud contamination in the Aqua MODIS record but indicated clear conditions in the Terra MODIS were used for evaluation. The percentage of cloud-contaminated and evaluated pixels were calculated for the entire Alaska domain, while the accuracy and misclassification were calculated as the percentage of the evaluated pixels. Both Terra and Aqua MODIS SCE data show similar accuracy, while only the Aqua MODIS results are shown here due to a higher percentage of cloud-contaminated pixels (available for evaluation) in the Aqua imagery.

Download Print Version | Download XLSX

Table 3Statistics of 1 km MERRA-2 snow depth data generated using different spatial interpolation schemes compared with in situ snow depth measurements at the Alaskan SNOTEL sites.

Download Print Version | Download XLSX

Compared with in situ snow depth measurements from the Alaskan SNOTEL sites, the 1 km MERRA-2 snow depth data generated using the new downscaling algorithm showed an overall improvement compared to the original spatial interpolation scheme used in Yi et al. (2018) (Table 3). The new 1 km snow depth data showed overall reduced RMSE and lower bias except in Interior Alaska at elevations between 400 and 800 m. At these elevations the USGS DEM used for spatial downscaling at the 1 km grid shows large deviations from the reported SNOTEL site elevations (Fig. S3), which may account for the relatively poorer performance of the new snow depth data set in this elevation band. In Arctic Alaska, the new snow depth product has a modest improvement over the Yi et al. (2018) product, with an RMSE of 0.16 m and bias of −0.01 m versus an RMSE of 0.18 m and bias of −0.03 m for the original data set. However, there are only eight SNOTEL sites in this region and only two sites at the Alaska North Slope. Compared with the Yi et al. (2018) product, the new MERRA-2 downscaled snow product captured more fine-scale details of spring snow melting and topographically varying snow distribution, especially in mountain areas (Figs. 5 and S4).

Figure 5Illustration of the snow data processing: the raw (a, d) and cloud-filtered (b, e) MODIS SCE images using an elevation-based spatial filter and downscaled MERRA-2 snow depth (in meters) data (c, f) using the filtered MODIS SCE and DEM data during spring snow melt (a–c: 23–30 April) and early snow accumulation (d–f: 30 September–7 October) period in 2007. In the MODIS images, snow-covered areas are shown in gray, while land and cloud-covered areas are shown in black and white, respectively.

Download

The snow offset and onset derived from the MODIS SCE and downscaled MERRA-2 snow depth records show very similar spatial patterns and trends over the 2001 to 2016 study period (Figs. S5 and S6). These results indicate that the downscaled MERRA-2 snow depth data generally capture the regional variability in snow cover conditions during the transitional season indicated from the MODIS observations. During the study period, the data sets show similar spring snow offset dates over Alaska (e.g., DOY 138±7 for MODIS versus DOY 140±7 for downscaled MERRA-2 snow depth), while there is an approximate 10-day difference in fall snow onset dates between the two data sets (MODIS: DOY 284±5; MERRA-2: DOY 273±5). The difference in mean snow onset is likely due to different methods used to determine snow onset for the two data sets. For the MODIS SCE record, snow onset was chosen as the composite period with more than three adjacent 8-day snow-covered periods within a 40-day moving window, while snow onset was defined from the downscaled snow depth record as the composite period with a mean snow depth above a 0.05 m threshold within a 24-day moving window. A higher snow depth threshold results in a later snow onset date in the MERRA-2 data set. However, the two records show similar snow onset patterns (Fig. S5) and interannual variability during the study period (R=0.79, p<0.01). Both data sets show a generally earlier snow onset trend in Arctic Alaska over the study period, which is discussed in Sect. 3.2.

Table 4Comparisons of model-simulated and in situ observed zero-curtain period at 0.35 m soil depth for sites along the DHN transect. The model-simulated mean zero-curtain period was calculated from 2001 to 2016. The zero-curtain period calculated using in situ soil temperatures at closely adjacent sites (e.g., the three Franklin Bluff sites, Table S1) were combined for a longer observational record.

Note that the asterisk denotes p<0.1, WD is West Dock, DH is Deadhorse, FB is Franklin Bluffs, SagMAT is Sagwon MAT, SagMNT is Sagwon MNT, HV is Happy Valley and GL is Galbraith Lake. Imnaviat 1 and GL were closely adjacent sites and were combined to form a longer time series.

Download Print Version | Download XLSX

Figure 6Comparisons of soil-freezing characteristics derived from GTN-P soil temperature measurements and model simulations along the DHN transect in northern Alaska: (a) interannual variations in zero-curtain period derived from in situ measurements, (b) variations of model versus in situ zero-curtain period at 0.35 m soil depth relative to MODIS SCF averaged for the northern part of the DHN transect, (c, d) changes in the soil freeze lag rate with depth at two DHN sites. Both sites have lower ALT in the earlier years shown here, ranging from 0.37 m in the earlier period to 0.43 m during later period of record for the WD1 site and from 0.47 to 0.57 m for the GL site.

Download

3.2 Sensitivity of active-layer freezing to snow cover

We first analyzed the sensitivity of active-layer freezing to seasonal snow cover using in situ observations, local-scale airborne P-band radar retrievals, MODIS snow cover observations and model simulations, focusing on the DHN transect in the Alaska North Slope. We then evaluated the regional sensitivity of model-simulated soil freeze onset and zero-curtain period within the active-layer to early snow accumulation indicated by the MODIS SCE record across Arctic Alaska.

3.2.1 Integrated analysis along the DHN airborne flight transect

For all sites along the DHN transect, both the model-simulated and in situ observed zero-curtain periods at 0.25 and 0.35 m soil depths showed significant positive correlations ($R=\mathrm{0.69}\pm \mathrm{0.14}$, p<0.1) with MODIS SCF. The zero-curtain period and SCF record showed similar interannual variability across the DHN transect, particularly north of Sagwon sites (>69.4^∘ N). Thus years with greater (less) snow cover are associated with a generally longer (shorter) zero-curtain period at these depths. The average zero-curtain period at 0.35 m soil depth for the DHN transect sites located from 70.4 to 69.4^∘ N is shown along with the corresponding MODIS SCF observations in Fig. 6b. These observations show large interannual variability in early season snow accumulation in this area, with the fall (September–October) SCF varying from 0.35 to 0.77 from 2001 to 2015. Correspondingly, both the in situ and model-simulated zero curtain show large variability throughout the study period; here, the in situ observations varied from 23.5 to 79.2 days at 0.25 m and from 24.7 to 90.3 days at 0.35 m, while the model simulations ranged from 29.3 to 76.0 and 44.0 to 90.7 days for the same soil depths.

The timing of early snow accumulation was the primary factor affecting the freezing process of the top soils, while ALT is more closely related to the length of zero-curtain period in deeper soils, particularly in areas with a shallower thaw depth. This can be seen from the delayed soil freezing below 0.30 m soil depth at two monitoring sites (West dock and Galbraith Lake) in years with larger ALT (Fig. 6c–d). Both monitoring sites show deeper ALT conditions during the later years of the study period, which are also associated with larger soil-freezing lag rates. Here, the soil-freezing lag rate is defined as the ratio of the soil freeze onset difference between two adjacent soil layers and the respective depth difference between the two layers. However, soil-freezing lag rates derived from both the in situ measurements and model simulations show large variability and are likely associated with large uncertainty, especially for deeper soil layers. The temporal coverage of the in situ observations used for estimating the soil parameter conditions is also less extensive than the model simulations, which may contribute additional uncertainty.

Figure 7Soil-freezing process indicated by the radar-retrieved (Chen et al., 2019) soil dielectric constant (ε1) of surface soils ($\sim <\mathrm{0.10}$ m) at the DHN transect in August (a) and October (b), and changes in ε1 in relation to MODIS SCF (c). The ε1 differences between August and October were binned to 0.1^∘ latitudinal bins, while SCF was calculated as the percentage of snow-covered pixels indicated by MODIS SCE data for each 0.1^∘ bin. The standard deviations of ε1 differences for each 0.1^∘ bin were shown as error bars.

Download

The airborne P-band radar retrievals over the DHN flight transect in early October (9 October in 2014, 1 October in 2015) showed a larger reduction in the surface soil dielectric constant (ε1) in areas with shallower snow cover during September than areas with deeper snow conditions (Fig. 7). In both years, the ε1 differences were negatively correlated with MODIS SCF across the DHN latitudinal gradient (2014: $R=-\mathrm{0.69}$, 2015: $R=-\mathrm{0.76}$, p<0.01, n=17). The reduction in ε1 retrievals in early October was more obvious in the northern part (>69.5^∘ N) of the transect (2014: $-\mathrm{10.35}\pm \mathrm{5.46}$, 2015: $-\mathrm{6.64}\pm \mathrm{1.94}$), which had a shallower seasonal snow cover (mean SCF of 0.22 in 2014 and 0.38 in 2015) and relatively early onset of frozen conditions than the southern part of the transect. Variation in vegetation cover may also contribute to the large contrast in the ε1 retrievals between the northern and southern portions of the transect, though the large apparent changes in ε1 between 2014 and 2015 in the northern portion of the transect provide a relatively robust indicator of frozen soil conditions. A similar relationship between radar-retrieved dielectric constant changes and MODIS SCF was observed over the Atqasuk (ATQ) flight transect. However, the opposite pattern was observed over the Ivotuk (IVO) transect, where ε1 changes increased with increasing seasonal snow cover at higher elevations (Fig. S9). The IVO transect differs from the other Alaska tundra transects by encompassing more variable topography (Fig. 2) and higher elevations (614±75 m) than both DHN (208±36 m) and ATQ (34±7 m).

Figure 8Model-simulated zero-curtain period at 0.35 m in relation to snow accumulation during the early snow season indicated by filtered MODIS SCE images for 2 selected years with later (2007: a–c) and earlier (2015: d–f) snowfall. In the MODIS images, snow was shown in dark gray, while land was shown in black and the areas masked out were shown in white.

Download

3.2.2 Model sensitivity analysis in Arctic Alaska

We first compared the model-simulated zero-curtain period over Arctic Alaska for 2 selected years with relatively late (2007) and early (2015) snowfall. These results indicate that snow accumulation during the early cold season is a primary control on the zero-curtain period within the upper (<0.4 m) soil layers (Fig. 8). The regional mean surface freeze onset and snow onset based on the MODIS LST and SCE observations was DOY 269±5 (surface freeze onset) and 281±9 (snow onset) in 2007, and DOY 259±8 (surface freeze onset) and 262±12 (snow onset) in 2015. The later snow cover establishment in 2007 resulted in an overall shorter zero-curtain period over most of the Arctic Alaska region, with a model-simulated mean zero-curtain period of 49.3±25.1 days at 0.25 m depth and 64.3±26.3 days at 0.35 m depth. In contrast, earlier snow accumulation in 2015 resulted in a longer zero-curtain period, ranging from 69.4±22.1 days at 0.25 m depth to 84.7±25.2 days at 0.35 m depth. The spatial pattern of the model-simulated zero-curtain period also corresponded well with the snow accumulation pattern indicated from the MODIS SCE observations, leading up to full snow cover conditions.

Figure 9Trends in model-simulated soil freeze onset (a) and zero-curtain period (b) at 0.35 m depth from 2001 to 2016. The fall snow onset trend for the same period derived from the MODIS SCE data is shown in (c).

Download

Figure 10Regional statistics of model-simulated zero-curtain period (a) and its sensitivity to MODIS snow onset and model-simulated ALT (b–c) from 2001 to 2016: (a) regional mean of model-simulated zero-curtain period at different depths; (b) changes in correlations between snow onset and zero-curtain period with depths; for both (a) and (b), the study area was divided into four groups: soil column froze below 0.33, 0.45, 0.55 and 0.7 m. The soil column for the majority of the study area froze below 0.7 m. (c) The proportion of pixels with a significant positive correlation between the zero-curtain period and ALT at different depths. The total number of unfrozen pixel was shown as “npixel”.

Download

During the study period, the spatial pattern of model-simulated soil freeze onset and zero-curtain period trends at 0.25 and 0.35 m is closely associated with the trends of MODIS snow onset (Fig. 9); areas with earlier snow onset generally show later freeze onset and a longer zero-curtain period. Further analysis indicates that early snow accumulation is the main control on the soil freezing for the upper ($\sim <\mathrm{0.4}$ m) active-layer, while the zero-curtain period of deeper soils is more related to ALT (Fig. 10). The model simulations show an overall longer zero-curtain period in areas with deeper ALT. Approximately 55 % of the Arctic Alaska has a maximum thaw depth between 0.45 and 0.55 m. In those areas, correlations between MODIS snow onset and zero-curtain period generally decrease with increasing soil depths below 0.4 m, which corresponds to an increasing positive correlation between the ALT and zero-curtain period. This reduction of the correlation between snow onset and zero-curtain period is more pronounced in areas with shallower ALT. Compared with the zero-curtain period, we found a much weaker correlation between model-simulated soil freeze onset in near-surface soils (<0.1 m) and MODIS snow onset but similar relationships for soil depths below (>0.2 m) (Fig. S10). This is due to a positive correlation between surface freeze onset derived from MODIS LST and MODIS snow onset (R=0.71, p<0.1) in Arctic Alaska during the study period; earlier surface freezing leads to cold underlying soil, while earlier snow onset generally leads to warm soil in this area. Therefore, soil freeze onset for near-surface soils has weaker correlations with MODIS snow onset compared with the middle of the active layer (∼0.2–0.4 m), while the soil freeze onset is subject to similar controls as the zero-curtain period with soil depth increases. These results indicate that the land-surface temperature or F–T status may be a relatively poor indicator of soil-freezing status and the zero curtain period in the deeper active layer for the Arctic region.

4.1 Sensitivity of active-layer freezing process to recent climate change

Our results show a strong correlation between the active-layer freezing process and snow accumulation during the early snow season in Arctic Alaska, especially within upper ($\sim <\mathrm{0.4}$ m) soil layers. Earlier snow onset and establishment of a complete snow cover generally delays active-layer freezing and promotes a longer zero-curtain period (Figs. 8 and 9). Previous studies have highlighted the decoupling of surface air and soil temperatures during the winter season in the northern high latitudes due to strong insulating effects of snow cover (Morse et al., 2012; Throop et al., 2012; Koven et al., 2013; Smith et al., 2016). Changes in the rate of accumulation, timing, duration, density and amount of snow cover play an important role in determining how soil F–T dynamics respond to surface warming (Zhang, 2005; Lawrence and Slater, 2010). The relationship between fall snow onset and soil warming may vary depending on the timing and magnitude of snowfall, and local climate conditions (Yi et al., 2015). Early snow onset may enhance thermal buffering of cold surface temperatures and promote soil warming in colder climate zones such as Arctic Alaska. A shorter snow season may cool the soil in colder areas due to less insulation from cold atmospheric temperatures but may warm the soil in warmer areas by promoting greater heat transfer into soils (Lawrence and Slater, 2010). The snow cover impact on soil F–T dynamics will also depend on differences between the timing of first surface freeze and snow cover establishment, especially for near-surface ($\sim <\mathrm{0.1}$ m) soils (Kim et al., 2015).

Our model simulations also show that the influence of snow cover on soil freezing is weaker for deeper soil layers ($\sim >\mathrm{0.5}$ m), where the freeze onset and zero-curtain period are more closely related to the summer maximum thaw depth (i.e., ALT, Figs. 10 and S10). This can be largely explained by a close link between ALT and the soil-freezing lag rate of deeper soils. During fall, active-layer freezing proceeds both downward from the surface and upward from the underlying permafrost table (Outcalt et al., 1990; Oechel et al., 1997; Zona et al., 2016). This can be seen from the negative soil-freezing lag rate (related to differences in soil freeze onset between two adjacent soil layers) at the bottom of the active layer at the GTN-P sites (Fig. 6c–d), indicating that the bottom of the active layer freezes first. Increases in ALT can lead to abrupt changes in the soil-freezing lag rate at the same soil depth, which can change from a negative value to a small positive value and result in abrupt changes in the zero-curtain period at the bottom active layer (e.g., at the WD1 and WDN sites in Fig. S8b). Previous studies have also reported a delayed soil freeze-up and thus longer zero curtain with increasing ALT (Morse et al., 2012; Euskirchen et al., 2017). Based on the GTN-P site measurements, deeper soils show a general delay in soil freeze onset relative to shallower active layer, with a mean lag rate of 0.79±0.52 days cm⁻¹ at 0.35 m depth; large variability in the soil-freezing lag rate is likely associated with different soil structure and variable active-layer soil moisture content (Throop et al., 2012). Therefore, a deepening active layer associated with climate warming will likely lead to a longer zero-curtain period in deeper soils.

The Arctic is expected to experience continued warming and precipitation increases under projected climate trends (IPCC, 2013), though the potential response of active-layer freezing to these changes may vary depending on changes in seasonal snow cover. Both surface warming and a changing precipitation regime can modify seasonal snow cover, leading to a nonlinear response of soil temperatures to warming (Lawrence and Slater, 2010; Yi et al., 2015). Increases in winter precipitation and snowpack deepening may enhance soil warming, while a reduced snowpack may promote soil cooling in colder climate areas. More frequent and intense rain-on-snow events during fall and early winter have been observed across the ABZ with recent warming trends (Ye et al., 2008; Langlois et al., 2017). Therefore, how these climate trends will affect soil moisture and thermal dynamics is a key challenge for accurately estimating soil F–T dynamics and potential carbon and climate feedbacks. In addition, with continued warming and ALT deepening, unfrozen conditions may persist in the bottom of the active layer, resulting in a perennially thawed subsurface soil layer or talik zone; once a talik forms, it can greatly accelerate permafrost degradation and result in large changes in surface hydrology and soil carbon decomposition (Yoshikawa and Hinzman, 2003; Parazoo et al., 2018).

Our model simulations are associated with large uncertainties, particularly regarding soil moisture effects on soil heat transfer during the F–T period. Changes in liquid water content during soil freezing vary for different soil conditions, while accurate simulation of this process is challenging due to complex processes controlling ice formation, liquid water movement and heat transfer in frozen soils (Outcalt et al., 1990; Romanovsky and Osterkamp, 2000; Schaefer and Jafarov, 2016). Our study used in situ soil dielectric constant (ε) measurements to parameterize the unfrozen water curve and determine the temperature threshold used to define soil freeze onset and calculate the zero-curtain period. However, our results also show large ε variability in response to freezing temperatures from the in situ SoilSCAPE measurements; the relationship between ε and liquid water content in organic-rich soils may also be different from mineral soils (Engstrom et al., 2005; Mironov et al., 2010). A reliable soil dielectric model characterizing the relations between unfrozen water content and ε for organic soils will help reduce the uncertainty in the estimated temperature threshold at soil freeze onset; a model sensitivity analysis using different b values in the unfrozen water curve (Eq. 11) may also help quantify uncertainties in the model-simulated zero-curtain period and its regional pattern. On the other hand, potential soil moisture redistribution with active-layer deepening is not accounted for in the current model, though this effect is likely small due to small ALT trends in this area during the study period (Yi et al., 2018). Increasing disturbance from thermokarst and wildfire are expected to alter microclimate and soil moisture conditions, vegetation cover and SOC stocks in the ABZ, which will also likely influence the dynamics of ground-ice evolution and permafrost degradation (Grosse et al., 2011; Liljedahl et al., 2016).

4.2 Potential use of remote-sensing to improve regional monitoring of soil F–T process

Large-scale satellite observations and global reanalysis data generally have difficulty capturing finer-scale snow cover variations and associated impacts on soil F–T dynamics. These limitations are exacerbated in the Arctic due to a paucity of regional climate stations and complex microclimate and snow cover conditions influenced by local topography, vegetation and winds (Liston and Sturm, 2002; Gisnås et al., 2016). Optical satellite remote sensing, including Landsat and MODIS sensors, can provide accurate local-scale information on the snow cover extent, though effective regional monitoring from these observations is constrained by persistent cloud cover and reduced solar illumination for much of the year. Moreover, these observations do not include snow depth or water content information, which are critical parameters for hydrologic and ecological applications (Brown et al., 2010; Painter et al., 2016). Snow-covered areas attenuate the emitted microwave radiation from the underlying surface, while the magnitude of microwave emissions and attenuation depends on the sensor frequency, snow liquid water content, snow grain size and shape. Thus, snow properties, including snow water equivalent (SWE), may be derived from passive microwave sensors, albeit at relatively coarse spatial scale (Kelly et al., 2003; Armstrong et al., 2005; Takala et al., 2011). However, its accuracy is limited in deep snowpack conditions, and its applicability is limited in forest areas and wet snow conditions (Frei et al., 2012). Compared with passive microwave sensors, active radars or scatterometers such as Ku band are capable of much higher spatial resolutions and can be particularly useful for regional snow mapping (Yueh et al., 2009). However, more studies are needed to clarify the multiple scattering effects from snow microstructure variations and contributions from other elements within the sensor footprint including vegetation, soil and open-water effects, to ensure accurate retrieval of snow properties (King et al., 2018). Airborne laser altimeters (lidar) also show strong potential for mapping snow depth patterns (Deems et al., 2013; Painter et al., 2016), while the recently launched ICESat-2 is expected to provide new capabilities of satellite lidar for regional snow mapping (Kwok and Markus, 2018). In the near term, significant improvements in acquiring geospatial information on snow properties in the Arctic will likely come from merging in situ and modeling data sets with multi-sensor snow products (Painter et al., 2016).

Another major challenge for regional permafrost modeling is the lack of information on subsurface properties, particularly for organic soils with distinct hydraulic and thermal properties. Current permafrost models generally use regional or global SOC inventory data to parameterize the SOC profile following an exponentially decreasing curve (Lawrence and Slater, 2008; Rawlins et al., 2013; Yi et al., 2018). However, large discrepancies are apparent from the available SOC inventory records in the ABZ (Liu et al., 2013; Hugelius et al., 2014). There is also large regional variability in the vertical SOC distribution due to multiple processes affecting the SOC distribution in cryoturbated soils (Mishra et al., 2013; Hugelius et al., 2014; Hossain et al., 2015). Long wavelength radar, including P-band (∼70 cm) and L-band (∼24 cm), is sensitive to surface vegetation structure, soil surface and subsurface dielectric properties (e.g., Dobson et al., 1985; Mironov et al., 2010; Tabatabaeenejad et al., 2015). Our model experiments and analysis using in situ dielectric constant measurements and a P-band radar-retrieved soil dielectric constant for soil-freezing studies, albeit simple, show the potential of longwave radar remote sensing in mapping of SOC, active layer F–T and moisture profiles. However, similarly to many other inversion problems, radar retrievals suffer from ambiguity in the inversion parameter definitions, mainly due to insufficient information on the subsurface profile (e.g., Tabatabaeenejad et al., 2015; Chen et al., 2019). Therefore, new methodologies are needed to address the underdetermined nature of the radar backscatter inversion and associated land parameter retrievals, by either including additional observations or other synergistic information from soil physical models to reduce parameter dimensions in the radar model (e.g., Sadeghi et al., 2016). The vegetation canopy also has a large impact on the radar backscatter, especially at L-band and shorter wavelengths, while separating the radar contribution of subsurface soils from the vegetation canopy remains a challenge. Radar measurements with shorter wavelengths (e.g., X, C, L-bands) can be also useful for subsurface retrievals by providing contributing information on surface snow, soil and vegetation conditions (Moghaddam et al., 2000; Yueh et al., 2009; King et al., 2018), which can be used to reduce the uncertainties in the longwave (e.g., P-band) radar soil parameter retrievals. However, more sophisticated modeling experiments capable of representing complex landscapes and multi-frequency radar backscatter characteristics are needed to fully clarify the value of multi-frequency observations. Additional airborne radar sampling targeting regional disturbance gradients may also provide the necessary information for the regional modeling framework to represent increasing disturbance regimes and associated impacts on active-layer F–T dynamics in the ABZ.