Fine-Scale SAR Soil Moisture Estimation in the Subarctic Tundra

In the subarctic tundra, soil moisture information can benefit permafrost monitoring and ecological studies, but fine-scale remote-sensing approaches are lacking. We explore the suitability of C-band SAR, paying attention to two challenges soil moisture retrieval faces. First, the microtopography and the heterogeneous organic soils impart unique microwave scattering properties, even in absence of noteworthy shrub cover. Empirically, we find the polarimetric response is highly random (entropies >0.7). The randomness limits the applicability of purely polarimetric approaches to soil moisture estimation, as it causes a tailor-made decomposition to break down. For comparison, the L-band scattering response is more surfacelike, also in terms of its angular characteristics. The second challenge concerns the large spatial but small temporal variability of soil moisture observed at our site. Accordingly, the Radarsat-2 C-band backscatter has a limited dynamic range (~2 dB). However, contrary to polarimetric indicators, it shows a clear surface soil moisture signal. To account for the small dynamic range while retaining a 100-m spatial resolution, we embed an empirical time-series model in a Bayesian framework. This framework adaptively pools information from neighboring grid cells, thus increasing the precision. The retrieved soil moisture index achieves correlations of 0.3–0.5 with in situ data at 5 cm depth and, upon calibration, root-mean-square errors of <0.04 m3m−3. As this approach is applicable to Sentinel-1 data, it can potentially provide frequent soil moisture estimates across large regions. In the long term, L-band data hold greater promise for operational retrievals.


I. INTRODUCTION
Subarctic tundra ecosystems are responding rapidly to climate change, as evidenced by widespread changes in their vegetation cover and permafrost conditions [1], [2]. These ecosystems are characterized by organic soils that host mosses, lichens, graminoids and forbs of limited biomass as well as, increasingly, higher-biomass shrubs. The permafrost has widely been observed to be warming, often associated with local permafrost degradation [2], thus threatening infrastructure. Soil moisture plays a central role in the changing ecology and permafrost conditions [3], [1], [4]. It exerts an important control on plant growth and a wide range of biogeochemical processes, and it alters the soil thermal dynamics through its control on the surface energy balance and the soil thermal properties.
Despite the clear need for frequent, spatially extended soil moisture information, we currently lack reliable operational satellite remote sensing solutions [5]. The coarse-scale passive microwave SMAP product is widely considered to be the most accurate global product, but it has been found to provide poor estimates in the subarctic tundra [5]. Alternatively, finescale SAR approaches seem promising, as the large spatial heterogeneity in land surface properties at the 100-200 m scale would make them attractive for a wide range of applications. However, the applicability of SAR data for soil moisture inversion in the subarctic tundra remains largely unknown. This also applies more generally to high-latitude regions, although pioneering studies over the high arctic tundra and boreal peatlands and fire scars give reason for hope [6], [7], [8], [9], [10], and recent efforts such as NASA's ABoVE and the German-Canadian PermaSAR (DLR) campaigns [11], [12] are paving the way for further algorithmic advances in highlatitude regions.
Manuscript received x Corresponding author: S. Zwieback (email: zwieback@uoguelph.ca). SAR remote sensing of soil moisture faces numerous potential challenges in the subarctic tundra, in terms of the hydrology and microwave scattering properties. The hydrological behavior is complex, owing partially to the meter-scale microtopographical features. These can take on many forms; in the subarctic, hummocks of ∼ 20 cm height are widespread [13]. They modify the soil moisture distribution not only due to their relief but also because they impart a sizable small-scale variability in soil properties. In general, the near-surface soils are largely organic in nature, characterized by large porosities, but their porosities and water retention characteristics vary greatly [13]. The implications for SAR retrievals are that soil moisture is highly heterogeneous on the sub-resolution scale [14], thus posing a challenge to validating the retrievals. Also, the temporal variability can be small by comparison, which again makes time-series analyses of SAR data difficult [5].
In terms of the complex microwave scattering, the most obvious challenge are the shrubs, as above-ground vegetation influences the observed signals in complex ways and generally reduces the sensitivity to soil moisture variations [15], [16], [17]. But above-ground biomass remains very limited in vast areas, and even the open tundra presents challenges of its own. First, the microtopography likely makes the surface appear very rough in terms of its polarimetric and angular backscattering behavior [18]. As the size and spacing of the hummocks vary, the variability needs to be accounted for when interpreting the SAR measurements. Second, the organic soils also complicate soil moisture retrieval. When these soils are dry, microwaves can penetrate easily into the subsurface, and subsurface scattering can contribute appreciably to the received signal [7], [19]. In the subarctic tundra subsurface scattering may, for instance, originate from the interface between the organic and the underlying mineral soil [13]. The microtopography and the organic soils likely impart unique scattering characteristics (polarimetry, angular dependence), but these have not been quantified, and dedicated scattering models are lacking.
In light of these complexities, the most expedient approach to soil moisture estimation remains an open question. Three general approaches are commonly applied, but none has been tested in the subarctic tundra. First, model-based backscatter approaches furnish soil moisture estimates by inverting a model of the backscatter magnitude. Existing models vary greatly in complexity, ranging from computer-intensive electromagnetic simulations to simple semi-empirical models [20], [21], [22]. However, there are no dedicated scattering models available for the hummocky organic soils in the subarctic tundra, and it is not clear whether existing ones can replicate the angular dependence and other scattering characteristics. Second, simple time-series approaches have been applied successfully to a wide range of land covers, even those with appreciable vegetation [23], [17], [24]. However, their flexibility comes at a price, namely that the estimates are not calibrated, i.e. they are only relative indicators of moisture content. Third, there are purely polarimetric methods that do not rely on the backscatter magnitude, chief amongst them being decomposition methods [16], [25]. They try to isolate the surface response whose observed polarimetric scattering mechanism can then be converted into an absolute soil moisture estimate. Numerous such approaches have been applied with some success in agricultural regions [16], [26], but most approaches cannot account for the variable microtopography. In light of their simplicity, polarimetric methods and simple empirical methods seem to constitute a promising first step, which could ultimately lead to bespoke backscatter-based models.
To address these open questions, and more generally provide a first overview of the limitations and opportunities of soil moisture estimation in the subarctic tundra, we pursue 3 closely linked objectives 1) to characterize key polarimetric properties of the hummocky tundra, and to develop a tailor-made Micro-Topography and Vegetation (MTV) decomposition to separate the surface (with variable microtopography) from the volume scattering contribution. 2) to characterize the backscatter's angular dependence, and to determine the applicability of semi-empirical backscatter models and of a linear empirical backscatter model. 3) to assess simple soil moisture inversion approaches over the open tundra: we study two polarimetric approaches and one magnitude-based approach based on the empirical time-series model.
We analyze two years (2014, 2016) of fully polarimetric C-band Radarsat-2 data at a well instrumented study site in the Northwest Territories, Canada. While C-band is rarely considered optimal for soil moisture estimation, such data are the most widely available, especially since Sentinel-1 started to frequently image northern Canada in late 2016. Instead of C-band, lower frequency data are generally considered preferable, and NISAR and TanDEM-L will provide frequent fine-scale data at L-band. To provide a first insight into the scattering properties at L-band, we also studied one airborne UAVSAR acquisition (objectives 1 and 2).
We first provide an overview of the study site in Sec. II, highlighting the large observed subresolution spatial variability in soil moisture, which contrasts with a low temporal variability. In Sec. III, we present the Radarsat-2 C-band and UAVSAR L-band data. We then address each objective in turn in Sec. IV-VI, describing first the methods and subsequently the results. Our polarimetric (objective 1) and angular (objective 2) analyses inform our investigation of Cband soil moisture estimation (objective 3). We focus on a simple linear time series model of the backscatter magnitude. In light of the low dynamic range of surface soil moisture observed at our site, we modify the standard pixel-level timeseries approach by embedding it in a Bayesian framework that adaptively pools information across adjacent pixels, thus reducing the noise level. We assess the soil moisture retrievals at two instrumented plots in the open tundra, for both the VV backscatter (available from Sentinel-1) and the surface backscatter extracted using the MTV decomposition from objective 1.  (Fig. 1a,b). Three main types of land cover can be distinguished [1]. The open tundra's sparse vegetation cover is dominated by lichen, bryophytes and graminoids (Fig. 1c) [13], [5]. Dwarf shrub tundra additionally supports erect dwarf shrubs of 50 cm height (Fig. 1d). Finally, tall shrubs (40 -200 cm) occur along water channels, in disturbed areas and, increasingly, in patches on hillslopes (Fig. 1d). The rolling morainal landscape is underlain by continuous permafrost [2]. During the short summer season, the soil thaws to a depth of 0.4-1.3 m, depending on, amongst other things, the soil profiles.

II. STUDY SITE AND
The soils and the microtopography reflect the area's glacial and periglacial history. The soils are organic cryosoils. The silty clay mineral soils are overlain by organic materials of up to 50 cm depth [13]. The organic materials are characterized by a large porosity (0.60 -0.95). Their nature and depth vary on small scales, as they are closely associated with the microtopography. This microtopography takes the form of hummocks, which rise up to 30 cm above the surrounding interhummocks (Fig. 1c). The hummocks, which are around 50-100 cm in diameter, consist of a mineral soil core overlain by a thin layer of organic soil. The layer is about 15 cm deep in Fig. 2a but can vary greatly in depth on short spatial scales. The interhummocks, by contrast, are characterized by thicker organic layers, mainly comprising living mosses and peat. Apart from the hummocky tundra, there are isolated patches of degraded ice-wedge polygons, where such hummocks are largely absent.
Soil moisture varies on a range of scales, as its distribution is shaped by the microtopography, topographic position and vegetation cover. Ideally, SAR remote sensing should capture the hillslope-scale variability (100-300 m). However, the subresolution variability has to be considered in the following, as it was observed to be large. Large spatial variability in tundra landscapes is associated with the microtopography and variations in thaw depths [14]. To quantify it, we conducted a manual sampling campaign in September 2017 in four plots with variable shrub cover (OT-MM, DS-1, DS-2, TS-2, see Fig.  1b). Each plot comprised two parallel transects of 20 samples each, spaced 1 m apart. The measurements were made using a Stevens Hydra Probe II capacitance probe (5 cm prongs, inserted vertically). Once calibrated to the soils at TVC [27], this probe was previously found to achieve an accuracy <0.05 m 3 m −3 . Our measurements shown in Fig. 2b) revealed a large plot-level interquartile range of 0.08-0.29 in the four plots, corresponding to up to 70% of the median soil moisture. The temporal variability in surface soil moisture is subdued compared to the spatial variability. To capture the temporal dynamics of soil moisture, we employed permanent soil moisture probes during two summer seasons (2014, 2016; coincident with the Radarsat-2 data). Each plot was located in hummocky open tundra and was instrumented with 3-6 near-surface probes distributed across hummocks and interhummocks. The probe locations were chosen to represent the spatial variability within the 100-200 m plot, but owing to the limited number of probes per plot, representativeness errors remain. The probes we used were Stevens Hydra Probe II sensors at 5 cm depth. A location-specific calibration was used to derive soil moisture (spatial scale of <10 cm) from the measured dielectric constant [5]. We instrumented two plots, OT-MM and OT-UP, where OT indicates the tundra vegetation (open tundra with very little above-ground biomass). OT-MM was located on a plateau close to the camp's main meteorological station, and OT-UP was located on the Upper Plateau, 10 km to the southwest of the main study area (Fig.  1). To obtain a representative volumetric surface soil moisture value, we aggregated the probe-level measurements within each plot. As shown in Fig. 2c, the individual probe-level measurements indicated a limited temporal variability. At any point in time, the spatial variability of the five probes is considerably larger than the change in the spatial mean over the entire summer season. For all plots, the aggregated plot-level observations exhibited a lower temporal variability than the spatial variability we observed in our dense transects. Similar findings have been reported from Alaska, where surface soil moisture was highly variable but also showed a clear but subdued response to precipitation events [14], [28].

III. REMOTE SENSING DATA
Our analyses focused on a dense time series of fully polarimetric Radarsat-2 data [29], taken in the stripmap Fine Quad mode with a native slant range and azimuth resolution of ∼ 8 m. The data were acquired during the summers of 2014 and 2016 (01 July to 31 August). During these time periods, there was no snow present, and the soils were thawed to a depth of at least 25 cm. We thus considered the soil as a thawed medium for the purpose of C-band remote sensing. The dense temporal sampling facilitated the interpretation of the data (every 3-5 days: 13 acquisitions in 2014, 16 in 2016). However, not all of these acquisitions covered the OT-UP plot, due to the limited swath width. The incidence angles (IA, θ) varied from low (≈ 22 • ) to high (≈ 45 • ), with an approximately equal split into low and high IA acquisitions.
We applied standard processing steps to convert the disseminated Single Look Complex data into multilooked, calibrated, georeferenced images of the radar backscatter [30], [7]. The multilooking was implemented in two ways: a) a mediumresolution (50 m at low IA, 40 looks) and b) a low-resolution (100 m at low IA, 160 looks) version. The advantage of the lower resolution is a doubling of the precision with which the backscatter from a distributed target can be measured (0.34 vs 0.69 dB) [31]. The resulting multilooked images were stored as polarimetric covariance matrices C (in the Pauli basis) [18]. We adjusted the orientation angle [18], subtracted the thermal noise power obtained from the Radarsat-2 metadata, and calibrated the backscatter magnitude to σ 0 using a tangentplane approach [32] with elevation information derived from the TanDEM-X DEM (12 m resolution) [33]. These corrected covariance matrices served as the basis for all analyses of and tall hillslope shrubs (TS-2) reveal a large spatial variability. c) Limited temporal variability observed in the open tundra at the permanent OT-UP plot, as the five probes (gray) and their mean (black) vary less over time than the spread between the probes at any given moment. the polarimetric and backscatter characteristics. Finally, the results were geocoded to a geographic reference system using the TanDEM-X DEM.
An airborne L-band (24 cm wavelength) image complements the higher-frequency C-band data. The fully polarimetric image was acquired by NASA's UAVSAR system on 13 Sept 2017. The incidence angle varied from 22 • at near range to more than 60 • . Unfortunately, the vicinity of the MM station was only imaged partially, at the near-range edge of the image; to characterize the tundra's backscatter response at higher incidence angles, we also analyzed areas to the west of our study sites (W in Fig. 1). To provide a fair comparison between L-and C-band data, we artificially degraded the resolution of the UAVSAR data so that the resolution matched that of the Radarsat-2 data. Otherwise, the radar processing paralleled that of the C-band data.

IV. POLARIMETRIC RESPONSE OF THE TUNDRA
Polarimetric information can help to disentangle the scattering contributions from the vegetation (volume) from the underlying soil (surface), thus facilitating soil moisture estimation [16]. A prerequisite for the separation of the two components is an understanding of the scattering response across a gradient of shrub density. To account for the spatially variable microtopography, we introduced a tailor-made decomposition method.

1) Scattering response across shrub gradient
To characterize the scattering behavior of the landscape at C and L-band, we derived the Cloude-Pottier entropy H and alphaᾱ parameters from the eigen-decomposition of the polarimetric covariance matrices C [18]. The entropy reflects the randomness of the scattering process, with H = 0 corresponding to a deterministic target and H = 1 to one of maximum polarimetric randomness, which contains no further polarimetric information. The entropy is expected to increase with both surface roughness and vegetation density [16]. The mean alphaᾱ parameter characterizes the scattering mechanisms. Low values ofᾱ ≈ 0 are indicative of surface scattering and are thus expected for barely vegetated terrain.

Conversely, intermediate values ofᾱ ≈ π
4 are indicative of dipole volume scattering [18]. We focused on the three most important land cover types: open hummocky tundra, dwarf shrubs and tall hillslope shrubs ( Fig. 1).
2) X-Bragg model for rough surface scattering To interpret the observed polarimetric response of the bare open tundra, we turned to the X-Bragg model. It can predict the polarimetric response from a surface with microtopography [18], [34], [16]. Its great strength is its ability to reproduce the depolarization observed over rough surfaces [16]. It is a two-scale surface scattering model: the micro-scale roughness h of a small patch gives rise to zero-entropy Bragg scattering (small perturbation model, SPM), but the combination of many such patches on a hummock-interhummock sequence results in polarimetrically random scattering. In the Pauli basis [18] (1) f s controls the overall backscatter magnitude, and κ describes the scattering mechanism of an individual patch. κ depends, according to the SPM, only on the dielectric constant, which in turn varies with soil moisture as described by a mixing model like that of Mironov [35]. This makes κ a potentially useful quantity for soil moisture estimation: as the soil moisture increases, the real part Reκ becomes increasingly negative, corresponding to larger α angles [16]. The microtopography is represented by ψ, the maximum slope along a hummock, which is an indicator of the largescale surface roughness. With increasing ψ, the scattering becomes successively more random (H deviates increasingly from 0).
To assess the applicability of the X-Bragg scattering model over the open tundra, we compared the observed H andᾱ values to those that the model can predict across all feasible values of ψ (0 to π 2 ) and all values of κ, which are simulated using the Mironov mixing model for Arctic organic soils [18], [35]. A priori, the validity is restricted by the microscale roughness h < λ 2π . As this limit is approached or exceeded, higher-order scattering becomes increasingly important. Ballester-Berman et al. [36] showed how the higher-order scattering can largely be accounted for by a volume-like highentropy return that is commonly used to represent vegetation like shrubs.

3) MTV decomposition
To isolate the surface contribution from volume contributions, we introduced a novel decomposition. It is very similar to the approach by [36] in that it explicitly estimates the variable large-scale roughness ψ, here associated with the hummocky terrain sketched in Fig. 3. It accounts for the rough surface and volume scattering by combining the X-Bragg surface model (1) with a random volume of dipoles of magnitudef v [18], [16]: In this model the relative importance of surface scattering is given by However, attributing the volume component to vegetation scattering is difficult in the presence of higher-order surface scattering, as the higher-order contribution that is not included in the X-Bragg model is polarimetrically indistinguishable from vegetation scattering in the framework by [36] (see Sec. IV-A2). The surface contribution was, on the other hand, shown to be amenable to soil moisture estimation [36]. A priori, a random volume model seems appropriate as the shrubs do not exhibit conspicuous canopy orientation, and the observedᾱ ≈ 45 • over tall shrubs at C-band are also consistent with random dipoles. Note that there is no dihedral scattering mechanism because the lowᾱ 90 • do not indicate an important double bounce contribution at C-band. All these assumptions are subject to reservations. Testing the assumptions with observations is, however, difficult because the inversion is not overdetermined. An observed covariance C matrix provides five constraints (assuming reflection symmetry), and we want to determine the five parameters summarized in Fig. 3: f s , the surface scattering parameter κ (real and imaginary part), ψ, and the volume powerf v .
To solve for all MTV model parameters, it is expedient to rearrange the matrix expressions into a five-dimensional vector of real observables x: We split the estimation in two parts. The first part solved ] T by formulating it as a least-squares problem subject to the physical constraints ψ ≤ π 2 , |κ| < 1, f s &f v ≥ 0. We applied the trust-region quasi-Newton approach by Voglis and Lagaris [37], using initial values obtained by assuming ψ 0 = π 8 and solving for the other parameters. This yielded an estimate of u. In the second part of the estimation, the remaining parameter arg(κ) was computed from x 5 .
To assess the accuracy of the inversion, we conducted a simulation study. For each scenario, we simulated 5000 observable covariance matrices C 0 subject to speckle (40 and 160 looks) and estimated the parameters for each C 0 . Three scenarios served to study the accuracy as a function of the prescribed scattering characteristics. The first scenario Vegetation: low IA was meant to illustrate the impact of an increasing vegetation contribution at low incidence angles. To this end, we varied its relative contribution η from (3) while keeping the prescribed ψ 0 small ( π 10 ) and κ = −0.1, a typical value for low incidence angles. The simulation results in Fig. 4 show that the estimation of Reκ and f s is accurate for small η (Fig. 4). However, it quickly deteriorates as η and with it H increase, as relative errors of a factor 2-5 are common, irrespective of the number of looks. The errors are also large in the second scenario Vegetation: high IA, whose value of κ = −0.2 is more characteristic for higher incidence angles (Fig. 4). This scenario formed the baseline for all further scenarios. The third sensitivity scenario, Roughness, was intended to elucidate the dependence of the retrieval accuracy as a function of the microtopographic roughness expressed by the angle ψ 0 for moderate vegetation (η = 0.5). Fig. 4 illustrates the small impact of varying ψ 0 on the accuracy in this scenario.
In summary, the biggest problem is high-entropy scattering, which we will see is dominant in the tundra at C-band. The reason for the difficulties is that both a surface component with large |ψ| and the vegetation return induce high-entropy scattering that is (nearly) azimuthally symmetric. Consequently, for high entropies the attribution to these two terms is highly sensitive to the observed C, as our simulations suggest that relative errors in Re κ and f s of >100% occur just due to speckle (e.g. Fig. 4, Vegetation: low IA). Especially the low accuracy with which κ can be retrieved is a problem for soil moisture inversion. These problems are not restricted to the MTV estimation, as they have also been documented for related decompositions [25].
We also tested the robustness to the model parametrization by applying the MTV inversion to simulated data that were generated by extensions of the MTV forward model. When us-ing a raised-cosine instead of the X-Bragg uniform distribution for the microtopographic roughness in the forward simulations (scenario: Surface model in Fig. 4), the estimates based on a uniform distribution barely change. They are similarly barely affected (e.g. <10% for Reβ) in the scenario Heterogeneity: bare, in which we let dielectric heterogeneity induce a range of surface scattering mechanisms (within κ ± ∆κ; shown for a bare soil with η = 0). Next, however, a misspecified volume particle shape cannot be handled (scenario: Volume shape), as the inversion cannot disentangle the surface and volume contributions correctly. To this end, we replaced the random dipole volume model by a generic azimuthally symmetric volume model (normalized to span 1, so that the meaning of f v does not change, cf. [18]) The response is governed by the shape parameter R p , which varies between zero (dipole) and one (sphere). It corresponds to the ratio of the polarizability of the spheroidal particles in the Rayleigh limit. We find that, as the particles become more spherical (R p → 1), the estimation of both κ and f s breaks down. A similar problem occurs in the final scenario (Volume orientation), where we dropped the random volume assumption and replaced it by an oriented volume. To this end, we replaced the random dipole volume model by an oriented dipole model [18], where the orientation angle of the vegetation particles, ψ v , varies from 0 (perfectly aligned in the horizontal direction) to 0.5 π (random volume):

4) Sensitivity to soil moisture: time series analysis
To explore the dependence of polarimetric parameters on the incidence angle θ and on soil moisture v at C-band, we applied regression analysis. We expressed a radar observable y (Cloude-Pottier parameters: H andᾱ; MTV parameters: η, f s ,f v , Reκ; and VV σ 0 ): We thus accounted for the influence of θ (angular slope b θ ), v (surface soil moisture dependence b v ) and their interaction (b θ·v ). The parameters were estimated using ordinary least squares. As time series were required, we only applied regression analysis to the C-band data. Due to the additional requirement of soil moisture data, we applied the regression analysis at the OT-MM (open hummocky tundra) site, and we did so for the 2016 data because of the superior availability of Radarsat-2 data.
B. Results and discussion

1) Scattering characteristics
The open hummocky tundra is characterized by highentropy scattering at C-band, contrary to what would be expected for barely vegetated terrain. At the two open tundra plots, the H-ᾱ plots in Fig. 5 show H values of 0.6-0.8 ( Fig.  5: 2016-08-16 at 22 • , 2016-08-15 at 45 • ), increasing with incidence angle. These high values are due to the physical scattering behavior, rather than due to measurement noise (detailed analysis not shown). The mean alpha angleᾱ shows that the scattering is approximately azimuthally symmetric in the bare tundra [18]. This is indicated by values close to the lower boundary of feasibleᾱ for a given entropy (Fig. 5).
At L-band, the bare tundra looks more like a typical surface, with entropies of 0.1 to 0.2 (Fig. 5). Also the alpha angles are largely consistent with surface scattering. The comparison across incidence angles is complicated by the fact that the values shown in Fig. 5 derive from different areas. Specifically, the shallow incidence angle data were taken from areas with analogous land cover further west, as the main study area was only partially imaged (and at steep incidence angles). These indicators suggest that L-band data will be more amenable to soil moisture estimation using purely polarimetric methods than C-band data.
There is little contrast between the open tundra and shrub patches at C-band, underscoring the volume-like scattering of the open tundra. At C-band, the scattering at dwarf or tall shrubs is slightly more random than over the open tundra, with H > 0.7 (Fig. 5). At such high entropies, theᾱ angles are not very meaningful, as they are constrained to a small range. However, they appear to be on the low end of the feasible range, consistent with the MTV assumption of negligible dihedral scattering. The low contrast in polarimetric scattering characteristics between the open and the shrub tundra is also clearly evident in the RGB Pauli and entropy images of Fig.  6. They show a heterogeneous region 10 km west of the MM site (Fig. 1a), where both C-band and the L-band data are available at an incidence angle of ∼ 45 • . They underscore the high entropies at C-band across the landscape, apart from the lakes and, to a lesser extent, the polygonal fields. Conversely, at L-band the different shrub communities are more clearly distinguished from the open tundra sites, as the increase in entropy from shrub-free (H ≈ 0.1) to shrub sites (∼ 0.4) is much larger, Fig. 5. Theᾱ angles deviate markedly from the lower bound (azimuthal symmetry), as they move towards the dipole scattering region.
2) Applicability of the X-Bragg model The scattering from the open hummocky tundra cannot be reconciled with the X-Bragg surface model at C-band. The open tundra plots outside the range of potential soil responses in Fig. 5 (X-Bragg with Mironov mixing model; the lower/upper blue line represents a dry/saturated soil as ψ varies). The observed entropies are higher than the model predicts, whereas at L-band the model fits much better. What else could induce such highly random scattering? One possibility is that it is the variability in soil moisture and thus the polarimetric scattering mechanism κ across the microtopography, but our model results of Fig. 4 (Heterogeneity:bare) suggests that this mechanism cannot account for such high H. Alternatively, volume scattering from either the very limited above-ground biomass or from the subsurface could account for it. Finally, also higher-order surface scattering could contribute. We will   now turn to the MTV decomposition to tease out the surface scattering component.

3) MTV decomposition
Also the MTV decomposition highlights the volume-like scattering response over the tundra at C-band. While the estimated surface contribution decreases with increasing shrub cover in Fig. 6 and 7, there is limited contrast between shrubs and the open tundra. In particular, the surface contribution does not dominate (η ≈ 0.5) over the open hummocky tundra, despite the very limited shrub cover. The decomposition's restricted applicability at high H becomes evident by the very noisy appearance of η in Fig. 6. At C-band, dominant surface scatter is only observed over the road, the lakes, and certain polygonal fields. Conversely, at L-band, the surface contribution always dominates in the open tundra in Fig. 6. The L-band relative surface contribution η drops to values close to 0 over dense and tall riparian shrubs, with intermediate values of 0.5-0.7 over the hillslope tall shrubs and the dwarf shrubs. Overall, at C-band (in contrast to L-band) the high-entropy, volume-like scattering indicates a lack of information content for purely polarimetric approaches.

4) Sensitivity to soil moisture: time series analysis
Time series analysis at C-band suggests that the backscatter magnitude is more amenable to soil moisture estimation than purely polarimetric indicators. According to our regression analysis using (7), the surface backscatter f s is a good predictor of surface soil moisture (Tab. I). Also the VV backscatter and the volume powerf v increase with surface soil moisture, but the relationships are weaker. Conversely, the relative surface contribution η exhibits a very weak relation with soil moisture (Tab. I). The other polarimetric indicators, Reκ, as well as H andᾱ do not show a significant soil moisture signal either. The magnitude thus appears to be promising for soil moisture estimation, and understanding its angular dependence is crucial for this task.

V. ANGULAR CHARACTERISTICS
The angular dependence of the backscatter is another important indicator of the scattering behavior, and it is particularly useful for distinguishing surface from volume scattering. For surfaces, the backscatter decreases more quickly with θ than for volumes, especially when the surface is smooth. Further, understanding the angular dependence is of great value for soil moisture estimation, as satellite data are often acquired across a range of incidence angles.

A. Methods
To estimate the angular dependence at C-band, we used ordinary least squares (OLS) regression. As opposed to the regression model of (7), we only included the angular dependence i.e. y = a + b θ (θ − 30 • ), which allowed us to apply the approach to the entire study region. We thus implicitly assumed that there was no systematic relation between θ and the soil moisture. At L-band, we had to resort to comparing similar land cover types and near and far range to estimate the angular dependence.
Two semi-empirical surface scattering models were assessed in terms of their predicted angular characteristics. First, we studied the one-parameter Baghdadi model, a semi-empirical model that applies to a wide range of roughness conditions  [22]. It has been found to accurately represent the angular characteristics of bare mineral soils for a wide range of roughness conditions, especially in comparison to theoretical models of varying complexity. It is parameterized in terms of the dielectric constant, which in turn depends on soil moisture; here, we used the Mironov model for Arctic organic soils [35]. The Baghdadi model's only remaining parameter is the roughness height h, which we fitted to the observations over the open tundra by choosing h so that the angular dependence matched, to first order, the observed one [22]. Second, we replicated the analyses for the Oh model [21], a one-parameter (h) semi-empirical model that was derived from measurements of very rough soils. One crucial aspect of the angular dependence is the extent to which the dynamic range (due to soil moisture variations) varies with incidence angle. Many models fail to capture the very limited decrease of the dynamic range for steeper incidence angles, thus limiting their suitability for soil moisture estimation from multi-angular time series data [17]. Conversely, when the dynamic range stays essentially constant, a simple linear model is often more suitable.

B. Results and discussion
The backscatter decreases slowly with incidence angle (0.11 dB/deg at VV, Fig. 7c) over the open tundra. This decrease is much more subdued than what most previous studies have found over bare or sparsely vegetated (shrubs, crops, grass) areas (Fig. 7e). To be more precise, this holds for the hummocky terrain in most of the study area, rather than over the icewedge polygons (Fig. 7c). The angular slope is even smaller in magnitude over shrub patches, with values comparable to previous observations over forest [17], [40]. Similar to the polarimetric characteristics, at L-band the angular slope conforms much more to expectations. The results in Fig. 7c show that the angular decrease is larger for all land cover classes. It is largest for the open tundra, which acts much more like a surface scatterer than at C-band.
The angular characteristics of the open tundra are characterized by a slow linear decrease at a constant low dynamic range, which the semi-empirical models can only partially reproduce. According to the models, the observed backscatter magnitude and angular slope require that the soils be rough (large h). The Baghdadi model can reproduce them for a roughness height of h = 6 mm (Fig. 8). However, it wrongly predicts that the variability (due to soil moisture changes) decreases rapidly with θ. In the data, the dynamic range remains approximately constant. This is consistent with a near-constant soil moisture sensitivity because the soil moisture distribution at the acquisition times was essentially the same for θ above or below 30 • according to the in-situ data (difference in mean <1%, difference in variability ≈ 10%). It is also indicated by the small (and non-significant) interaction term b θ·v in the VV regression analysis (Tab. I). Conversely, the observations in Fig. 8 indicate that the Oh model can replicate these properties more accurately. At such high h ≈ 1 cm, increasing h leads to an increase in σ 0 that is almost independent of θ. It also more accurately captures the observed quasi-linear trend with incidence angle, whereas the Baghdadi model predicts a noticeably nonlinear dependence for wet soils.
The linear angular dependence and the constant variability over the open tundra lend themselves to a simple linear model that is amenable to soil moisture estimation.

VI. SOIL MOISTURE INVERSION
To estimate soil moisture at C-band, we adapted the widely employed linear empirical model, which is based solely on the backscatter magnitude [17], [23]. It has two key advantages: i) it does not require external soil moisture or vegetation data for calibration, and ii) it is very flexible in that it can be applied -at least in theory -to any surface irrespective of its roughness or vegetation cover, provided these characteristics remain constant. A disadvantage is that only a relative index of soil moisture can be retrieved. We adapted this model to account for the high noise (at 160 looks: 0.3 dB) relative to the limited dynamic range observed over our tundra site (<2 dB, Fig. 8). The limited dynamic range is likely associated with the low temporal variability of soil moisture (Fig. 1). To increase the signal-to-noise ratio, we embedded it in a Bayesian hierarchical framework, whereby neighboring grid cells are dependent. The dependence reflects the temporal similarity of the backscatter. While the large microtopographical variability of soil moisture cannot be resolved, the idea is to exploit the temporal correlation of soil moisture and backscatter on larger scales of several hundred meters in response to the atmospheric forcing [14]. The degree of correlation, and hence the amount of pooling of information across pixels, is estimated in an adaptive way, thus reducing the noise level while keeping all the advantages of the linear model.
The reasons for focusing on the empirical backscatter model are that the polarimetric and angular characteristics identified above render more physical approaches inappropriate. In particular, they illustrate why purely polarimetric approaches failed in the hummocky tundra. Instead of discussing the results obtained with two polarimetric approaches in detail, we only provide a quick summary of the problems. First, we tried the multipolarimetric Oh model inversion for bare surfaces [21], as the Oh model could replicate the most important angular characteristics (Sec. V-B). However, the inversions were not successful because of the model's inability to predict the HV backscatter (too small at large incidence angles). Second, the κ value of the MTV surface return, predicted to vary with soil moisture, was too noisy and commonly outside the physical range of the SPM/X-Bragg model (Fig.  9). The noisy nature was expected as the simulations in Sec. IV suggested that κ estimation is infeasible for high-entropy scattering (Fig. 4). The problem of observed Reκ values exceeding 0 could be partially mitigated by using a Fresnel instead of the SPM/X-Bragg model to estimate soil moisture from κ [36], but the inability to reliably estimate κ in the open tundra precluded such an approach. This is not to say that polarimetric information is useless: in our soil moisture retrievals using the linear model, we thus compare the MTV surface component with the VV backscatter.

A. Methods
The linear model relates the observed backscatter σ 0 ij at resolution cell i for acquisition j to the soil moisture v ij in a linear fashion [17], [23]. In our probabilistic framework, it reads It further assumes that the dependence on the incidence angle θ ij is linear with slope β i , and that the variability is dominated by speckle noise, whose variance s 2 L depends on the number of looks L. We focused on the VV channel, because it is widely available; in particular, the Sentinel-1 satellites started acquiring VV Interferometric Wide Swath Data across the Canadian tundra in late 2016. Additionally, we studied the MTV surface scattering component, not only because the residual influence of volume scattering is potentially reduced, but also because it had a higher R 2 with the soil moisture observations (Tab. I).
The forward model (8) rests on several crucial assumptions. First, it assumes that the backscatter at a given location increases linearly with soil moisture. Indeed, our preliminary regression analysis showed a positive, approximately linear relation (Tab. I). Second, the coefficients are assumed constant in time, which seems reasonable owing to our focus on July and August, the period after leaf out but before freeze-back. Third, the θ-dependence is assumed to be linear, as there is no indication for strong higher-order terms (Fig. 8, Sec. V). Finally, it assumes there are no interactions between θ and v. This implies that the dynamic range of the observed backscatter should be independent of θ. Fig. 8 shows that the observed dynamic range does not vary substantially with θ. An approximately θ-independent soil moisture sensitivity is also commonly found in temperate regions [17].
Inferring soil moisture using a time-series approach is plagued by measurement noise and by problems of nonuniqueness [30]. We attempted to address these problems by estimating all parameters and v ij in a one-step Bayesian procedure [44]. It estimates the posterior probability of the parameters and the soil moisture given the observed σ 0 ij using Hamiltonian Monte Carlo sampling. The Bayesian estimation reduces the impact of measurement noise, which is comparatively large due to the high resolution (limited multilooking) in comparison to the small dynamic range observed in the tundra (Fig. 8), by pooling information from neighboring resolution cells. The mathematical details are provided in the appendix. The non-uniqueness cannot be avoided, but we implemented two strategies to cope with it. First, the Bayesian approach imposed regularization via the prior distribution. It is, however, contingent on assumptions about the true soil moisture distribution; in particular, that v ij and θ ij are independent, enabling the identification of β i (see Appendix).  This assumption seemed reasonable as the association between the in-situ measurements and θ is small (Sec. V). Second, we used validation metrics that are insensitive to an overall shift or scaling of the soil moisture time series, neither of which can be identified uniquely. For instance, an overall increase in v ij has the same effect as an increase in µ i . We applied the Bayesian inversion procedure over two small (about 3 km 2 ) study regions shown in Fig. 1 (gray rectangles: MM and UP), at the center of each was one of the instrumented plots. The study regions were chosen so as not to include lakes, because the Bayesian model would likely have to be adapted to handle the lakes well. We validated the retrievals at these two instrumented plots in the open tundra by comparing the retrievals to the spatially aggregated in-situ observations. The validation metrics we employed were the correlation coefficient ρ and the calibrated root mean square error cRMSE, which is the RMSE between the in-situ surface soil moisture and the SAR estimates calibrated using ordinary least squares (offset, dynamic range).

B. Results and discussion
The soil moisture estimated using the VV time series approach tracks the in-situ measurements to some extent in Fig. 10a, but there are issues. First of all, the dynamic range of the in-situ measurements series is very small (<0.1 m 3 m −3 ), making it difficult to pick out the temporal agreement between the two data sets. The temporal sampling of the Radarsat-2 acquisitions compounds the issue, as some of the largest wetting/drying events are not adequately sampled. Second, there are a few conspicuous outliers, and these tend to be associated with precipitation events. The more common type is an underestimation by Radarsat-2 compared to in-situ data (three such occasions are marked by diamonds in Fig. 10a). The observation on 2016-07-21 illustrates the phenomenon (OT-MM), as the Radarsat-2 soil moisture indicates dry conditions despite recent rainfall. The low soil moisture estimate is due to low backscatter: it is 1.8 dB lower than 12 days previously, even though θ is the same and the in-situ soil moisture values are comparable (|∆v| < 0.01). While the change in VV backscatter is large compared to 2-3 dB dynamic range, the polarimetric indicators are almost constant (|∆ᾱ| < 1 • , ∆H < 0.01; not shown), illustrating the difficulty in detecting such an event based on time series data.
The validation metrics are 0.3-0.5 for the correlation ρ and 0.02 m 3 m −3 for the calibrated RMSE on average (Fig. 10b). When interpreting these metrics, it is important to keep the small dynamic range during the study period in mind. For both metrics, there is little difference between the results for VV and those for the MTV surface component f s , again on average. For both polarizations, there is some variability across years (2014, 2016) and sites (OT-MM and OT-UP). Finally, the incidence angle seems to have a minor impact on the quality of the retrievals, as there are no conspicuous differences in variability when plotting the retrievals against the in-situ data in Fig. 10c.
To explore the soil moisture estimates over shrub-covered surfaces, we had to compare them to the in-situ measurements over the open tundra rather than within the shrub patches, due to lack of in-situ measurements in the shrub patches. The estimates over the dwarf shrubs (DS-1) and the tall shrubs (TS-2) have lower correlations (0.26 and 0.31, respectively) with the in-situ measurements in the open tundra (OT-MM) than the estimates over the open tundra. Overall, the time series of retrieved soil moisture for the different shrub communities plotted in Fig. 11 are similar. For instance, the correctly estimated drying trend during July is common to all. However, the lack of in-situ measurements precludes a more in-depth assessment.

VII. DISCUSSION AND CONCLUSIONS
We analyzed the potential of C-band SAR remote sensing to estimate soil moisture in the subarctic tundra. We identified two challenges that soil moisture inversion in the subarctic tundra faces. One is to do with the complex scattering in terms of the polarimetric and angular response, the other one with the spatiotemporal variability of soil moisture. Regarding the latter, we observed large spatial variations in soil moisture on the sub-resolution scale (<100 m). Consequently, we required 3-6 in-situ measurements for validating SAR retrievals within a single resolution cell (∼ 100 m). Conversely, the observed temporal variability was small (<0.1 m 3 m −3 ).
The microwave scattering signal in the subarctic tundra is complex and difficult to interpret. The difficulties are not only to do with the shrubs, whose structure, density and biomass varies across the different shrub communities. After all, the polarimetric characteristics differ only marginally between dense shrubs and the open tundra at C-band. Despite the very low biomass in the open tundra, we observed high entropies at C-band (H ≈ 0.7, Fig. 5), which are considered unusual for barely vegetated surfaces [18]. They are likely related to the pronounced hummock-interhummock microtopography. However, the entropies are much larger than what would be expected theoretically (X-Bragg model) for a very rough surface, so that the organic soils (subsurface scattering) and the limited above-ground biomass also likely play a role. The very rough and volume-like scattering characteristics of the open tundra also became evident in terms of the incidence angle behavior, which was more akin to previous observations over forest than over bare soils (Fig. 7).
The high entropies made soil moisture estimation using purely polarimetric methods impracticable for our study site. We found that we could not reliably isolate the ground contribution using our tailor-made MicroTopography and Vegetation (MTV) decomposition approach that could account for the microtopography. This lack of robustness was due to the low information content (high H), according to our theoretical analyses and simulations. The low information content also limits the applicability and range of values of complementary polarimetric quantities such as the alpha angle or the anisotropy.
This left the backscatter magnitude as the main signal for soil moisture estimation. Indeed, the magnitude tracks the overall wetting and drying over the open tundra. To account for the limited dynamic range of both backscatter and soil moisture, we embedded a simple time-series method in a Bayesian framework within which the noisy observations from adjacent pixels were partially pooled to increase the precision. The retrievals achieved correlations of 0.3-0.5 (Fig. 10b) with in-situ surface soil moisture, for both the polarimetrically derived MTV surface power and the VV backscatter. These correlations are not as high as one would want, but they are not uncommon in practice [17], [16], [45], [20]. As the time-series approach can be applied to single-pol VV data, it appears to be promising for the future given the recent availability of frequent Sentinel-1 data at constant incidence angles. One open question in this context is the decrease in backscatter that we observed in response to precipitation events. This is a challenge for soil moisture estimation in general, as the decrease is contrary to the usual increase in backscatter with soil moisture. While the data set is small and this behavior was not universally observed, we conjecture that this may be caused by sub-surface scattering. As the surficial organic material becomes wet, it attenuates the return from the underlying denser organic materials and mineral soil, even though it also becomes a stronger scatterer [19]. To better understand the hummocky tundra's backscatter response to soil moisture changes, we recommend dedicated observations of moisture and backscattering profiles, and tailor-made scattering models.
For the future, L-band SAR seems more promising. The polarimetric and angular characteristics we observed are more in line with what would be expected from sparsely vegetated terrain (Fig. 5, 6). However, the limited data at our disposal do not rule out several of the problems identified at C-band. In particular, the subsurface contributions may also be critical because of the increased penetration into organic soils at larger wavelengths. Overcoming these challenges is critical, because fine-scale soil moisture estimates can make vital contributions to ecology, biogeochemistry and permafrost research.

APPENDIX
The Bayesian time series approach exploits the spatial similarity across the study region to increase the precision of the soil moisture estimates. It does this by pooling information across resolution cells and acquisitions. The pooling was induced by postulating that the model parameters (e.g. µ i ) and the v ij were similar in space, where a cell is denoted by i and a Radarsat-2 acquisition by j. The degree of similarity is part of the model structure and was estimated within the inference. The model parameters µ i , β i and γ i were each assumed to be distributed according to a normal distribution N (ϕ mean , ϕ 2 std ), with each of the ϕ population parameters being themselves assigned a weakly informative prior distribution (Tab. II). In the inference, these parameters were also updated: for instance, by adjusting the standard deviation φ µ,std , the degree of similarity of the µ i could be learned from the data.
For the soil moisture, we pooled information across space and time by decomposing v ij into a spatial average w j and a cell-specific anomaly u ij : where Φ = 0.8 is the porosity and π i describes the propensity of cell i to reflect the regional-scale soil moisture. The factor π i as well as u ij , w j can take on values between 0 and 1 and were assumed to be distributed according to Beta(φ π,alpha , φ π,beta ). The same restriction and an analogous population distribution applied to the degree of saturation values w j and u ij , see Tab. II. We estimated the posterior distribution by Hamiltonian Markov Chain Monte Carlo (MCMC) as implemented in the pymc3 package [46]. Following common practice, we simulated four chains of 2000 samples each and discarded the first 1000 samples.