A model to characterize soil moisture and organic matter profiles in the permafrost active layer in support of radar remote sensing in Alaskan Arctic tundra

In northern circumpolar regions, it is estimated that 70% of the total carbon (1000 Pg) is stored in the top 3 m of ground underlain by permafrost. (Tarnocai et al 2009, Hugelius et al 2014). Due to the current warming trends and the resulting impact on both anaerobic and aerobic soil respiration, the sequestered soil organic carbon (SOC) in permafrost could be released into the atmosphere as CH₄ and CO₂ (Mishra et al 2013, Schuur et al 2015, Abbott et al 2016). In addition to the impact of increasing temperatures, soil moisture variability also plays a major role in the loss of carbon from the permafrost SOC to the atmosphere (Natali et al 2015). In a recent study, land surface modelers prioritized SOC and soil moisture as key variables for studying Arctic ecosystem dynamics in a changing climate (Fisher et al 2018). Within US soils alone, observational SOC estimated from soil pedon data indicates that Alaskan soil accounts for up to half of the stored total US soil carbon (Bliss and Maursetter 2010). However, the vastness, remoteness, and harsh environment in the Arctic region have led to sparse sampling and an often incomplete and poorly constrained carbon inventory in Alaska (Ping et al 2008).

In permafrost regions, almost all biogeochemical activities occur in the active layer, the soil above permafrost that freezes and thaws seasonally (Hinzman et al 1991). The permafrost table restricts vertical water flow and maintains a shallow water table and a highly saturated soil across the active layer in the Arctic. Furthermore, cold temperatures and shallow water tables lead to a low SOC decomposition rate and subsequently cause an accumulation of high organic matter (OM) content in active layer soils (Chen et al 2019a). Pedon delineation and sampling processes designed for permafrost affected soil often collect samples from representative horizons (Ping et al 2013, O'Connor et al 2020). This horizon-based sampling results in just a few samples within the active layer with coarse sampling depths. Therefore, in part, some of the critical processes associated with variations of soil OM (SOM) properties over short vertical distances through the active layer soil profile might be overlooked (Harden et al 2012, Manies et al 2020). Consequently, a more detailed study of soil properties is of interest to model the fine-scale vertical profile and, more specifically, the SOM profile through the permafrost active layer (Jobbagy and Jackson 2000).

Soil moisture variations in the active layer are a primary driver of Arctic carbon exchange, active layer thermal conductivity, and energy transfer affecting permafrost stability (Lawrence et al 2015), yet in-situ soil moisture data in the Arctic and Boreal regions of North America is quite sporadic. In-situ observations have been collected either at permanently installed sites with near-real-time soil moisture data at several depths within the active layer, or manually during field campaigns (Nicolsky et al 2009, Biskaborn et al 2015). While the permanently installed sensors provide sufficient temporal resolution, they suffer from lack of spatial distribution and a limited number of sensors per profile that do not capture active layer fine-scale dynamics (Dafflon and Torn 2017, Chen et al 2018, Romanovsky et al 2020). The second observation method, employed during limited field campaigns, uses a soil moisture probe, ground penetrating radar, and/or gravimetric sampling of soil water content. While these samples can be taken with high spatial density; the total spatial coverage remains rather limited. These field campaign measurements are temporally sparse and comprise single (or multiple) point(s) measurements in time and space. Furthermore, even the largest available soil moisture dataset collected within the Alaskan tundra provides only the average values of soil moisture from probes inserted vertically into the ground in conjunction with GPR measurements (Bourgeau-Chavez et al 2010, Clayton et al 2021). The challenge still remains, in that there is a lack of fine-scale vertical profile characterization of soil moisture throughout the active layer.

A major user of fine-scale vertical profile models characterizing SOM and soil moisture behavior through the active layer is the radar remote sensing community. For example, soil moisture profile models and detailed soil texture maps available for temperate sites within the contiguous US have been used for retrieval of subsurface soil moisture from airborne low frequency radar observations (Tabatabaeenejad et al 2015, Tabatabaeenejad et al 2020, Cuenca et al 2016). However, in microwave remote sensing of permafrost active layer soils, neither soil moisture nor SOM profile models have been quantified (Chen et al 2019b, 2019c). Such models could also be incorporated into land surface models, replacing current approximate SOC profile models being used in numerical simulations, such as the exponential form used by Yi et al (2018).

The work reported here directly supports the NASA airborne campaigns within the Arctic Boreal Vulnerability Experiment (ABoVE) (Miller et al 2019), by providing a detailed experimentally derived model to describe soil moisture and OM distribution in the active layer soil profiles. We further develop a pedo transfer function (PTF) that relates soil mineral texture and OM to bulk density and porosity and can be used to fill the data gap in existing databases, in which critical information such as bulk density is often lacking (Johnson et al 2011, Mishra and Riley 2012). The dataset used in this paper was obtained in August 2018, through a series of field experiments conducted at five tundra sites located along the Dalton Highway in northern Alaska (Bakian-Dogaheh et al 2020). The sites are located within two of the ABoVE Airborne Campaign (AAC) flight lines, commonly referred to as the 'Deadhorse' and 'Toolik' lines.

2.1. Site description

The Deadhorse and Toolik flight swaths of the AAC (figure 1) are among the most widely studied areas in the North Slope of Alaska. The vicinity of these lines along the Dalton Highway and the Toolik Field Station are densely sampled areas and a benchmark for Arctic ecosystem studies. During the 2018 field surveys, we selected five study sites along the Dalton Highway to acquire field measurements within the airborne synthetic aperture radar (SAR) footprints of the 2017 AAC. The sites include Franklin Bluffs (FB), Sagwon (SGW), Happy Valley (HV), and Ice Cut (ICC), which are located within the Deadhorse flight line, and Imnavait Creek (IMN) within the Toolik flight line. The sites were chosen after close coordination with the ABoVE field scientists. FB, SGW, and ICC were chosen because there had been multiple soil coring experiments the year before by other ABoVE field scientists. The SGW site was chosen based on preliminary results from the ABoVE radar team suggesting an unexpectedly deep active layer thickness (ALT) in the radar retrievals. The HV site was chosen because of its proximity to the University of Alaska Fairbanks borehole temperature monitoring system (Wang et al 2018), the University of Southern California (USC) Arctic Soil Moisture Sensing Controller And oPtimal Estimator (SoilSCAPE) site (Chen et al 2018), and Alaska Circumpolar Active Layer Monitoring sites (Brown et al 2000). Also, there had been multiple soil coring efforts in previous years around the area. The presence of multiple eddy covariance flux towers in the IMN site was another factor in choosing the site.

Zoom In Zoom Out Reset image size

Further consideration and deciding factors for the site (figure 1 and table S1 available online at stacks.iop.org/ERL/17/025011/mmedia) selections involved capturing a wide range of land cover types and permafrost soil features specifically related to SOM content. The National Land Cover Database (NLCD) map and in-situ observations show that the IMN site is dominated by tussock graminoid tundra (tussock sedge, dwarf shrub, moss tundra). The HV and ICC sites are covered mostly with erect dwarf shrub tundra (low shrub tundra), while the SGW site is characterized by non-tussock graminoid (non-tussock sedge and dwarf shrub) or wetland (wet sedge) tundra (Stow et al 2004, Walker et al 2005, 2012, Epstein et al 2012, Raynolds et al 2012, Gagnon et al 2019).

2.2. Field experiment design

2.2.1. Soil moisture measurement

The in-situ determination of soil moisture content at the study sites was made using direct and indirect methods. Most conventional indirect methods for in-situ soil moisture measurement are electromagnetic techniques that use either capacitive measurement or time domain reflectometry (Dorigo et al 2011). These methods, which exploit the relationships between soil dielectric permittivity, soil moisture, and other physical properties of soil, always need further calibration (via a soil dielectric mixing model) to translate the measured capacitance or time delay of the permittivity measurement to the soil water content. On the other hand, the direct methods of soil moisture observation rely on field sampling of moist soil and finding the soil moisture value via over drying (this is often referred to as gravimetric water content and is the gold standard).

In this work, we used both measurement types. For the direct soil moisture measurement, two side-by-side (adjacent) soil samples were extracted from each layer (figure 2), and moist soil mass was recorded for each sample in the field to measure the in-situ soil moisture content ($\theta$). The gravimetric soil moisture content was calculated via oven drying (table S2). For the indirect method (dielectric probing) the values of dielectric constant were not translated into soil moisture because of the need for an organic soil dielectric model, which will be reported in our future work.

Zoom In Zoom Out Reset image size

For dielectric probing, the previous measurement method (Clayton et al 2021, Schaefer et al 2021) did not capture the profile behavior and only provided average information about near-surface or total active layer soil moisture. Our dielectric measurements were taken along the soil profiles using METER Group's TEROS 12 sensors (TEROS 12 manual 2020) and the ProCheck data logger. These sensors operate at 70 MHz and measure the dielectric permittivity and electrical conductivity. Measurements were taken at multiple discrete points along the active layer soil profile. Soil pits were dug, and sensors were inserted horizontally into the soil to collect the measurements (figure 2). This horizontal sampling was performed along the inner wall of the soil pits, with a sampling interval of 2.5 cm from the surface down to the top of the permafrost table, located at the base of the active layer up to 55 cm below the surface at the tundra sites. This approach provides the opportunity to capture the vertical variations of dielectric permittivity at a fine resolution throughout the active layer profile, which is not possible in alternate methods where probes are inserted vertically or at an angle. Furthermore, while vertical sampling may be faster in a field survey, it only measures the effective dielectric properties along the length of the probe from the surface, which cannot capture the subsurface features that are observable by low-frequency radars such as P-band (Bourgeau-Chavez et al 2010).

2.2.2. Soil property measurement

The soil sampling procedure was designed to obtain a full characterization of the subsurface profile. Conventional field protocols and soil sampling procedure in the Arctic Tundra region often consists of single coring or opening a soil pit and sampling from representative horizons, which results in a coarse characterization of the active layer (2–3 samples). (Jobbagy and Jackson 2000, Ping et al 2008, 2013, Tarnocai et al 2009, Agus et al 2010, Johnson et al 2011, Mishra and Riley 2012, Michaelson et al 2013, Hossain et al 2015, O'Connor et al 2020). The majority of the previous studies were conducted by researchers focused on Earth system and/or hydrology models. There is a fundamental difference between the way subsurface is characterized in Earth system models and physics-based radar remote sensing models. In the latter, a coarse horizon-based soil sampling may not capture the transition and details of the soil property variations in the active layer, where, despite its relatively shallow depth, its soil moisture and OM vary substantially. While we mainly follow the sampling protocols suggested by previous methods, the main difference in our method is fine resolution sampling to acquire more sampling points in the subsurface to characterize the active layer soil profile model. Otherwise, the transition points between each horizon will be ignored. Therefore, after opening each soil pit, two side-by-side tin cans of soil samples were harvested from the sidewall. Usually, the samples closer to the surface, which typically have more organic, are done at depths of 7–8 cm, and the deeper layer that have more minerals content are sampled at intervals of around 5–6 cm (table S3).

Furthermore, some pieces of information, such as bulk density (${\rho _{\text{b}}}$), porosity ($\varphi$), and root biomass (RB), are often missing in datasets. Sample containers with the same volume (V = 251 cm³) were used so that ${\rho _{\text{b}}}{ }$ and $\varphi$ could be measured in the lab. Wet soil mass was recorded in the field to measure the in-situ soil moisture content ($\theta$). The root layer decomposed organic layer, and the mineral layer are visually distinguishable in the two soil pits excavated at HV and ICC (figures 3 and S5).

Zoom In Zoom Out Reset image size

2.3. Lab experiment design

3.1. Pedotransfer function

Numerical coupled thermal-hydrologic processes used in land surface models and electromagnetic scattering models used in radar remote sensing usually represent the subsurface as a multi-layered structure. The subsurface thermal, hydraulic, and dielectric properties need to be parametrized for numerical simulations, and these properties are strongly dependent on soil moisture variation and texture (mineral, organic, root layer) (Tabatabaeenejad and Moghaddam 2006, Lawrence and Slater 2008, Bircher et al 2016). Generalized models for subsurface parameterization often rely on soil properties, such as mineral texture (sand (S) and clay (C) fractions), OM (which encompasses SOM and RB), bulk density (${\rho _{\text{b}}}$), and porosity ($\varphi$) (table S1). These soil properties serve as inputs for these parametrization models. PTFs leverage the strong correspondence between different soil properties and provide a means of reducing the complexity of these parametrization models by reducing the number of input data layers. The goal here was to find a function that estimates ${\rho _{\text{b}}}$ and $\varphi$ from soil texture (OM, sand, and clay fractions) through simple empirical modeling.

Important points that must be made are the definitions for organic and mineral soil. Our field observations indicate 35% OM (equivalent to 20% carbon content) as distinguisher threshold between mineral and organic soils (Agus et al 2010, O'Connor et al 2019, Manies et al 2020). Above this threshold, the properties of soil are governed only by OM (Manies et al 2020). This is consistent with our measurement observations, that for organic samples (OM > 35%) even combinations of the two adjacent samples did not have sufficient mineral constituents and was dominated by organic particles. Therefore, mineral texture of soil with OM beyond this range has minimal influence on the pedotransfer modeling. Subsequently, in our modeling of samples with OM > 35% we estimated the texture with the mean value of each component (sand and clay) due to almost-uniform behavior of mineral texture throughout the profile.

Our dataset provided 102 samples from eight soil pits including replicates, which are referred to as sample 'a' and 'b' and can be treated as an independent sample for PTF development (figure 2, table S3). We combined our data with another recently published dataset collected within the Toolik flight line (O'Connor 2019), which resulted in an expanded set of 214 samples covering a wide range of OM distribution (figure S1), from mineral to organic soil. It is important to note that O'Connor et al (2019) sampling follows a coarse characterization of the active layer; nevertheless, the data can be used for PTF model development. Unlike conventional PTFs that are based on data-driven regressions, herein we focused on correlation-based models. Here, a data-driven regressions PTF refers to multi-variable models that use machine learning techniques such as Random Forest to find a PTF model. Correlation-based models are parsimonious models with great explanatory predictive power. For instance, in the following we estimate different soil properties from soil textures.

Analysis of our data indicates that the bulk density for mineral fraction can be found as follows:

$$\begin{equation}{\rho _\text{b}}^{\text{Min}} = {e^{0.0056S + 0.0081C}}\,\,\,\,\,\,\,\,\,\,\,\,\end{equation} \tag{ 1 }$$

where ${\rho _{\text{b}}}^{{\text{Min}}}$ is the mineral bulk density associated with mineral texture fraction. Sand and clay mass fractions (%), which serve as input variables are denoted by (S) and (C) respectively.

Furthermore, knowing that pure organic soil often shows very small bulk density, to find the total bulk density we used an exponential function as follows:

$$\begin{equation}{\rho _{\text{b}}}\, = \,\rho _{\text{b}}^{{\text{Min}}}{e^{ - 0.0314{\text{OM}}}}{ }{\text{}}.\end{equation} \tag{ 2 }$$

Accordingly, once the value of ${\rho _{\text{b}}}$ was found based on OM, S, and C, to find the soil sample porosity we adopted a linear relationship between $\varphi$ and ${\rho _{\text{b}}}$ (O'Connor et al 2020).

$$\begin{equation}\varphi \, = \, - 0.3771{\rho _{\text{b}}}\, + \,0.9247{ }{\text{}}.\end{equation} \tag{ 3 }$$

The above relationship was confirmed from our site soil measurements, with parameters calculated from data.

3.2. Organic matter profile model

3.2.1. OM profile

Both soil moisture and OM vary substantially along the active layer profile. Therefore, an accurate empirical characterization to describe this profile behavior is necessary. For the purpose of profile characterization, we used only our fieldwork soil samples; the O'Connor et al (2020) samples were not included due to their sporadic sampling points within the profile. The OM distribution for depth intervals of 5 cm throughout the profile (figure S2(a)) and the median of OM content at each corresponding depth (figure S2(b)) show that the OM profile behavior is effectively represented using a logistic sigmoid (S-shaped) function.

$$\begin{equation}{\text{OM}}\left( z \right)\, = \,{\text{O}}{{\text{M}}_{z{\text{0}}}}\, + \,\frac{{{\text{O}}{{\text{M}}_{\text{M}}} - {\text{O}}{{\text{M}}_{z{\text{0}}}}}}{{1 + {e^{ - \beta \left( {z - m} \right)}}{ }}}{ }{\text{}}.\end{equation} \tag{ 4 }$$

Equation (4) explains the behavior of the OM(z) profile, where $0\, \leqslant \,{\text{O}}{{\text{M}}_{\text{M}}}\, < \,35\%$ is the OM fraction in the mineral horizon, $0\, \leqslant \,{\text{O}}{{\text{M}}_{z0}}\, \leqslant \,100$ is the surface OM, $\beta > 0$ is the decay rate, and $m > 0$ indicates the depth where the maximum decay rate occurred. Based on the 35% OM criterion, the organic layer thickness (${z_{{\text{OLT}}}}$) can be found for different soil pits. This S-shaped profile is different from previous observations, in which the OM was expressed as an exponential function (Meersmans et al 2009, Hossain et al 2015). The difference in profile functional forms reflects the more detailed sampling procedure used in this study, which provides finer resolution of profile gradients than more conventional bulk sampling from representative soil horizons (Chen 2019) (figure S3).

3.2.2. RB and SOM profile

Aboveground vegetation biomass plays a crucial role in governing subsurface hydraulic, thermal, and dielectric characteristics (Hinzman et al 1991). While between 70% and 80% of the vascular plant biomass in Arctic tundra is located below ground, there is still only limited knowledge about the below-ground responses of tundra ecosystems to climate change (Mokany et al 2006, Poorter et al 2012, Iversen et al 2015, Wang et al 2016). One conventional approach to quantifying the RB distribution is finding a relationship between roots (underground) biomass and shoot (aboveground) biomass, namely, the root:shoot ratio. Species dependent root:shoot and other allometric equations are poorly constrained for tundra (Mokany et al 2006, Wang et al 2016). Our detailed soil analysis based on particle dimensions also provides detailed information on the RB profile, RB(z), and SOM profile, SOM(z). As described in equation (5), the total subsurface OM profile OM(z) at each sample site is associated with different subsurface constituents, including RB(z), SOM(z), gravel fraction (GF(z)), and mineral fraction of the soil through mass fraction (table S1).

$$\begin{align}{\text{OM}}\left( z \right) = {\text{RB}}\left( z \right) + {\text{SOM}}\left( z \right)\left( {1 - \frac{{{\text{RB}}\left( z \right) + {\text{GF}}\left( z \right)}}{{100}}} \right){ }{\text{}}.\end{align} \tag{ 5 }$$

RB(z) and SOM(z) can be derived from the in-situ observations. The best profile model for RB(z) is an exponential function, as shown in equation (6), where $0\, \leqslant \,{\text{R}}{{\text{B}}_{z0}}\, \leqslant \,100\%$ is the surface RB and $\kappa \leqslant 0$ is the exponential decay factor. SOM(z) follows an S-shaped behavior similar to OM(z), where the corresponding parameters are also similar to what was described in equation (4).

$$\begin{equation}{\text{RB}}\left( z \right) = {\text{R}}{\text{B}_{z0}}{e^{\kappa z}}\,\end{equation} \tag{ 6 }$$

$$\begin{equation}{\text{SOM}}\left( z \right)\, = \,{\text{SO}}{{\text{M}}_{z0}}\, + \,\frac{{{\text{SO}}{{\text{M}}_{\text{M}}} - {\text{SO}}{{\text{M}}_{z0}}}}{{1 + {e^{ - \beta \left( {z - m} \right)}}}}{ }{\text{}}.\end{equation} \tag{ 7 }$$

Note the two different expressions for OM(z) in equations (4) and (5). In reality, we measure SOM(z), RB(z), GF(z), and find the total OM(z) based on equation (5). However, one can also fit an S-shaped behavior for total OM as described in equation (4). Among the 8 soil pits, except for SGW-1, the rest show negligible gravel fraction (figure S4). In equation (4), which we refer to as model 1, for convenience we ignore the gravel fraction (GF(z) = 0). A comparison between these two models is depicted in the results section (figures 5 and 6). For the purpose of the radar retrievals, we use equation (4) as it requires fewer parameters to characterize the total OM profile.

3.3. Soil moisture profile model

Because the porosity ($\varphi$), or maximum saturation level, substantially changes within different active layer soil horizons, the saturation water fraction (${\text{SW}}$) is a better representation for soil moisture ($\theta$), where it can be denoted as a normalized soil moisture level. This also allows for estimating the water table depth (${z_{{\text{WT}}}}$), where ${\text{SW}}$ reaches 1. The soil moisture profile $\theta \left( z \right)$ can be found as the product of SW fraction and porosity:

$$\begin{equation}\theta \left( z \right)\, = \,{\text{SW}}\left( z \right) \times \varphi \left( z \right){ }{\text{}}.\end{equation} \tag{ 8 }$$

The porosity profile $\varphi \,\left( z \right)$ is found by substituting equation (4) into equation (2), and knowing the soil mineral texture; then the resulting ${\rho _{{\text{b }}}}\left( z \right)$ profile is applied to equation (3) to find $\varphi \,\left( z \right)$. Our field observations show the best fit to characterize the ${\text{SW}}\left( z \right)$ is represented by a quadratic function:

$$\begin{align}{\text{SW}}\left( z \right)\, = \,\left\{ {\begin{array}{*{20}{c}} {1 - \left( {1 - {\text{S}}{{\text{W}}_{z{\text{0}}}}} \right){{\left( {\frac{z}{{{z_{{\text{WT}}}}}} - 1} \right)}^2}}&{z\, \leqslant \,{z_{{\text{WT}}}}} \\ 1&{z\, > \,{z_{{\text{WT}}}}} \end{array}} \right.{ }\end{align} \tag{ 9 }$$

where ${\text{S}}{{\text{W}}_{z0}}$ indicates the surface SW fraction, which varies between 0 and 1, and ${z_{{\text{WT}}}}$ is the water table depth.

Soil mineral texture and organic properties serve as inputs of the PTF model (figure 4). We assumed an average value for the mineral texture, with ${S_{{\text{avg}}}}\, = \,39\%$, and ${C_{{\text{avg}}}}\, = \,25\%$. The model showed close agreement with the measurements and was able to estimate ${\rho _{b\,}}\left( z \right)$ with a root mean square error (RMSE) of 0.15 (${\text{g}}\,{\text{c}}{{\text{m}}^{ - 3}}$) (figure 4(a)). Once the bulk density was found, it was used in equation (3) to reconstruct soil porosity with RMSE of 0.05 (${\text{c}}{{\text{m}}^3}\,{\text{c}}{{\text{m}}^{ - 3}}$) (figures 4(b) and (c)). The FB samples were outliers in the PTF development, the major feature of the FB soil is a thin RB layer with a homogeneous mineral texture profile having an average sand fraction around $60\% - 70\%$ and OM less than $35\%$. Therefore, the PTF model developed in section 3.1 and the corresponding results suggest that the model is overestimating ${\rho _{b\,}}$ for mineral soil that contains a high sand fraction (figures 4 and 5(e) FB-1).

Zoom In Zoom Out Reset image size

Next, we show the in-situ measurement results and profile model behavior for all of the soil pits in tundra sites (figure 5), which are the OM, ${\rho _{\text{b}}}$, and $\varphi$.

Zoom In Zoom Out Reset image size

At each depth, the measured value for adjacent samples 'a' and 'b' (figure S2) is reported. In most cases, the two measurements are consistent. However, due to microtopography variation within a site, profile behavior at each soil pit is distinct from other neighboring pits at the same site. Here, the${\text{ OM}}\left( z \right)$ profile follows the suggested S-shaped (sigmoid function) behavior in both soil pits (figures 5(a) and (d)). For the SGW-1 and HV-1 (compare SGW-1 (red) and HV-1 (black) in figure 5(a)) locations, the organic layer is shallower, and ${\text{OM}}\left( z \right)$ rapidly decreases from the surface to the mineral layer. However, the SGW-2 and HV-2 (compare SGW-2 (blue) and HV-2 (magneta) in figure 5(a)) and IMN-1, ICC-1 and ICC-2 (compare IMN-1 (blue), ICC-1 (black), and ICC-2 (magneta) in figure 5(d)) locations include a thicker organic layer of about 20–30 cm. FB shows a fairly homogenous mineral soil throughout the profile (see FB-1 (red) in figure 5(d)). The corresponding ${\rho _{\text{b}}}$, and $\varphi$, profiles were acquired by applying the ${\text{OM}}\left( z \right)$ profile from equation (4) and inserting it into the PTF models (figures 5(b), (c), (e), and (f)). In almost all soil pits except FB, which was outlier, the porosity and bulk density profile behavior also follows the measured values.

The determination of FB measurements as an outlier was in part visually and also based on the soil type analysis, and not based on any statistical test for outlier detection. For FB, the soil type is sandy clay loam, and therefore higher bulk density is expected (because it should have larger but fewer pore spaces). However, this behavior is not observed in our measurement (table S6). Compared to SGW-1, which is the same soil type and follows a higher bulk density and lower porosity, we considered FB-1 to be an outlier. Although it is also important to note that SGW-1 shows a relatively smaller porosity and higher bulk density due to the presence of gravel fraction. Nevertheless, the bulk density in FB-1 is smaller than expected within the range for sandy clay loam.

As for ICC, the model also cannot estimate the porosity well; however, the profile behavior still follows the measurements (figures 5(d)–(f)). Therefore, to summarize, we refer to FB site as an outlier, because not only the model shows error in estimating the absolute values, it also cannot show the profile behavior.

Comparing ${R^2}$ and the ${\text{RMSE}}$ values shows that the total OM, found from equation (5), does not show better performance in comparison to the OM profile that is based solely on equation (4), except at the IMN site (figures 6(a), (d), and table S4). Furthermore, one must note the difference between SOM and OM. The conventional LOI methods exclude particles with a dimension greater than 2 mm prior to the LOI process. Our method distinguishes between RB and SOM based on the sample particle dimensions as described in section 2.3. RB is not considered as soil; therefore, conventionally, it is being excluded from the LOI, resulting in the underestimation of the total OM content. However, this layer is in the subsurface and could substantially change subsurface thermal, hydraulic, and electromagnetic properties. The suggested exponential profile in equation (6) effectively captures the RB profile (figures 6(b) and (e)), while the S-shaped function shows similar utility in representing the SOM profile (figures 6(c) and (f)). Overall, the OM profile behavior performs well when it is reconstructed from the S-shape function, which is recommended for polarimetric SAR (PolSAR) and interferometric SAR (InSAR) retrievals. With the way that we distinguish between RB and SOM one can develop an empirical approach to finding RB from the total OM content. Such an approach has potential utility in estimating RB and the root layer from radar remote sensing but requires further research.

Zoom In Zoom Out Reset image size

Following soil texture characterization, the behavior of soil dielectric permittivity, electrical conductivity, soil moisture, and the corresponding SW fraction is of interest (figure 7). At the surface, due to faster evapotranspiration and lower soil moisture (figures 7(c) and (g)), and a highly porous RB layer, almost all of the samples exhibit a low dielectric constant (figures 7(a) and (e)). Among these sites, SGW-2, IMN-1, and FB-1 show a higher surface dielectric value because of higher moisture content (compare SGW-2 (blue), IMN-1 (blue), FB-1 (red) in figures 7(c) and (g)). As we go deeper into the active layer, the dielectric constant in all soil pits (except FB-1) reaches a maximum value (saturated organic layer) and then decreases as it reaches the deeper mineral layer. The transition depth varies based on the organic layer thickness (compare SGW-2 (black) and ICC-1 (blue) in figures 5(a), (d), 7(a), and 7(e)). The dielectric variation for the mineral layer is relatively constant due to saturation.

Zoom In Zoom Out Reset image size

The SGW-1 samples exhibit an unexpected low dielectric constant, and low soil moisture even though the mineral layer is fully saturated (figures 7(a), and (c), (d) SGW-1 (red)). The SGW-1 and SGW-2 soil samples are mostly sandy clay loam types (table S6). However, the measured porosity at the mineral horizon for the SGW-1 location is much lower than the expected porosity of sandy clay loam soils (compare SGW-1 (red) and SGW-2 (blue) in figure 5(c)). Furthermore, the mineral layers in SGW-1 and SGW-2 are fully saturated, but the dielectric constant and the soil moisture in SGW-1 are almost half of what they are in SGW-2 (figures 7(a) and (c)). The underlying reason for this behavior in SGW-1 is a substantially higher gravel fraction, leading to a lower porosity (figure S4).

The electrical conductivity, which determines the imaginary part of the soil dielectric constant and accounts for the lossy behavior of soil, is higher for the mineral layer compared to the surface organic layer. The abundance of cations in the mineral layer explains a higher electric conductivity (figure 7(b) SGW-2 (blue), HV-2 (magenta)) and, in turn, a higher loss for the mineral soil layer, which is consistent in all soil profiles (figures 7(b) and (f)). This is in agreement with the fact that the radar signal attenuation is greater in deeper soil layers, and radar sensitivity decreases accordingly (Chen et al 2019b).

While one can use various data-driven regression-based methods to develop a model to study the inter-relationships among soil physical properties, we focused on deriving a correlation-based PTF model that would be useful for further parametrization and prediction of Arctic soil hydraulic and dielectric permittivity behavior (Bakian-Dogaheh et al 2019, Chen et al 2019c). The developed PTF suggests that, except for sandy soils, the two essential soil physical properties ${\rho _{\text{b}}}{ }$and $\varphi$ can be found solely based on mineral soil texture and OM. This is of particular interest for developing a soil dielectric model for the permafrost active layer.

The detailed analysis and results findings suggest fewer independent variables as input parameters for soil dielectric modeling, which can avoid unnecessary modeling complexity and the need for including too many unknowns in the radar retrievals. We leave the topic of soil dielectric modeling to a future paper because of the need for extensive dielectric measurement for a wide range of organic soil and full soil moisture ranges from saturation to oven dry.

We showed that the PTF model was able to accurately find ${\rho _{\text{b}}}$ for the relatively high OM content soils in the near-surface layer. However, the model underestimates ${\rho _{{\text{b }}}}$ for the deeper mineral soil layers. This underestimation also propagates to the estimation of soil porosity depends on the mineral texture. However, the ${\text{OM}}\left( z \right)\,$ profile model can reconstruct a realistic soil physical model, in particular, by estimating the porosity profile $\varphi \left( z \right)$; where $\varphi \left( z \right)$ may have additional utility for estimating surface subsidence from permafrost degradation using InSAR remote sensing methods (Zwieback et al 2015, Hu et al 2018, Michaelides et al 2019, Chen et al 2020).

As depicted in the results, soil moisture and OM variations can explain the soil dielectric behavior throughout the active layer profile. These profile models are of particular interest for radar retrieval algorithms using both PolSAR and InSAR methods (figure 8). In a physics-based PolSAR retrieval, the subsurface parameters (${\bar X^t}$) serve as unknowns (for each radar pixel). At each retrieval iteration, an estimation of parameters reconstructs the subsurface profile (soil SW fraction and OM). Accordingly, these subsurface profiles translate into the dielectric profile by employing an organic soil dielectric model. Subsequently, the subsurface discretized into a multi-layer structure with N-layer and a depth interval of $\Delta z\, = \,\frac{{{z_{{\text{ALT}}}}}}{N}$. This discretized multi-layer dielectric structure along with the soil roughness height (h) feed into a multi-layer electromagnetic scattering model to find the simulated backscattered response of subsurface, where eventually these simulated values compared with the measured radar backscattered signal (figure 8).

Zoom In Zoom Out Reset image size

In summary the physics-based radar retrieval using PolSAR data requires three models: (a) A realistic subsurface soil properties profile model $\left( {{\text{SW}}\left( z \right),{\text{ OM}}\left( z \right)} \right)$ (b) A well-established organic soil dielectric model ($\tilde \epsilon \left( z \right)$), and (c) A multi-layer electromagnetic scattering model (figure 8).

The subsurface parameters can be categorized into three subgroups: (a) OM profile parameters $\bar X_{{\text{OM}}}^t\, = \,\left[ {{\text{O}}{{\text{M}}_{z0}},{\text{O}}{{\text{M}}_M},\beta ,m} \right]$, (b) soil saturation profile parameters $\bar X_{{\text{SW}}}^t\, = \,\left[ {{\text{S}}{{\text{W}}_{z0}},{z_{{\text{WT}}}}} \right]$, and (c) active layer parameters $\bar X_{\text{ALT}}^t = \left[ {h,\,{z_{\text{ALT}}}} \right]$, which are surface roughness height and ALT respectively. The distribution analysis of OM profile model parameters suggests that the deep OM content in the mineral layer ${\text{O}}{{\text{M}}_{\text{M}}}$ can be assumed to vary conservatively over the range from 5 to 7 (ignoring cryoturbation). Furthermore, the decay parameter $\beta$ can also be assumed to be nearly constant, at around $0.5$ (figure 9). In a PolSAR retrieval algorithm, these profile parameters serve as unknowns. The above analysis reduces the number of unknowns for ${\text{OM}}\left( z \right)$ from four to two, which is crucial for the PolSAR retrieval algorithm, given very few radar measurements (usually 2) are available.

Zoom In Zoom Out Reset image size

There are four primary sources of uncertainty in our measurements: (a) compressing highly organic porous samples, which results in over-estimating soil moisture and bulk density, (b) over-saturation while adding water to reach saturation, which may lead to overestimating porosity, (c) crushing/rubbing samples and passing excessive amount of OM through 2 mm sieve, which can lead to underestimating RB, and (d) the apparent dielectric permittivity measured by TEROS 12 sensors has an accuracy of ±1 for soil in the relative permittivity range of $1 - 40$ and 15% for soil with a relative permittivity range of 40–80.

Estimating the error and uncertainty of items (1)–(3) above will be qualitative because the error value is variable for different soil types, and since such an analysis will require a larger data set, we will leave this analysis to future work and a more comprehensive field campaign.

However, given that we have extracted multiple samples at each sampling depth, an estimation of the measurement precision (figure 10) can be calculated. Assuming the two adjacent samples (figure 2) should exhibit a similar soil physical property, we can find a distribution of correlation coefficient based on equation (11).

$$\begin{equation}a = \left[ {\rho _{\text{b}}^a,{ }{\phi ^a},{\text{R}}{{\text{B}}^a},{\text{O}}{{\text{M}}^a},{\text{SO}}{{\text{M}}^a},{S^a},{\text{S}}{{\text{i}}^a},{C^a},{\theta ^a},{\text{S}}{{\text{W}}^a}} \right]{ }\end{equation} \tag{ 10a }$$

$$\begin{equation}b = \left[ {\rho _{\text{b}}^b,{ }{\phi ^b},{\text{R}}{{\text{B}}^b},{\text{O}}{{\text{M}}^b},{\text{SO}}{{\text{M}}^b},{S^b},{\text{S}}{{\text{i}}^b},{C^b},{\theta ^b},{\text{S}}{{\text{W}}^b}} \right]{ }\end{equation} \tag{ 10b }$$

$$\begin{equation}{\boldsymbol{{\rho }}}\left( {a,b} \right) = \frac{{{{\mathop \sum \nolimits }}{{}}\left( {{a_i} - \bar a} \right)\left( {{b_i} - \bar b} \right)}}{{\sqrt {{{\mathop \sum \nolimits }}{{}}{{\left( {{a_i} - \bar a} \right)}^2}{{\left( {{b_i} - \bar b} \right)}^2}{ }} { }}}{ }{\text{.}}\end{equation} \tag{ 11 }$$

Zoom In Zoom Out Reset image size

In this paper, we presented a detailed study of soil properties in the active layer profile of Alaskan tundra. A fine-resolution dielectric (${\epsilon _r}^{\prime}$, and $\sigma$) profile characterization was provided, along with soil physical properties (${\text{OM}},\, {\text{RB}},\, {\text{SOM}},\,{\rho _\text{b}},\,\varphi ,\,\theta ,\, {\text{SW}}$) throughout the profile. The interrelationship between different soil textures (${\text{OM}},\,S,\,$and $C$) and soil physical properties (${\rho _{\text{b}}},{ }\varphi$) was used to develop a new PTF, which could be used in further modeling, especially to develop a soil organic model with a minimum number of parameters.

Furthermore, we showed that ${\text{OM}}\left( z \right)$ and ${\text{SW}}\left( z \right)$ follow independent S-shaped and quadratic profile distributions within the active layer. The parameters used to characterize these models were also studied. The resulting profile model can in turn be used to develop more realistic organic soil dielectric models for tundra, and the resulting dielectric profile can feed into a multilayer electromagnetic scattering forward model for predicting the radar backscattering coefficients (Tabatabaeenejad and Moghaddam 2006). Subsequently it improves the assessments and monitoring of permafrost active layer from airborne and spaceborne SARs.

Also, the study of the cryoturbation effect, in which the surface organic layer relocates to a deeper layer below the mineral horizon, has not been discussed in this work, although some of our observations were able to capture this phenomenon. This will be further studied in the future. The method used in this paper is a baseline study, and more extensive fieldwork and datasets could provide more insight into the fine resolution profile behavior of active layer soils.

The data that support the findings of this study are openly available at the following URL/DOI: https://doi.org/10.3334/ORNLDAAC/1759.

This work was supported by the NASA Earth and Space Science Fellowship (NESSF) Grant No. 80NSSC18K1410, and Arctic and Boreal Vulnerability Experiment (ABoVE), a NASA Terrestrial Ecology project, under Grant No. NNX17AC65A. Authors would like to greatly appreciate and thank the two anonymous reviewers whose valuable comments greatly improved this manuscript.