Mapping Surface Organic Soil Properties in Arctic Tundra Using C-Band SAR Data

Surface soil organic carbon (SOC) content is among the first-order controls on the rate and extent of Arctic permafrost thaw. There is a large discrepancy in current SOC estimates in Arctic tundra, where sparse measurements are unable to capture SOC complexity over the vast tundra region. Synthetic aperture radar (SAR) data are sensitive to surface vegetation, roughness, and moisture conditions, and may provide useful information on surface SOC properties. Here, we investigated the potential of multitemporal Sentinel-1 C-band SAR data for regional SOC mapping in the Arctic tundra through principal component analysis (PCA). Multiple in situ SOC datasets in the Alaska North Slope were assembled to generate a consistent surface (0–10 cm) SOC and bulk density dataset (n = 97). The radar VV backscatter shows a strong correlation with surface SOC, but the correlation varies greatly with surface snow, moisture, and freeze/thaw conditions. However, the first principal component (PC1) of radar backscatter time series from different years shows spatial consistency representing dominant and persistent surface backscatter behavior. The PC1 also shows a strong linear correlation with surface SOC concentration (R = 0.65, p<0.01), and an exponential relationship with bulk density (R = −0.65, p<0.01). The resulting predicted SOC maps show much lower soil bulk density and higher SOC concentration in the southern shrub tundra area than in the northern coastal region, consistent with in situ data. Our analysis shows that it is possible to separate the effects of different factors on the radar backscatter response using PCA and multitemporal SAR data, which may lead to more effective satellite-based methods for Arctic SOC mapping.


I. INTRODUCTION
T HE northern permafrost region contains a vast amount of soil organic carbon (SOC), which is potentially vulnerable to accelerated decomposition losses with continued warming and permafrost thawing [1], [2]. Highly organic soil, characteristic of tundra, can act as an effective insulator to protect underlying permafrost from surface warming due to its low bulk density and thermal diffusivity [3]. Surface SOC is particularly important in modulating active-layer thaw dynamics [4], [5]. Due to important controls of SOC on soil physical properties, land surface models or permafrost models rely on regional soil datasets for soil parameterization. Several soil parameterizations for organic soils have been developed for model applications in permafrost areas [6], [7], [8], [9]. Field and laboratory measurements show that soil bulk density has a strong influence on almost all soil physical properties [10], [11]. Soil bulk density decreases exponentially with increasing SOC concentration, with a well-established relationship between these two variables [11], [12]. Therefore, accurate estimates of soil bulk density or SOC concentration are essential for reliable model projections of the permafrost response to future climate trends.
Large discrepancies exist among current estimates of SOC properties in the northern high latitudes. A recent study reported more than a twofold difference in the SOC estimates from commonly used soil datasets including SoilGrids, the Harmonized World Soil Database (HWSD), and the Northern Circumpolar Soil Carbon Database, in the areas between 50°N and 80°N [13]. These datasets generally rely on geospatial analysis or statistical approaches to upscale in situ data to provide grid-cell mean estimates of SOC [2], [14]; however, the extremely sparse soil pedon measurements available in the northern permafrost region fail to capture the strong spatial variability of organic soils and contribute to the large discrepancy in regional SOC estimates [15]. In addition, regional SOC datasets targeting the permafrost region generally only provide estimates of the total SOC stock [1], [2], whereas soil bulk density or SOC concentration information has greater utility for modeling use. There is also a lack of high-resolution (<∼100 m) SOC datasets in the permafrost region, which are more suitable for representing the large characteristic tundra landscape variability.
Remote sensing data, particularly from optical remote sensors, have been widely used in regional soil mapping. Earlier studies have estimated SOC based on spectral reflectance mainly using regression methods. These methods may have limited accuracy in areas with significant vegetation cover, or areas This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ with large data gaps due to frequent cloud cover [16], [17]. Synthetic aperture radar (SAR) is sensitive to variations in soil moisture, vegetation structure, and surface roughness, which can provide useful information on soil properties, including SOC. SAR also has other advantages over optical remote sensing in the Arctic with its all-weather, day, and night imaging abilities. Therefore, SAR offers strong prospects for mapping diverse soil properties, including soil moisture [18], texture [19], [20], and chemistry [21]. Previous studies have successfully used X-and C-band backscatter to retrieve soil texture due to the strong control of soil texture on soil hydraulic conductivity and soil moisture dynamics [19], [20]. These studies have generally involved mineral soils, whereas similar applications involving organic soils are lacking, despite the large impact of SOC on soil hydraulic conductivity. However, a few studies have found close correlations between C-or L-band SAR backscatter and soil carbon concentration or bulk density in various ecosystems including agricultural lands, peatlands, and burned forest/shrub areas [22], [23]. The correlations were generally attributed to the sensitivity of SAR scattering processes to surface soil conditions. Moreover, lower-frequency SAR data generally show stronger sensitivity to SOC variations in more densely vegetated areas due to their greater penetration depth and stronger contribution from surface scattering processes [24]. More recent studies have explored the use of satellite-based SAR data in highresolution digital SOC mapping using statistical approaches or data mining methods, in addition to using multispectral or hyperspectral imageries [25], [26], [27], [28]. These studies show that including SAR images generally helps to improve model accuracy for SOC estimation than using optical images only, particularly for wetlands, although the importance of the SAR data may vary for different regions.
Relatively few studies have investigated SAR applications for regional soil mapping in the northern permafrost region, which may be due to the complexity of the radar scattering processes caused by variations in surface snow and vegetation cover, soil freeze/thaw transitions, and moisture changes [29], [30], [31], [32]. It is possible to separate the effects of different factors on the radar responses using multitemporal analysis based on statistical approaches [33]. Principal component analysis (PCA) provides a good way to reduce the redundancy in multiband or multitemporal remote sensing imagery to reveal the essential information contained in multidimensional data. After the principal component (PC) transformation, the first few PC in the new axes account for most of the variance in the original dataset and one or more factors may show more prominent effects on each PC. PCA has been widely used for change detection [34], land cover classification, and wetland mapping [35], [36], particularly based on optical remote sensing, whereas few applications have involved radar remote sensing for mapping soil properties and soil moisture variations [30], [33], [37].
SOC in the Arctic region shows strong spatial variability, accompanying similar heterogeneity in surface moisture, vegetation, and roughness conditions. C-band radar backscatter is sensitive to the above factors and may be effective in characterizing surface SOC patterns across vast tundra regions where soil measurements are lacking, but better soil information is needed  [38], [39]. Although many studies have used satellite SAR data in digital soil mapping mainly through data-driven methods, relatively few studies have investigated the sensitivity of satellite SAR data acquired in different seasons to variations in SOC properties, especially in the Arctic region. In this study, we investigated the sensitivity of multitemporal C-band radar backscatter time series to multiple factors including SOC and soil bulk density in Arctic tundra, and its potential for regional SOC mapping. There is generally a lack of in situ soil measurements in the Arctic for robust analysis. We, therefore, combined several in situ soil datasets to derive a comprehensive surface SOC and bulk density dataset in the Alaska North Slope as our study area. We used the PCA method to separate the effects of different factors on the C-band backscatter and evaluated the sensitivity of the first few PCs of backscatter to surface SOC properties and other confounding factors. The PC or backscatter showing the highest correlation with the SOC variables was used in the regional mapping.

A. Study Area and in situ Data
We selected a tundra subregion of ∼2°×4°located on the Alaska North Slope (≥ 68°N, Fig. 1) as our study area, which encompasses the Beaufort Coastal Plain and Brooks Range Foothills ecoregions [40]. The study area is within the core Arctic Boreal Vulnerability Experiment (ABoVE) domain [41] and has more in situ SOC data than other tundra regions in Alaska and elsewhere [42], [43], [44]. According to the National Land Cover Database (NLCD) [45], the Alaska North Slope is dominated by two major tundra types, including sedge/herbaceous and scrub/shrub tundra. Beaufort Coastal Plain in the north supports extensive lowland tundra plant communities, such as sedges, whereas the warmer Brooks Range Foothills in the south support tussock tundra, shrub tundra, and mixed tundra communities. Wetlands are also common in the northern coastal plains, due to poor drainage and waterlogged soils associated with underlying permafrost. Based on the NLCD map, the herbaceous tundra, shrub tundra, and wetlands account for 36.2%, 40.8%, and 4.9% of the study area, respectively, with the remaining area identified as barren land (∼9.3%), mainly located in the Mountains of the Brooks Range.
We used three in situ SOC datasets in our analysis (see Fig. 1 and Table I), including two published datasets [10], [42], [43]. The first dataset combined soil pedon data from the USDA-NRCS National Cooperative Soil Survey Soil Characterization Database and data from the University of Alaska Fairbanks northern soils research program. The dataset includes a total of 658 soil pedons for the Alaska region and is one of the most extensive soil C and N databases used in Alaska soil mapping. The dataset also includes soil classification, horizon designations, horizon boundary depths, 100°C oven-dried bulk density, C and N percentage by high-temperature combustion, and several other soil parameters [42]. There are about 124 pedons in the Arctic Alaska region, with an average sampling depth of 126 cm. We extracted data from sites with horizons spanning more than 5 cm in the top 10-cm soil layer, and with a horizon boundary no deeper than 15 cm below the surface. If multiple measurements were available for a single soil pedon, they were weighted based on the horizon depths. After processing, there were 54 sites with surface (0-10 cm) SOC and/or bulk density measurements in our study area.
The O'Connor dataset represents three watersheds in the Arctic Foothills region, with soil samples from more than 250 locations collected in 2017 and 2019. However, the soil organic matter (SOM) data were not reported at every sampling location [43]. The topography in this area is more variable with moderate to steep rolling hills. Soil stratification was determined for each site by measuring the vertical thickness of three main strata (acrotelm, catotelm, and mineral soil) typically found in organicrich or peat soils in our study area [see Fig. 1(c)]. For each site, an approximately 30 × 30 cm 2 section of soil was sampled for laboratory analysis, including soil bulk density and hydraulic properties. The SOM content was obtained from the samples by using loss on ignition (LOI) analysis and was converted to SOC concentration using a coefficient of 0.58 [46]. For this dataset, the SOC and bulk density were reported mainly for different soil strata. Therefore, the surface (0-10 cm) SOC data were calculated based on the weighted averages of SOC or bulk density of different soil strata spanning the surface layer. After processing, there were 38 sites with both SOC concentration and bulk density available.
In addition to the two above datasets, we also collected ∼112 soil samples from 5 sites along the Dalton Highway from field surveys conducted in 2018 [see Fig. 1(a)] [11], [44]. At each site, we took soil samples from 1 to 2 pits, with multiple soil samples taken along the soil active layer at ∼5 to 8 cm intervals for better characterization of the active-layer soil properties. Soil texture, including sand and clay fraction, root biomass, and soil bulk density were measured following standard protocols. The SOM content was determined using the LOI method, which heats the soil samples in a drying oven at 360°C for more than 2 h to combust the SOM. Mineral texture analysis was performed using the hydrometer method. Similar to the O'Connor dataset, we used the coefficient of 0.58 to convert the SOM measurements to SOC concentration.

B. SAR Data Processing and Analysis
Sentinel-1 satellite C-band (5.4 GHz) Ground Range Detected (GRD) scenes (path 94) from 2017 to 2018 were downloaded from the NASA SAR data archive at the Alaska Satellite Facility, and used for backscatter analysis in relation to surface SOC properties. For each acquisition, three scenes (frames 220, 225, and 230) were processed and mosaicked to generate a full coverage of the study area. There are 24 and 21 images in total with VV/VH polarized channel data acquired in years 2017 and 2018, respectively. In 2017, one image from late October was dropped due to a missing frame (220) in the southern region. In 2018, there were no GRD data released for the granule during April and May. Other granules were not used to fill this gap due to the different viewing geometries of the images. The GRD data processing includes thermal and radiometric calibration, speckle noise filtering, geolocation, and terrain calibration using the 5-m Alaska IFSAR DEM data. The output images were multilooked and resampled to 100-m resolution. This intermediate resolution was chosen to reduce computational demand, and also due to the relatively coarse resolution of the soil datasets (SoilGrids, 250 m) used in this study. The GRD data were processed using the SNAP-Python package v7.0.
Incidence angle has a relatively large impact on the radar backscatter response. Previous studies have shown that a reasonable linear fit between backscatter (at dB scale) and incidence angle can be obtained for incidence angles between 30°and 60° [47], [48]. In the study area, the local incidence angle of the SAR images mainly ranges from 36°to 42°with a mean value of ∼39°. We used a simple equation to normalize the local backscatter distribution from any incidence angle (θ) to the reference incidence angle where β is an empirical parameter, and falls mostly within the range of −0.1 to −0.2 dB/°for the tundra area [47]. In this study, we used the value of −0.13 dB/°following [48]. This value has a relatively small impact on the following PCA and correlation analysis since there is limited variation in the local incidence angle in most of the study area except for the southern mountain areas of the Brooks Range. However, we masked out the above areas with complex mountain topography, identified as "barren land" in the land cover map prior to the PCA. We applied PCA to the backscatter time series to analyze the potential relations between the PCs and SOC properties. Preliminary analysis showed an overall better correlation between the regional SOC data (SoilGrids) and backscatter at VV polarization than with VH polarization, particularly during winter. Therefore, only VV backscatter was used in the following analysis. The backscatter data show different spatial variations during both summer and winter in 2017 and 2018. The year 2017 had a drier spring season with only a few millimeters of total precipitation in June and a later autumn freeze onset, than the year 2018. Therefore, we applied PCA to the backscatter time series in 2017 and 2018 separately. The regional mean of the backscatter was first removed from each image, and the difference images were used as inputs to the PCA. Mountain areas in the southeast portion of the study area have lower incidence angles and much higher backscatter than other regions due to their greater topographic complexity, whereas large uncertainties may also exist in the terrain correction of the radar imagery in these areas. These characteristics of complex terrain could dominate the initial PCs of the larger study region and introduce uncertainties in determining regional soil properties. Therefore, we masked out mountain and ocean areas within the domain prior to the PC analysis to minimize the influence of these areas on the estimated SAR soil properties.

C. SOC Estimation Using Radar Backscatter and PCs
We mainly used correlation analysis to assess the relations between the SOC data and radar backscatter imagery (including the PCs). Linear or exponential relationships were tested between the in situ SOC data (including SOC concentration and bulk density) and extracted radar backscatter data (including the PCs). The established empirical relationship was then used to predict regional bulk density based on the radar backscatter imagery. PC1 was used instead of the raw backscatter data due to the large variability in backscatter time series. We used the following cross-validation method to evaluate the robustness of the empirical equations derived from the above analysis. We first randomly selected 80% of grouped PC1 and SOC data (bulk density or SOC concentration) for the model training, whereas the remaining 20% of samples were used for validation. A linear fit was used for estimating SOC concentration, and an exponential fit model was used for the bulk density estimates. Due to a limited set of observations (∼90 independent sites), we repeated the above experiment 10 times and calculated the mean root-mean-square error (RMSE) and bias from the results to evaluate the performance of the empirical equations.
Our preliminary analysis showed large variability in the relations between SOC and backscatter data for individual grid cells. Inconsistency in the methods for SOC estimation of the surface layer (0-10 cm) among the three in situ datasets, particularly for the O'Connor dataset, may contribute to the variance. Therefore, we used the following equation to estimate the SOC concentration (g/g): where ρ b (g/cm 3 ) is the predicted soil bulk density, and ρ b,min (g/cm 3 ) is the soil bulk density associated with the mineral texture fraction, which is calculated based on the sand and clay fraction from SoilGrids data. The above exponential equation was derived using the in situ soil samples from the North Slope region [11], and shows stronger correspondence (R 2 = 0.91) and lower RMSE values particularly for soil samples with higher carbon content (SOM>0.4 g/g). An exponential relationship between SOC and soil bulk density has also been reported in previous studies [12].
In addition to the above in situ SOC datasets, we also used SoilGrids (version 1 [49] and version 2 [50]) in the regional analysis. SoilGrids was chosen mainly because it shows relatively high accuracy [14] and has a much finer spatial resolution (∼250 m) than other regional and global soil inventory records, which are generally coarser than 1 km. SoilGrids also provides estimates for both SOC concentration and bulk density at multiple soil depths, whereas other local datasets tailored for the permafrost region such as NCSCD (Northern Circumpolar Soil Carbon Database) and a regional upscaling dataset [2] only provide estimates of the bulk SOC content (kgC/m 2 ) at coarser spatial and/or vertical resolution. However, the two SoilGrids data versions show different SOC patterns and levels over the North Slope region (not shown). In the study area, a visual comparison of these data also indicates a large discrepancy in the bulk density data at the 0-5 cm interval, whereas some areas show unrealistically high values in both SOC concentration and bulk density in the SoilGrids version 2 data (see Fig. S1). The SoilGrids version 1 estimates of surface SOC and bulk density are more consistent with the in situ SOC data for the study region (SOC: R = 0.78, bias = −4.7%, and RMSE = 0.12 g/g; bulk density: R = 0.79, bias = 0.19, and RMSE = 0.22 g/m 3 ) than the version 2 data (SOC: R = 0.53, bias = −0.81%, and RMSE = 0.13 g/g; bulk density: R = 0.59, bias = 0.14, and RMSE = 0.23 g/m 3 ) (see Fig. S2). Therefore, we mainly used the SoilGrids version 1 data in our regional analysis.

D. Other Ancillary Datasets
Other ancillary data used in the radar data analysis included the 2016 30-m NLCD land cover data, 5-m Alaska IFSAR DEM data, and Daymet surface meteorology data [51]. All datasets were resampled to 100-m resolution, consistent with the processed radar data. Mean elevation and slope were generated from the Alaska DEM dataset and used for the analysis. Daymet surface meteorology data were used to evaluate the backscatter response to changing surface freeze/thaw and moisture conditions. We calculated the accumulated precipitation in the prior 12-day period for each satellite acquisition date. Maximum and minimum air temperature data were used to define the surface freeze/thaw onset, which was defined as the date when both daily air temperature records and their 7-day moving averages drop below (or rise above) 0°C. The precipitation and freeze/thaw onset estimates from the Daymet data were evaluated using additional weather station data from three SNOTEL sites (SNOwpack TELemetry, http://www.wcc.nrs.usda.gov) located in the study area [see Fig. 1(a)]. Overall, Daymet provides a reasonable estimate of the precipitation pattern, which generally increases from north to south within the study region (see Fig. S3).

A. Principal Component Analysis of C-Band VV Backscatter
Different land cover types show distinct C-band VV backscatter characteristics (see Fig. 2). The vegetated areas generally show higher backscatter (∼3 to 4 dB) during summer than during winter, whereas barren land has much higher backscatter overall throughout the year than other areas, with reduced differences between summer and winter. The summer backscatter is also more variable than the winter backscatter. Among all vegetated areas, the scrub/shrub tundra shows the highest backscatter throughout the year, and wetland areas have the lowest backscatter, particularly during the winter. In 2017, a large decrease (∼2 dB) in radar backscatter was observed during the spring thaw (around DOY 143), due to thaw-induced surface inundation. We were unable to determine the backscatter response to the spring thaw in 2018 due to missing SAR data during the seasonal transition in April and May. A gradual decrease in the radar backscatter was observed during June and July in 2017 due to dry conditions from very low precipitation, whereas the backscatter during the thaw season of 2018 is more variable. Increases in backscatter were generally observed in response to wetter conditions from high antecedent precipitation in 2018, but not in 2017. The backscatter decreases after the freeze onset, and the spatial pattern of the winter backscatter is different between 2017 and 2018, likely associated with variations in winter snowpack distribution.
Although the winter and summer backscatter of years 2017 and 2018 show quite different spatial patterns, the PCA produces similar results for both years. The first three PCs account for 56%, 19%, and 7% of the total input variance in 2017 and 48%, 26%, and 9% of the total variance in 2018. The first PC (PC1) generated from the backscatter time series in 2017 and 2018 shows a very similar spatial pattern between years [see Fig. 3  . From the PC1 image, we can also see that shrub tundra areas generally have positive values, whereas other tundra areas have mostly negative values. With increasing PC1 values, the surface soil bulk density generally decreases, wheras the SOC concentration generally increases, indicated by the SoilGrids data. It should be noted that the SoilGrids version 1 data generally show higher SOC concentration and lower bulk density at 5 cm depth than at the surface (0 cm). The SoilGrids estimates at 5 cm depth also show an overall higher correlation and lower bias compared with in situ soil measurements for the Alaskan domain, particularly for soils with high SOC content (not shown). Therefore, we used the SoilGrids estimates at 5 cm depth in the regional comparison with our SOC estimates for the surface layer, which was defined as the interval of 0-10 cm, to be consistent with the assembled in situ dataset.
The PC2 and PC3 can be explained by the variations in regional precipitation, spring thaw onset, and winter snowpack (see Fig. S4 and Fig. 4). In 2017, the images acquired on DOY 227 and DOY 251 have the largest contribution to PC2, whereas, in 2018, the acquisitions for DOY 186 and DOY 210 have the largest contribution to PC2, based on the PCA loadings. These dates generally correspond to large antecedent precipitation (from the prior 12-day period) according to the Daymet and SNOTEL precipitation data, particularly in 2017. In 2018, the in situ precipitation data at the three SNOTEL sites show large spatial variability, compared with the more consistent regional drying/wetting patterns in year 2017; this may contribute to temporal differences between the PC2 loadings and the precipitation data in 2018 (see Fig. S4). The regional variability in the spring thaw onset also likely contributes to large changes in the PC2 loadings in June. As for PC3, there are much larger differences in the backscatter between areas with positive and negative PC3 values during winter (∼3 dB) than during summer (<1 dB). Areas with higher PC3 values seem to correspond with higher snowpack according to the Daymet SWE data (see Fig. S5). In addition, the loadings have larger negative values (−0.2 to −0.3) for PC3 during winter than the loadings of the first two PCs, which indicates a higher contribution of winter acquisitions to PC3. However, the Daymet data are unable to capture the strong variability in winter snowpack conditions due to the sparse regional in situ weather station network, which prevents a quantitative comparison.

B. Regional SOC Mapping
The radar backscatter shows a significant (p<0.05) correlation with the in situ SOC data. However, the correlation varies greatly (R = ∼0.3-0.64) for different seasons in different years (see Fig. S6). Changes in winter snowpack, summer moisture and vegetation conditions may all contribute to variations in the radar backscatter, and thus affect the correlation between the backscatter and surface SOC properties. In 2017, the radar backscatter in late winter shows the highest correlation with the in situ SOC data; whereas the radar backscatter in the early thaw season of 2018 shows the highest correlation with in situ SOC data. Backscatter acquired on DOY 131 (2017) and DOY 210 (2018) show the strongest correlations with surface SOC properties, including SOC concentration (R = 0.64, p<0.01) and soil bulk density (R = −0.62, p<0.01) among all images (see Fig. 5). The PCA can separate the effects of different factors, and the transformed PCs may provide more reliable estimates of surface SOC. As shown above, PC1 derived from different backscatter time series show a consistent spatial pattern and may be used to predict surface SOC properties. The correlations between PC1 and surface SOC properties are comparable or slightly better than correlations between the in situ SOC data and backscatter acquired during the winter in 2017 and the early thaw season in 2018 [see Fig. 5(b) and (d)]. A favorable linear correlation was found between PC1 and SOC concentration (R = 0.65, p<0.01), whereas an exponential equation provides a better fit for the relation between PC1 and soil bulk density (R = −0.65, p<0.01). The cross-validation analysis obtained similar exponential or linear relationships explaining the soil bulk density and SOC concentration patterns (see Table S1). The cross-validation results also indicate that the exponential fit between soil bulk density and PC1 works well, with mean RMSE of 0.24±0.09 g/cm 3 and 0.20±0.05 g/cm 3 for years 2017 and 2018, respectively (see Fig. 6). However, there is more scatter in the model predicted and in situ SOC concentrations. Therefore, (2) was used for the SOC estimation as described in Section II-C. Fig. 7 shows the predicted regional SOC maps and uncertainty estimates using PC1 of the radar backscatter data. The soil bulk density was estimated using the empirical relationship shown in Fig. 5(d), whereas the SOC concentration was derived from the estimated bulk density and (2) due to a larger scatter between the in situ SOC data and PC1 as shown in Fig. 5(b). The uncertainty estimates were defined as the 95% confidence range. The southern areas of the domain dominated by scrub/shrub tundra generally show higher SOC concentration and lower soil bulk density compared with other tundra areas. Larger uncertainties were mostly observed in the northern coastal region. The estimated SOC and bulk density parameters show quite different spatial  patterns from SoilGrids (see Fig. S1). The scrub/shrub tundra shows much lower bulk density and higher SOC concentration than the two SoilGrids datasets, with more consistency between the in situ data and the estimated values from this study (see Fig. 8). For herbaceous tundra, our results show similar SOC but lower soil bulk density than the two SoilGrids datasets, and large variability against the in situ data. There is also a large difference between the in situ data and the other three soil datasets in the wetland area. However, there are very limited observations and greater uncertainty associated with this land cover type (3 sites), compared with the other two tundra types (scrub/shrub: 54 sites; herbaceous: 31 sites).

IV. DISCUSSION
Our results show strong correlations between Sentinel-1 Cband VV backscatter and SOC properties affecting active layer and permafrost conditions in tundra. These results indicate the potential utility for satellite SAR-based SOC mapping across the pan-Arctic. However, the relationship between radar backscatter and soil properties can vary greatly for different seasons and years, depending on heterogeneous surface conditions including winter snowpack and summer moisture levels indicated from this study, and vegetation growth indicated from previous studies [32]. A strong correlation between C-band winter backscatter and SOC storage has previously been reported for Arctic tundra, but not for summer backscatter [31]. This is likely because the backscatter under frozen conditions can better represent surface roughness and vegetation structure, which are interactive with SOC and other soil properties. Some prior studies have also found that compared with summer backscatter, winter C-band backscatter is more useful for land cover and surface characterization due to enhanced sensitivity to surface roughness and vegetation structure [29], [47]. However, the presence of seasonal snow cover may increase volume scattering and thus affect the C-band backscatter response [52], limiting the use of winter backscatter data for surface characterization. For example, we found a much lower correlation between surface SOC properties and winter backscatter in 2018; that was different from the previous year (2017) and was likely due to varying snow conditions between years (see Figs. S5 and S6). We also found that the VV backscatter data show an overall better correlation with SOC properties than the VH backscatter, especially during winter. This is likely because volume scattering contributes more to the VH backscatter; therefore, the VH backscatter is more sensitive to surface vegetation than soil conditions in the tundra area [23].
How the radar backscatter responds to precipitation or soil wetness depends on the land surface physical properties. Soil texture and SOC content mainly determine soil hydraulic properties; therefore, multitemporal radar backscatter can potentially be used to map soil physical properties including SOC. A few studies have tried to use multitemporal radar data to retrieve soil texture in semiarid regions [19], [20]. Another study also demonstrated the usefulness of SAR data for soil drainage mapping in an agricultural field using a maximum likelihood classifier, but with less accuracy than using hyperspectral reflectance data [53]. However, comparatively few studies have used SAR data to retrieve SOC or soil moisture in Arctic tundra. The relatively long Sentinel-1 revisit period and strong C-band radar backscatter sensitivity to vegetation may limit such applications in tundra [23]. The heterogenous microtopography and highly organic soils typical of Arctic tundra impart unique microwave scattering properties and the polarimetric response can be highly random, which also limit the utility of SAR data for mapping soil properties in Arctic tundra [18]. Our study did not show a consistent C-band backscatter response to antecedent precipitation. The C-band radar backscatter sensitivity to soil moisture is limited to a very shallow surface layer (<5 cm depth). Our recent study [9] documented surface layer soil drydown characteristics in the study area, and found a more rapid drydown (∼5 to 10 days) in the southern organic-rich soils; that was generally shorter than the satellite revisit period (i.e., 12 days). The use of longer wavelength (e.g., L-band) observations and data fusion methods able to exploit complementary information from active and passive sensors could enable more effective soil retrievals and a better understanding of the microwave sensitivity to different soil properties [23], [54].
The shallow penetration depth of C-band SAR and inconsistent data sampling among the in situ datasets may add considerable uncertainties to our SOC estimates. The direct sensitivity of C-band radar is largely confined to the top few centimeters of the surface, and may reveal less information about deeper soil layers, especially in areas with denser tundra shrub covers. Previous studies have shown stronger L-band microwave sensitivity to SOC and soil bulk density in forest and peatland areas than Cband, consistent with enhanced soil sensitivity at lower frequencies [21], [22]. Polarimetric decomposition using quad or dualpolarized SAR data can effectively separate the contributions from volume and surface scattering and may provide additional insights on whether the shallow C-band penetration depth limits its use of representing surface conditions in Arctic tundra [23].
Our study did not consider temporal changes in SOC properties. However, our study period was relatively short, during which limited SOC changes are expected to occur in the undisturbed tundra [55]. On the other hand, we assembled three different in situ SOC datasets, and the associated differences in the soil sampling process may add large uncertainties to the assembled surface (0-10 cm) SOC dataset. In particular, the O'Connor dataset only provides measurements of SOC concentration or bulk density for different soil strata (e.g., acrotelm and catotelm), whose depths vary greatly depending on microtopography and vegetation types within the study area [10].
A number of studies have used various statistical approaches or data mining methods with inputs from multisensor remote sensing data to map regional SOC patterns for peatlands and other landscapes [25], [26], [27], [28]. Multitemporal SAR data have been found to provide useful information for such applications [25], [26]. One study shows that including SAR data at longer wavelengths such as L-band helps distinguish different peatland classes from nonpeatlands, but with less importance than C-band [25]. The different studies also report different conclusions on the relative importance of optical remote sensing versus SAR data for regional SOC characterization [25], [26], [27], which likely depends on the regional differences in topography, vegetation, and surface wetness conditions represented from the various studies. The strong sensitivity of radar backscatter to multiple factors including land cover, surface snow, and moisture changes affected by heterogeneous freeze/thaw and precipitation events may add additional challenges for regional mapping. The PCA can separate these different factors and extract essential information from the multitemporal radar data, which may provide useful inputs to other statistical approaches or data mining methods such as Random Forest [30]. For example, in our study area, the PC1 seems to show dominant and persistent backscatter features of the land surface, closely linked to vegetation type and SOC characteristics (see Fig. 3). Previous studies have drawn similar conclusions that PC1, together with the other first few PCs, can effectively distinguish surface features pertaining to various vegetation, topography, and roughness conditions [30], [33], [35], [36]. Therefore, a combination of PCA and other data mining methods may prove to be a useful approach for satellite SAR-based SOC mapping in tundra [34].

V. CONCLUSION
In this article, we investigated the potential of multitemporal C-band (5.4 GHz) SAR data for regional SOC mapping in the Arctic tundra. We combined multiple in situ soil datasets to generate consistent SOC and bulk density estimates for the surface (0-10 cm depth) soil layer over an Alaska North Slope subregion. Our analysis shows that both winter and summer Cband backscatter show a strong correlation (R>0.6, n = 97, and p<0.01) with in situ surface SOC data. However, the correlation between the backscatter and SOC data depends on additional radar sensitivity to surface snow, moisture, and freeze/thaw conditions, which can vary greatly for different seasons and years. Through PCA using multitemporal Sentinel-1 SAR data, we separated the effects of these different factors on the radar response. The first PC showed a consistent spatial pattern across different years, and represents persistent backscatter behavior of land surface, which is closely correlated with SOC and bulk density in surface soils. The predicted SOC maps show much lower soil bulk density and higher SOC concentration in the southern shrub tundra area than in the northern coastal region, and are more consistent with the in situ data, compared with the SoilGrids products. Our study shows the potential of satellite C-band SAR data for regional SOC mapping in the Arctic. Future work will further investigate the sensitivity of polarimetric SAR data at different frequencies to surface soil conditions. We also plan to combine PCA with other data mining methods to generate high-resolution and high-quality SOC maps in the Arctic using multitemporal and multifrequency SAR data.