Close
About
FAQ
Home
Collections
Login
USC Login
Register
0
Selected
Invert selection
Deselect all
Deselect all
Click here to refresh results
Click here to refresh results
USC
/
Digital Library
/
University of Southern California Dissertations and Theses
/
Radar remote sensing of permafrost landscapes and active layer soil properties
(USC Thesis Other)
Radar remote sensing of permafrost landscapes and active layer soil properties
PDF
Download
Share
Open document
Flip pages
Contact Us
Contact Us
Copy asset link
Request this asset
Transcript (if available)
Content
Radar Remote Sensing of Permafrost Landscapes and Active Layer Soil Properties by Richard H. Chen A Dissertation Presented to the FACULTY OF THE GRADUATE SCHOOL UNIVERSITY OF SOUTHERN CALIFORNIA In Partial Fulllment of the Requirements for the Degree Doctor of Philosophy (Electrical Engineering) December 2019 Doctoral Committee: Professor Mahta Moghaddam, Chair Professor B. Keith Jenkins Professor Andreas F. Molisch Assistant Professor Kelly T. Sanders Research Assistant Professor Alireza Tabatabaeenejad To my dad and my sister. To Ying and Renee. ii TABLE OF CONTENTS DEDICATION : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : ii LIST OF FIGURES : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : v ABSTRACT : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : xii Chapter I { Introduction 1.1 Background and Motivation . . . . . . . . . . . . . . . . . . . . . 1 1.2 Retrieval Process of Subsurface Geophysical Properties . . . . . . 6 1.3 Objectives . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 1.4 Contributions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8 Chapter II { Retrieval of Permafrost Active Layer Properties: A Time- Series Approach 2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 2.2 In-Situ Data and Time-Series Considerations . . . . . . . . . . . . 11 2.3 Radar Forward Scattering Model . . . . . . . . . . . . . . . . . . . 13 2.4 Radar Sensing Depth . . . . . . . . . . . . . . . . . . . . . . . . . 18 2.5 Sensitivity of Radar Backscatter to Subsurface Variables . . . . . . 21 2.6 Inverse Problem . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22 2.6.1 Cost Function . . . . . . . . . . . . . . . . . . . . . . . . . 23 2.6.2 Largest Possible Depths . . . . . . . . . . . . . . . . . . . . 25 2.6.3 Time-Series Retrieval . . . . . . . . . . . . . . . . . . . . . 30 2.7 Retrieval of Active Layer Properties . . . . . . . . . . . . . . . . . 34 2.7.1 Retrieval Results from Pre-ABoVE IDS Flights . . . . . . . 36 2.7.2 Retrieval Results from 2017 ABoVE Flights . . . . . . . . . 44 2.8 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49 Chapter III { Harmonized Parametrization of Soil Properties for Multi- Physics Approach 3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51 3.2 Soil Parametrization Scheme . . . . . . . . . . . . . . . . . . . . . 53 3.2.1 Incorporating Fiber Content into Parametrization . . . . . 53 3.2.2 Parametrization of Soil Hydraulic Properties . . . . . . . . 56 iii 3.3 Soil Dielectric Mixing Model . . . . . . . . . . . . . . . . . . . . . 61 3.4 Model Validation . . . . . . . . . . . . . . . . . . . . . . . . . . . 65 3.5 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67 Chapter IV { Retrieval of Soil Moisture and Organic Matter Proles in the Arctic-Boreal Region 4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68 4.2 Prole Representation of Permafrost Soils . . . . . . . . . . . . . . 69 4.2.1 Organic Matter Proles . . . . . . . . . . . . . . . . . . . . 70 4.2.2 Soil Moisture Proles . . . . . . . . . . . . . . . . . . . . . 73 4.3 Radar Retrieval Process . . . . . . . . . . . . . . . . . . . . . . . . 74 4.3.1 Forward Model . . . . . . . . . . . . . . . . . . . . . . . . . 74 4.3.2 Inverse Problem . . . . . . . . . . . . . . . . . . . . . . . . 76 4.4 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77 4.4.1 Synthetic Data . . . . . . . . . . . . . . . . . . . . . . . . . 77 4.4.2 AirMOSS Data in 2017 ABoVE Airborne Campaign . . . . 82 4.5 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84 Chapter V { Joint Retrieval of Active Layer Proles Using PolSAR Backscat- ter and InSAR Deformation 5.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85 5.2 Forward Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 86 5.2.1 ReSALT Subsidence Model . . . . . . . . . . . . . . . . . . 87 5.2.2 PolSAR Backscatter Model . . . . . . . . . . . . . . . . . . 91 5.3 Combined PolSAR-InSAR Retrieval . . . . . . . . . . . . . . . . . 93 5.4 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 95 5.4.1 Synthetic Data . . . . . . . . . . . . . . . . . . . . . . . . . 95 5.4.2 UAVSAR and AirMOSS Data in 2017 ABoVE Airborne Campaign . . . . . . . . . . . . . . . . . . . . . . . . . . . 97 5.5 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 100 Chapter VI { Conclusion and Future Work 6.1 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 101 6.2 Future Work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 103 6.2.1 Implementing a joint radar scattering and soil process inver- sion framework . . . . . . . . . . . . . . . . . . . . . . . . . 103 6.2.2 Estimation of Vegetation Structure Using SAR and Lidar . 104 BIBLIOGRAPHY : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : 107 iv LIST OF FIGURES Figure 1.1 Denition of permafrost and the active layer and the soil features related to seasonal freezing and thawing. Redrawn from Woo 2012 [1]. 2 1.2 AirMOSS radar swaths and permafrost observation sites over North- ern Alaska. CALM: Circumpolar Active Layer Monitoring program [www2.gwu.edu/calm]. UAF GIPL: Geophysical Institute Permafrost Laboratory, University of Alaska Fairbanks [permafrost.gi.alaska.edu]. 4 1.3 Flight transects own during the 2017 ABoVE airborne campaigns. 5 1.4 Retrieval process of subsurface geophysical properties. . . . . . . . . 6 2.1 In-situ data from the GTN-P site operated by UAF GIPL in Dead- horse, Alaska [2]. (a) Soil moisture and temperature proles taken by Stevens HydraProbe sensors at depths of 15, 30, and 60 cm. The data are extracted from the real-time data portal at http://permafrost.gi. alaska.edu/site/dh1. (b) Isotherms interpolated from the thermistor measurements at depths of 0.01, 0.03, 0.11, 0.18, 0.25, 0.34, 0.41, 0.57, 0.72, and 0.93 m. The vertical dashed lines mark the AirMOSS ights on August 29 and October 1 in 2015. . . . . . . . . . . . . . 12 2.2 Geometry of a multilayer dielectric structure (left) and the three-layer model (right) used in the retrievals of active layer properties. The multilayer structure shows a more general case of N L layers that can be used to represent more sophisticated stratigraphic features in the active layer. The three-layer model includes the surface active layer ( 1 ), the bottom saturated active layer ( 2 ), and the permafrost ( f ). The top layer dielectric constant 1 and thicknessz 1 can vary between mid-late August and early October to account for the changes in soil moisture or freeze-thaw state in the near-surface layer. . . . . . . . 14 v 2.3 Dielectric properties of major soil classes on the Alaska North Slope. (a) Real ( 0 ) and imaginary ( 00 ) part of soil dielectric constant and (b) the ratio of 00 to 0 over the entire volumetric soil moisture (m v ) range based on the soil porosity, which is equivalent to the soil moisture at saturation. The soil classes are extracted from the soil texture data by De Lannoyet al. [3] and there are ve major soil classes on the Alaska North Slope (#68, #75, #224, #245 for the mineral soil classes with low (SOM < 0.49 %) to moderate (1.1 % SOM < 15 %) SOM content, and #253 for the organic soil class with high SOM content (SOM 15 %)). The soil dielectric models by Peplinski et al. [4] and by Bircher et al. [5] are used to compute the dielectric constant for the mineral and organic soil classes, respectively. . . . . . . . . . . . 16 2.4 Simulated radar backscattering coecient of a two-layer dielectric structure with varying rst-layer interface (z 1 ). The real parts of layer dielectric constants are set to be ( 0 1 , 0 2 ) = (30, 5) in this example. 17 2.5 Example of calculating radar sensing depth s for a top dielectric layer of 0 1 = 30 with a target dielectric substrate 0 2 = 5 at z 1 . The inset shows a zoom-in view of how s is determined. n represents the radar calibration accuracy and is set to 0.5 dB in this case. . . . 19 2.6 P- and L-band radar sensing depth s and penetration depth p with varying top-layer dielectric constant 0 1 . The calibration accuracy n of the radar systems is assumed to be 0.5 dB. The sharp peak at 0 1 = 5 is because of the perfect match between the two layers and therefore lack of a scattering interface. Note that the sensing depth is considerably larger than the penetration depth for both L-band and P-band. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20 2.7 P-band radar backscattering coecients 0 as a function of top-layer dielectric constant 0 1 and depth of the top layerz 1 . Dashed-dot line: corresponding radar sensing depth s . . . . . . . . . . . . . . . . . . 21 2.8 P-band radar backscattering coecients 0 as a function of bottom- layer dielectric constant 0 2 and depth of the top layerz 1 . Dashed-dot line: corresponding radar sensing depth s . . . . . . . . . . . . . . . 22 2.9 Cost function values of one-layer (half-space) problems with dierent rst-layer dielectric constants 0 1 . The dashed line is the error criterion nominally set to 2 n , the variance of calibrated radar backscattering cross section. Vertical axis is truncated at L(X) = 1 for display simplicity. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25 2.10 Cost function values for a two-layer problem with ( 0 1 , 0 2 ) = (30, 5) and one unknown z 1 . (a) z true 1 = 0.4 m. (b) z true 1 = 0.32 m. In (b), red circles and black crosses are the local minima that occur at z 1 z true 1 and atz 1 >z true 1 , respectively. Note that there are two red circles located closely around z true 1 = 0.32 m. The dashed line is the error criterion nominally set to 2 n , the variance of calibrated radar backscattering cross section. . . . . . . . . . . . . . . . . . . . . . . 26 vi 2.11 (a) Simulated radar backscattering coecient of a two-layer dielectric structure with varying rst-layer interface z 1 and ( 0 1 , 0 2 ) = (30, 5). The regions where z 1 can be solved by the LPD approach are high- lighted. (b) z min [n], the distances between adjacent local minima of L(z 1 ) in Fig. 2.10(b). Red circles and black crosses correspond to the local minima in Fig. 2.10(b) that occur at z 1 z true 1 and at z 1 >z true 1 , respectively. Red star corresponds to n q in (2.5). . . . . . 27 2.12 Comparison of error performance of z est 1 . (a) LPD approach at n = 0.5 dB. (b) Global minimum approach at n = 0.5 dB. (c) LPD approach at n = 0.3 dB. (d) LPD approach at n = 0.1 dB. The gray shaded regions represent one standard deviation around the average values of z est 1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28 2.13 Example demonstrating cost function values of a two-layer problem with one unknown z 1 , using two snapshots (L(X(t 1 )) and L(X(t 2 ))) and the time-series combined (L(X(t 1 ;t 2 ))). ( 0 1 (t 1 ); 0 2 ) = (30; 5) at t 1 and ( 0 1 (t 2 ); 0 2 ) = (15; 5) at t 2 . The true value of unknown z 1 is set to z true 1 = 0:4 m. The vertical pink dotted lines correspond to the solutions estimated by each time snapshot LPD approach. The vertical black dashed-dot lines are the search space for the LPD- assisted time-series approach. . . . . . . . . . . . . . . . . . . . . . 31 2.14 Comparison of error performance ofz est 1 . (a) LPD-assisted time-series approach using t 1 and t 2 . (b) Global minimum approach using t 1 andt 2 . (c) LPD-assisted time-series approach usingt 1 ,t 2 andt 3 . (d) Global minimum approach using t 1 , t 2 , and t 3 . ( 0 1 (t 1 ); 0 2 ) = (30; 5); ( 0 1 (t 2 ); 0 2 ) = (15; 5); ( 0 1 (t 3 ); 0 2 ) = (25; 5). The gray shaded regions represent one standard deviation around the average values of z est 1 . n is set to be 0.5 dB. The vertical dashed lines are the expected sensing depths s for the snapshots t 1 , t 2 , and t 3 . . . . . . . . . . . . 33 2.15 Root-mean-squared error (RMSE) of the LPD-assisted time-series approach. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34 2.16 Retrieved maps of layer dielectric constants i , layer vertical positions z i , ALT, and surface roughness h using the AirMOSS data acquired for the Deadhorse ight line over the northern section of Dalton High- way in 2015. Similar maps have been generated for 2014 observations but not shown for the purpose of brevity. . . . . . . . . . . . . . . . 38 2.17 Uncertainty maps of retrieved layer dielectric constants i , layer ver- tical positions z i , ALT, and surface roughness h using the AirMOSS data acquired for the Deadhorse ight line over the northern section of Dalton Highway in 2015. . . . . . . . . . . . . . . . . . . . . . . . 39 vii 2.18 Distribution of the 2014 and 2015 AirMOSS retrievals (ALT, 0 2 , and h) for the Deadhorse ight line across dierent NLCD land cover types: 31-Barren land, 95-Emergent Herbaceous Wetlands, 72- Sedge/Herbaceous, 51-Dwarf Scrub, 52-Shrub/Scrub. The red cen- tral mark in each box indicates the median, and the bottom and top edges of the box indicate the 25th and 75th percentiles, respectively. The whiskers are extended to the 5th and 95th percentiles. . . . . . 40 2.19 Retrieved maps of bottom layer dielectric constants 0 2 and ALT for the ight lines own in 2015 over Atqasuk, Barrow, and Ivotuk. The East-West gradient in (a) is an artifact related to the wide range of incidence angles within the AirMOSS imagery, which is expected to improve in the future when enhanced calibration of antenna patterns becomes available. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41 2.20 AirMOSS ALT retrievals and CALM in-situ ALT comparison. The red solid (blue open) symbols represent the sites where the in-situ ALT is larger (smaller) than the radar sensing depth s . The error metrics are shown in the parenthesis after the site name in the legend for each site. The rst number in the parenthesis is the retrieval bias; the second number is the 2 statistic introduced in [6]. . . . . . . . 42 2.21 Retrieved ALT for the legacy ight lines and the ones over Toolik Lake and Teller, Alaska own during the ABoVE Airborne Campaign in 2017. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45 2.22 Interannual comparisons of August surface soil moisture and ALT maps retrieved from AirMOSS data acquired over the Deadhorse legacy line in 2014, 2015, and 2017. . . . . . . . . . . . . . . . . . . 46 2.23 AirMOSS ALT retrievals and CALM in-situ ALT comparison. (a) Validation plot of 2017 data. The purple (blue) symbols represent the sites where the in-situ ALT is larger (smaller) than the radar sensing depth s . The error metrics are shown in the parenthesis after the site name in the legend for each site. The rst number in the parenthesis is the retrieval bias; the second number is the 2 statistic introduced in [6]. (b) Validation plot of all acquisition years (2014, 2015, 2017). . . . . . . . . . . . . . . . . . . . . . . . . . . . 48 3.1 (a) FC vs. bulk density. (b) Bulk density vs. OM. (c) FC vs. OM. The sand fraction is 35% for this example. The curves are NOT derived from the USDA-NCSS data, which are shown here merely for comparisons. Note that the ber content in USDA-NCSS data is dened as a percentage by volume, not by mass. . . . . . . . . . . . 55 3.2 (a) Sub-phase volumetric fractions (f v ) of mineral, humus, and - brous materials and (b) Soil porosity () as a function of organic matter content (OM). The sand fraction is 35% for this example. The porosity at 100% and 45% of OM correspond to the values of bric and sapric peat reported in Letts et al. (2000), respectively. . 56 viii 3.3 Soil hydraulic properties as a function of organic matter content (OM). (a) Air-entry potential ( s ). (b) Exponent b of the Camp- bell water retention curve. (c) Saturated hydraulic conductivity (k s ). The mineral soil texture is to (sand;clay) = (35%; 10%) for this ex- ample. The values of bric and sapric peat are from the ones reported in Letts et al. (2000) with minor adjustments for b and k s to avoid too slow drainage in well-decomposed organic soils. . . . . . . . . . 60 3.4 Soil dielectric constant as a function of soil moisture and organic matter content (OM) at P-band (430 MHz). . . . . . . . . . . . . . 63 3.5 Volumetric fractions of soil water phases in the total volume of soil (v =V x =V soil ). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64 3.6 Soil dielectric constant as a function of soil temperature and organic matter content (OM) at P-band (430 MHz). The saturation fraction (S w == s ) is set to 70% in this example. . . . . . . . . . . . . . . 64 3.7 Comparison of in-situ and model soil dielectric constants. The num- bers in the parentheses in the legend are the Root Mean Square Error (RMSE) of each soil pit prole. . . . . . . . . . . . . . . . . . . . . 66 4.1 Comparison of in-situ organic matter proles and two prole repre- sentations: exponential and sigmoid functions. The numbers in the parentheses in the legend are the RMSE of each tted prole. . . . 71 4.2 Comparison of in-situ bulk density ( b ) proles and two prole repre- sentations: exponential and sigmoid functions. The numbers in the parentheses in the legend are the RMSE of each modeled prole. . . 72 4.3 Comparison of in-situ saturation fraction (S w ) proles and the prole representation in (4.4). The numbers in the parentheses in the legend are the RMSE of each tted prole. . . . . . . . . . . . . . . . . . . 73 4.4 Example of prole representation of soil compositional properties (or- ganic matter content (OM), porosity (), bulk density ( b )), satura- tion fraction (S w ) and soil moisture (), and the corresponding soil dielectric prole (). The prole parameters are set to (A, M, S w0 , WTD) = (0.7, 0.2 m, 0.2, 0.2 m) for this example. . . . . . . . . . . 75 4.5 Geometry of a multilayer dielectric structure (left) representing the active layer and underlying permafrost, where (z) are discretized with a spacing of 1 cm and z ALT is the ALT. . . . . . . . . . . . . . 75 4.6 Prole models used in the retrieval of active layer soil properties. The organic matter prole is assumed to be time-invariant. For satura- tion proles, the surface saturation fraction (S w0 ) and water table depth (WTD) can vary between mid-late August and early October, whereas the bottom of the active layer remains fully saturated. . . . 77 4.7 Comparison of true and estimated proles of organic matter, poros- ity, saturation fraction, and soil moisture for the case ofM = 0:05m. Other prole parameters are set to (A, S w0 (t 1 ), WTD(t 1 ), S w0 (t 2 ), WTD(t 2 )) = (0.8, 0.3, 0.2 m, 0.2, 0.3 m). The shaded regions rep- resent one standard deviation around the mean estimated values for each parameter. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78 ix 4.8 Comparison of true and estimated proles of organic matter, poros- ity, saturation fraction, and soil moisture for the case of M = 0:1m. Other prole parameters are set to (A, S w0 (t 1 ), WTD(t 1 ), S w0 (t 2 ), WTD(t 2 )) = (0.8, 0.3, 0.2 m, 0.2, 0.3 m). The shaded regions rep- resent one standard deviation around the mean estimated values for each parameter. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79 4.9 Comparison of true and estimated proles of organic matter, poros- ity, saturation fraction, and soil moisture for the case of M = 0:2m. Other prole parameters are set to (A, S w0 (t 1 ), WTD(t 1 ), S w0 (t 2 ), WTD(t 2 )) = (0.8, 0.3, 0.2 m, 0.2, 0.3 m). The shaded regions rep- resent one standard deviation around the mean estimated values for each parameter. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80 4.10 Comparison of retrieved and in-situ soil dielectric proles at the SoilSCAPE site in Happy Valley, Alaska. . . . . . . . . . . . . . . . 83 4.11 Retrieved soil moisture and OM proles at the SoilSCAPE site in Happy Valley, Alaska. . . . . . . . . . . . . . . . . . . . . . . . . . . 83 5.1 Schematic diagrams of active layer conditions throughout a thaw sea- son. The shaded blue regions represent the part of the active layer that is still frozen. is the seasonal subsidence that occurs when the maximum thaw of the active layer is reached. . . . . . . . . . . . . . 87 5.2 Surface deformation in units of meters with respect to varying water saturation proles (S w0 , WTD) and thaw depths (TD). Parameters of organic matter prole are xed to A = 0.8 and M = 0.15 m. . . . 90 5.3 Surface deformation in units of meters with respect to varying or- ganic matter proles (A, M). . . . . . . . . . . . . . . . . . . . . . . 91 5.4 P-band HH radar backscatter 0 HH in units of dB with respect to varying water saturation proles (S w0 , WTD) and thaw depths (TD). Parameters of organic matter prole are xed to A = 0.8 and M = 0.15 m, and the surface roughness (h) is assumed to be 2 cm in this example. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92 5.5 P-band radar backscatter 0 in units of meters with respect to varying organic matter proles (A, M). (S w0 , WTD, TD, h) is set to (0.5, 0.2 m, 0.5 m, 0.02 m). . . . . . . . . . . . . . . . . . . . . . . . . . 92 5.6 Cost function values of a active layer soil prole with varying water saturation characteristicsS w0 (-) and WTD (m) and thaw depths TD (m). (a) InSAR only (L ). (b) PolSAR only (L ). (c) Joint PolSAR- InSAR (L joint ). The true values of unknowns are set to S w0 = 0.5, WTD = 0.1 m, TD = 0.5 m, h = 0.02 m. . . . . . . . . . . . . . . . 94 5.7 Comparison of true and estimated proles of organic matter, porosity, saturation fraction, and soil moisture. True parameters are set to (A, M S w0 , WTD, TD) = (0.8, 0.2 m, 0.3, 0.2 m, 0.5 m). The shaded regions represent one standard deviation around the mean estimated values for each parameter. . . . . . . . . . . . . . . . . . . . . . . . 96 x 5.8 Comparison of true and estimated thaw depth (TD). Other true pa- rameters are set to (A, M S w0 , WTD) = (0.8, 0.2 m, 0.3, 0.2 m, 0.5 m). The error bars represent one standard deviation around the mean estimated TD values. . . . . . . . . . . . . . . . . . . . . . . . 97 5.9 Jointly retrieved maps of ALT and surface water saturation fraction S w0 using P-band AirMOSS PolSAR data and L-band UAVSAR In- SAR deformation acquired over the Yukon-Kuskokwim (YK) Delta in 2017. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 98 5.10 Jointly retrieved maps of ALT and surface water saturation fraction S w0 using P-band AirMOSS PolSAR data and L-band UAVSAR In- SAR deformation acquired over the Anaktuvuk River in 2017. . . . 99 6.1 Diagram of the joint modeling framework . . . . . . . . . . . . . . . 104 6.2 Estimation of canopy heights for radar-only pixels in Yukon Flats ecoregion in Alaska . . . . . . . . . . . . . . . . . . . . . . . . . . . 105 xi ABSTRACT Radar Remote Sensing of Permafrost Landscapes and Active Layer Soil Properties by Richard H. Chen Chair: Mahta Moghaddam This dissertation investigates radar remote sensing techniques for monitoring per- mafrost landscapes and estimating active layer soil properties, particularly the active layer thickness (ALT) and soil moisture prole, in the Arctic-boreal region. Synthetic aperture radar (SAR) polarimetry has been widely used to estimate soil moisture and freeze/thaw state; this is possible because radar backscatter is sensitive to soil dielec- tric properties, which are strongly correlated with the unfrozen water content and soil composition. Radar retrieval process can be partitioned into 1) dielectric struc- ture estimation based on a radar scattering model and 2) soil dielectric model that translates desired geophysical parameters into dielectric constant. We rst propose a baseline retrieval algorithm using time-series P-band polarimetric SAR (PolSAR) ob- servations to estimate the active layer as a three-layer dielectric structure. To resolve the ambiguity of the retrieved layer thicknesses, an approach of nding the largest possible depths (LPD) is combined with time-series observations, where the ALT is assumed time-invariant between the maximum thaw and early freezing season. In permafrost landscape with extensive organic soils, soil organic matter (SOM) and its degree of decomposition greatly aect soil hydrological and thermal characteristics, xii as well as soil dielectric properties. We thus propose a harmonized parametrization scheme for various soil physical properties, along with a soil dielectric model that can be applied to full range of organic matter content. With this development and a proper prole representation of permafrost soils, we are able to extend the base- line algorithm to estimate ALT and proles of soil moisture and organic matter from time-series P-band PolSAR observations. We also incorporate the seasonal subsidence measured by the interferometric SAR (InSAR) technique into the retrieval process, which allows us further advance the retrieval accuracy and capabilities for permafrost and active layer monitoring. The methods proposed in this dissertation have been applied to the P-band Air- borne Microwave Observatory of Subcanopy and Subsurface (AirMOSS) data ac- quired in 2014 and 2015 under a NASA Interdisciplinary Science (IDS) project, and in 2017 as part of NASA's Arctic-Boreal Vulnerability Experiment (ABoVE) airborne campaigns. The surface deformation derived from the L-band Uninhabited Aerial Ve- hicle Synthetic Aperture Radar (UAVSAR) interferograms is also used in the joint PolSAR-InSAR approach. We have shown, from our retrieval results, that radar re- mote sensing techniques can serve as a tool for monitoring permafrost and active layer dynamics that respond to warming climate and wildre disturbance. xiii CHAPTER I Introduction 1.1 Background and Motivation Climate change has caused air temperature in high-latitude regions to rise at a rate nearly twice faster than the global average [7]. The sustained warming will accelerate the degradation of permafrost and can have strong impacts on ecosystems, hydrologic and landscape processes, and the carbon cycle in the Arctic and Boreal Region (ABR) [8{10]. Permafrost, or perennially cryotic ground, is ground that remains at or below 0 °C for at least two consecutive years (Fig. 1.1). In regions underlain by permafrost, the active layer is the top soil layer subject to seasonal freeze-thaw conditions. Active layer soils have unique hydrological characteristics due to the presence of underlying impermeable permafrost. As active layer thaw deepens in the thaw season, groundwater is accumulated on top of the frost table and forms a water table (upper limit of the saturated zone) within the active layer. The vertical prole of soil moisture (unfrozen water content) is then controlled by both the thawing progress and the water uxes of the active layer [1]. The active layer thickness (ALT) is dened as the annual maximum thaw depth observed at the end of the summer, and serves as an important indicator for monitoring permafrost degradation. Both soil moisture prole and ALT are key variables in understanding the permafrost active layer dynamics that respond to the summer climate and land 1 Figure 1.1: Denition of permafrost and the active layer and the soil features related to seasonal freezing and thawing. Redrawn from Woo 2012 [1]. surface conditions. Ground-based ALT measurements are usually conducted by inserting a metallic rod as a probe into the soil to measure the depth to the frost table [11]. Other approaches include using frost or thaw tubes to measure the ice-water interface and measuring soil temperature prole to determine the depth of the freezing point isotherm [11] [12]. By denition, only when the active layer reaches the maximum thaw in late summer can ALT be measured; otherwise what is measured is referred to as thaw depth. The Circumpolar Active Layer Monitoring (CALM) program has established more than 200 sites in both hemispheres to provide valuable in-situ ALT data and active layer conditions since the 1990s [11]. The CALM ALT data, how- 2 ever, are limited in terms of spatial sampling rate and coverage. The ALT probing measurements at most of the CALM sites are taken in either 1-hectare or 1-km 2 grid (with a probing spacing of 10 m or 100 m, respectively), and the sites are mainly distributed along the Arctic coastline in Alaska and Central Siberia. To account for the spatial variability of active layer dynamics, remote sensing tech- niques are favored for their large spatial coverage as most permafrost areas are remote and inaccessible. Empirical models driven by ground measurements and remotely sensed datasets have been applied to extrapolate the in-situ ALT measurements to generate regional maps of ALT and near-surface permafrost distribution [13{16]. De- tailed soil process models that account for the physics of coupled heat and moisture transfer are used to further improve the modeling accuracy and the prediction of active layer dynamics [17] [18]. Radar interferometry has also been used to estimate ALT by correlating the observed interferometric phase, hence surface deformation, to the seasonal ground subsidence caused by the thawing of the active layer soils [19] [6]. Synthetic aperture radar (SAR) polarimetry has demonstrated the capability of retrieving surface and subsurface prole soil moisture from radar backscattering mea- surements [20{22]. This is possible because the measured backscattering coecients are directly related to the dielectric properties of imaged scenes, and soil dielectric constant is directly related to soil moisture. Therefore, soil moisture retrievals are equivalent to retrieving the dielectric constant (or relative permittivity) of soil. Soil dielectric constant is also a function of other soil physical properties including soil tex- ture, bulk density, temperature, and freeze-thaw conditions [4] [23]. So, soil dielectric models and knowledge of other soil parameters are required for the conversion of di- electric constant to soil moisture, which can take place either during or after the radar scattering inversion (or retrieval) process. NASAs Airborne Microwave Observatory of Subcanopy and Subsurface (AirMOSS), an Earth Venture 1 (EV-1) mission from 2012 to 2015, has proven its capability of using P-band SAR to map root-zone soil 3 Figure 1.2: AirMOSS radar swaths and permafrost observation sites over Northern Alaska. CALM: Circumpolar Active Layer Monitoring program [www2.gwu.edu/calm]. UAF GIPL: Geophysical Institute Permafrost Laboratory, University of Alaska Fairbanks [permafrost.gi.alaska.edu]. moisture (RZSM) proles at various major North American biomes with good delity [16]. As a low-frequency radar, the AirMOSS P-band (420-440 MHz) radar is also expected to be able to penetrate into the active layer and reach the permafrost ta- ble. Since permafrost has a much lower dielectric constant compared to thawed soils, the dielectric contrast around permafrost table is sharp enough to result in a strong scattering of electromagnetic waves, providing a potential signature in backscattering measurements for directly retrieving ALT. In support of a NASA Interdisciplinary Science (IDS) project, the AirMOSS P- band radar ew and acquired ve time-series SAR measurements over a 2,500 km transect in northern Alaska since August 2014, covering a number of existing per- mafrost observation sites, as shown in Fig. 1.2. The acquisitions of AirMOSS ights are designed to capture three active layer conditions, including maximum thaw (mid- late August), partially frozen (early October), and fully frozen (early-mid April). The 4 Figure 1.3: Flight transects own during the 2017 ABoVE airborne campaigns. AirMOSS fully polarimetric SAR imagery are provided at 0.5 arcsec and 3 arcsec post- ings, and are radiometrically calibrated to 0.5 dB [24]. As a low-frequency radar, the AirMOSS P-band (420-440 MHz) radar is expected to be able to penetrate into the active layer and reach the permafrost table, providing an opportunity for measuring ALT directly. This AirMOSS dataset provides the rst-ever P-band SAR imagery across various permafrost zones in Alaska, from which the majority of developments in this thesis was inspired and motivated. To better understand the vulnerability and resilience of socio-ecological systems in the ABR to environmental and climate change, the NASA Terrestrial Ecology Program has initiated the Arctic-Boreal Vulnerability Experiment (ABoVE), a large- scale eld campaign currently conducting in Alaska and Western Canada, for 8 to 10 years, starting in 2015. Multiple remote sensing sensors, including airborne SAR sys- tems (P-band AirMOSS and L-band Uninhabited Aerial Vehicle Synthetic Aperture Radar (UAVSAR)), waveform lidar, and hyperspectral imaging systems, have own and acquired data over the ABoVE study domain during the airborne campaigns in 2017, as shown in Fig. 1.3. The ight transects span over major vegetation (tun- 5 Figure 1.4: Retrieval process of subsurface geophysical properties. dra to boreal forest), permafrost (continuous to sporadic) and climate (polar Arctic to boreal) gradients, and also encompass a diverse range of disturbance conditions, including thermokarst, roads and other human infrastructure development, and wild- re burns. These datasets provide a valuable opportunity for the fusion of dierent physics which could potentially enable a more accurate and integrated framework for monitoring permafrost active layer dynamics. 1.2 Retrieval Process of Subsurface Geophysical Properties The retrieval process employed in this thesis is forward-model-centric, in which optimization is carried out to search for the solutions that minimize the dierence between radar observations and forward model simulations. The overall retrieval pro- cess is illustrated in Fig. 1.4, where both soil dielectric and radar scattering models are considered as part of the forward model. Soil dielectric model is used to calcu- late the dielectric constants of soils with certain geophysical properties (e.g., soil moisture, freeze/thaw state) and to model the entire soil column as a layered dielec- 6 tric structure. Radar scattering model is used to simulate the radar backscattering coecients 0 from this dielectric structure. The unknown geophysical variables X can be estimated by a global optimization scheme which iteratively adjusts X for the forward model simulations, until a specically dened cost function is minimized. 1.3 Objectives As more remote sensing sensors and techniques have been deployed to observe the Arctic-boreal landscapes, the objectives of this research are not only to develop and formulate newer radar remote sensing techniques for retrievals of active layer properties, but also to link the physics of various sensor measurements and soil process models for better monitoring of permafrost active layer dynamics. The conventional methods for retrievals of surface/subsurface soil features using radar usually solve the problems from the aspects of electromagnetics, and do not consider other physics that also govern the dynamics in the soil, such as hydrology and thermal dynamics. That being said, the objectives of this thesis can be categorized into the Electromagnetic Inverse Modeling and Joint-Physics Approach: For Electromagnetic Inverse Modeling, the objectives are: 1. Identify the dielectric properties of permafrost soils as well as the subsurface features that can be retrieved from radar observations. 2. Explore the potentials of time-series and multi-frequency radar observations that are available in the ABR. 3. Develop a baseline radar retrieval algorithm which can estimate the dielectric structure of permafrost soils. For Joint-Physics Approach, the objectives are: 7 1. Develop a parametrization scheme that harmonizes the dielectric, hydraulic, and thermal properties of permafrost soils. 2. Generalize the baseline algorithm using the harmonized soil parametrization to improve the retrieval results that can have more direct impacts on other science disciplines. 3. Advance the polarimetric SAR (PolSAR)-based approach to be combined with other remote sensing techniques, as well as the modeling of soil hydrologic pro- cesses. 1.4 Contributions 1. First-ever demonstration of direct retrieval of active layer properties, including ALT and soil moisture proles, in permafrost landscapes using time-series low- frequency SAR observations 2. Quantication of radar sensing depth that accounts for both the attenuation and scattering eects caused by the subsurface dielectric discontinuities and the calibration accuracy of the radar system. 3. Demonstration of harmonized parametrization of soil physical properties that is applicable for full range of organic matter content and degree of decomposition in permafrost soils. 4. Improved soil prole representation of active layer that enables simultaneous estimation of organic matter and soil moisture proles from radar observations. 5. Advanced ALT and soil moisture estimation by incorporating the surface de- formation derived from interferometric SAR (InSAR) into the PolSAR-based approach. 8 6. Enhanced understanding and monitoring of permafrost and active layer dynam- ics that respond to climate warming and wildre disturbance. 9 CHAPTER II Retrieval of Permafrost Active Layer Properties: A Time-Series Approach 2.1 Introduction The simultaneous retrieval of multiple geophysical parameters can be dicult be- cause the combined eects of these parameters are often ambiguous in radar backscat- tering measurements. For example, a smooth surface, dry soil moisture, or frozen soils can all contribute to a low value of backscattering coecient and their eects are therefore dicult to separate in the inversion process. To resolve the ambiguity among the unknown parameters, time-series radar observations can be used if some of the parameters do not change signicantly with time and can be considered con- stant within the time period of interest [25]. A typical example of such time-invariant parameter is surface roughness, with the assumption that there are no major human activities or land surface processes within the time window considered [26]. In the context of ALT retrievals, there are also some physical properties of active layer soil that can be considered as nearly time-invariant between the maximum thaw in late August and the partially frozen condition in early October. The retrieval algorithm developed in this work aims at retrieving the dielectric prole of the active layer using time-series P-band SAR observations. The time-series 10 utilized in the retrieval includes any observation acquired after the maximum thaw and before the bottom frost table starts to rise signicantly so that the ALT would be observable for the entire time window. Moreover, the soil dielectric constant in the saturated zone of the active layer is also assumed constant in time, because the water content depends only on the porosity of the host medium. From the retrieved dielectric prole, the active layer properties, such as ALT, soil moisture, and freeze- thaw conditions, can be derived. 2.2 In-Situ Data and Time-Series Considerations To observe seasonal active layer dynamics, in-situ soil moisture and temperature proles taken at the Global Terrestrial Network for Permafrost (GTN-P) site [2] operated by the Geophysical Institute Permafrost Laboratory (GIPL), University of Alaska Fairbanks (UAF) in Deadhorse, Alaska are shown in Fig. 2.1. At the onset of the thaw season, the active layer starts to thaw from the surface and the soil moisture increases rapidly as the soil water content transitions from the frozen to the thawed state. Most of the thawed active layer soils remain saturated until the beginning of the freezing season, making the corresponding soil moisture variations nearly constant except that the surface layers might occasionally be under unsaturated conditions. ALT is observed by the thaw depth at the maximum thaw in late August. The freeze-up of the active layer begins after the time of the maximum thaw when the ground surface temperature drops below 2 °C but is still above 0 °C [27] [28]. The ground freezing occurs from both sides of the active layer, with the upward freezing starting rst from the permafrost table, followed by the downward freezing from the surface when the ground surface temperature becomes negative. When the freezing of active layer soils starts, the soil moisture drops from the values of saturation due to the decreasing unfrozen water content. The soil temperature does not decrease signicantly below 0 °C until the freezing of that soil layer is mostly completed and 11 May Jun Jul Aug Sep Oct Nov Dec Jan 0 0.1 0.2 0.3 0.4 0.5 Soil Moisture [m 3 /m 3 ] 15 cm 30 cm 60 cm May Jun Jul Aug Sep Oct Nov Dec Jan 2015 / 2016 -10 -5 0 5 Soil Temperature [ ° C] (a) May Jun Jul Aug Sep Oct Nov Dec Jan 2015 / 2016 0 0.2 0.4 0.6 0.8 Depth [m] -10 -10 -8 -8 -8 -8 -6 -6 -6 -6 -4 -4 -4 -4 -4 -2 -2 -2 -2 -2 -2 -2 -2 -2 -2 -1 -1 -1 -1 -1 -1 -0.5 -0.5 -0.5 -0.5 -0.5 -0.5 -0.5 -0.5 -0.3 -0.3 -0.3 -0.3 -0.3 -0.3 -0.3 -0.3 -0.1 -0.1 -0.1 -0.1 -0.1 0 0 0 0 0 0 0 2 2 2 2 2 2 2 2 2 2 4 4 4 4 4 4 4 6 6 6 6 6 6 8 8 8 (b) Figure 2.1: In-situ data from the GTN-P site operated by UAF GIPL in Deadhorse, Alaska [2]. (a) Soil moisture and temperature proles taken by Stevens HydraProbe sensors at depths of 15, 30, and 60 cm. The data are extracted from the real-time data portal at http://permafrost.gi.alaska.edu/site/dh1. (b) Isotherms interpolated from the thermistor measurements at depths of 0.01, 0.03, 0.11, 0.18, 0.25, 0.34, 0.41, 0.57, 0.72, and 0.93 m. The vertical dashed lines mark the AirMOSS ights on August 29 and October 1 in 2015. the unfrozen water content has gone below a certain level [29]. For example, the soil moisture values at 60 cm starts to decrease from around 0.4 m 3 /m 3 at -0.05 °C, representing the start of freezing. As the freezing proceeds progressively, the soil moisture drops to 0.14 m 3 /m 3 at -0.3 °C, from which the soil temperature starts to decrease signicantly. The soil layer is not fully frozen until the soil moisture reaches 12 the minimum value around 0.05 m 3 /m 3 at -6 °C. As suggested in [30], a threshold of -0.3 °C is chosen to indicate when the majority of liquid water has been frozen for a soil layer, so the -0.3 °C isotherm in Fig. 2.1(b) can be seen as the the depths to the two freezing fronts. The upward freezing from the permafrost table is seen to start at the beginning of September, and about 3.5 cm of the bottom active layer has frozen upwards by early October. The downward freezing front is observed to deepen at a faster rate after mid September and the active layer becomes entirely frozen around mid November. There is a time period when the upward freezing front has not risen signicantly from the permafrost table and the active layer is not entirely frozen. During such time period the depth to the bottom frost table is still close to the maximum thaw depth, namely the ALT (e.g., Fig. S2 in [31]); therefore, ALT can still be observable within this time window. 2.3 Radar Forward Scattering Model Heterogeneities of subsurface soil properties including soil moisture proles, vege- tation roots, small pebbles and air voids, are known to have signicant eects on radar backscatter if radar signals can reach the depths where these heterogeneities reside and have features that are a signicant fraction of the radar wavelength. When using low-frequency radars such as P-band (430 MHz center frequency), soil moisture pro- les and any layer boundary contrast in dielectric constant (water tables, freeze/thaw boundaries) will be the dominant features observable in backscattering measurements. The subsurface soils, therefore, should be represented by a layered dielectric structure with a top rough surface and homogeneous dielectric layers. The dielectric constant and thickness of these layers construct a dielectric prole for describing the soil condi- tions throughout the soil column. The P-band AirMOSS mission 1 has already proven 1 The AirMOSS mission here is referred to the NASA Earth Ventures 1 (EV-1) mission from 2012 to 2015; the AirMOSS mentioned elsewhere in this paper is referred to the P-band radar instrument. 13 Figure 2.2: Geometry of a multilayer dielectric structure (left) and the three-layer model (right) used in the retrievals of active layer properties. The multilayer struc- ture shows a more general case of N L layers that can be used to represent more sophisticated stratigraphic features in the active layer. The three-layer model in- cludes the surface active layer ( 1 ), the bottom saturated active layer ( 2 ), and the permafrost ( f ). The top layer dielectric constant 1 and thickness z 1 can vary be- tween mid-late August and early October to account for the changes in soil moisture or freeze-thaw state in the near-surface layer. its capability of mapping root-zone soil moisture (RZSM) proles at various major North American biomes based on these assumptions [22]. The layered dielectric struc- ture is also applicable for modeling the active layer and underlying permafrost, where both soil moisture and freeze/thaw state can be represented by proper dielectric con- stant values. The geometry of a multilayer dielectric structure is illustrated in Fig. 2.2, where each layer i has a complex dielectric constant i = 0 i +i 00 i and is separated by an interface at z = z i . The top rough surface (z = 0) separates the free space (layer 0, 0 ) from the subsurface soil layers and has roughness represented by root-mean- square (rms) heighth and correlation lengthl c . Both soil composition (sand/silt/clay fractions and bulk density) and geophysical conditions (water content and freeze/thaw state) can aect soil dielectric constant. In this work, we assume that the real part of the dielectric constant of active layer soils is relatively high ( 0 > 10) when thawed 14 and very low ( 0 5) when frozen. This assumption is based on the fact that liquid water content is the major contributor to the values of soil dielectric constant and the signicant dierence between the dielectric constants of liquid water and ice, which are approximately 80 and 3.2 in microwave region, respectively [32]. The dielectric contrast between thawed and frozen soils results in a strong scattering of electromagnetic waves and therefore provides a potential signature in backscattering measurements for freeze/thaw boundary detection. The imaginary part 00 of a complex dielectric constant, or the dielectric loss fac- tor, is associated with the power loss of electromagnetic waves propagating in a lossy medium. Both 0 and 00 of soil dielectric constant are required to accurately sim- ulate the scattering process within the soil medium, and their relationship can be described by soil dielectric models along with associated soil parameters from exist- ing soil databases. In the Arctic regions, however, the soil texture data is usually limited by the available detailed soil surveys and has coarse spatial distributions of soil properties. Moreover, the soils of the Arctic usually have high soil organic matter (SOM) content, for which both a high-quality soil organic carbon map and a soil dielectric model accounting for a wide range of organic matter content are needed. To observe the soil dielectric properties for the Arctic soils, we extracted all ve soil classes on the Alaska North Slope from the soil texture data by De Lannoy et al. [3], where four mineral soil classes are classied with low (SOM < 0.49 %) to moderate (1.1 % SOM < 15 %) organic matter content and one soil class for peat soils with high organic matter content (SOM 15 %). Fig. 2.3 shows the dielectric properties of the ve soil classes on the Alaska North Slope using the dielectric models by Peplin- ski et al. [4] and by Bircher et al. [5] for mineral and organic soil classes, respectively. The ratio of 00 to 0 ranges from 0.04 to 0.23 over the entire volumetric soil moisture range. The vertical distribution of SOM content throughout the soil column is also required to determine a proper soil dielectric model to be used for each soil layer, and 15 0 0.2 0.4 0.6 0.8 1 m v [m 3 /m 3 ] 0 20 40 60 80 Soil Class #68 Soil Class #75 Soil Class #224 Soil Class #245 Organic Soil 0 0.2 0.4 0.6 0.8 1 m v [m 3 /m 3 ] 0 2 4 6 8 Soil Class #68 Soil Class #75 Soil Class #224 Soil Class #245 Organic Soil (a) 0 0.2 0.4 0.6 0.8 1 m v [m 3 /m 3 ] 0 0.05 0.1 0.15 0.2 0.25 / Soil Class #68 Soil Class #75 Soil Class #224 Soil Class #245 Organic Soil (b) Figure 2.3: Dielectric properties of major soil classes on the Alaska North Slope. (a) Real ( 0 ) and imaginary ( 00 ) part of soil dielectric constant and (b) the ratio of 00 to 0 over the entire volumetric soil moisture (m v ) range based on the soil porosity, which is equivalent to the soil moisture at saturation. The soil classes are extracted from the soil texture data by De Lannoy et al. [3] and there are ve major soil classes on the Alaska North Slope (#68, #75, #224, #245 for the mineral soil classes with low (SOM< 0.49 %) to moderate (1.1 % SOM< 15 %) SOM content, and #253 for the organic soil class with high SOM content (SOM 15 %)). The soil dielectric models by Peplinski et al. [4] and by Bircher et al. [5] are used to compute the dielectric constant for the mineral and organic soil classes, respectively. such information is not available in existing soil databases. Since we currently do not have a good soil texture and SOM data suitable for the needs of this study, we simply assume that 00 = 0 = 0:15 for the soils discussed in this work. This assumption can be improved if more detailed soil texture and soil organic carbon maps and accurate soil dielectric models for the Arctic soils become available in the future. To calculate radar backscattering coecients of a layered dielectric structure, a small perturbation method (SPM) up to the rst order developed in [33] is used as the forward model in this work. Fig. 2.4 shows a simulation example for co- polarized P-band backscattering coecients ( 0 hh and 0 vv ) from a two-layer dielectric structure with a varying top-layer depth (z 1 ). The top layer has a dielectric constant of 0 1 = 30, representing the thawed active layer with high unfrozen water content; the bottom layer represents the frozen active layer or permafrost with a low dielectric 16 0 0.2 0.4 0.6 0.8 z 1 [m] -28 -26 -24 -22 -20 -18 -16 -14 -12 -10 -8 0 [dB] 0 hh 0 vv Figure 2.4: Simulated radar backscattering coecient of a two-layer dielectric struc- ture with varying rst-layer interface (z 1 ). The real parts of layer dielectric constants are set to be ( 0 1 , 0 2 ) = (30, 5) in this example. constant of 0 2 = 5. The depth z 1 is therefore the thaw depth of the active layer, or ALT if the maximum thaw is reached. The scattered wave within the top layer can interfere constructively and destructively due to multiple scattering, causing 0 hh and 0 vv to behave as an oscillating function of the thaw depth z 1 , with a peak-to- peak distance approximately equivalent to half of the eective wavelength ( eff;1 =2 = 0 =2 Re( p 1 )). Both soil layers are considered as a lossy medium with a positive imaginary part ( 00 1 = 4:5 and 00 2 = 0:75) of dielectric constant. As a result, the amplitude of 0 variations decays asz 1 moves toward larger depths, and 0 converges to the value expected for a half-space structure with the top-layer dielectric constant 1 . The oscillating behavior of the backscattering coecient 0 versus thaw depth (z 1 ) poses challenges for thaw depth or ALT retrievals, as the retrieved thaw depth values can be ambiguous for a given value of backscattering coecient. 17 2.4 Radar Sensing Depth Conventionally, when determining how deep the radar signal can reach and \sense" the subsurface soil properties, penetration depth p is the most commonly used metric, dened as the depth at which the power of radar signal attenuates to 1=e of its original level right beneath the surface. However, penetration depth only accounts for the attenuation as the electromagnetic wave travels within the lossy medium, but does not account for the scattering eects caused by dielectric discontinuities. It also does not account for radar system sensitivity characteristics, or for radar calibration uncertainties. If the soil column is a homogeneous half space or has very weak dielectric discontinuities, the scattering within the subsurface will not be strong enough to contribute to the total backscattering coecient observed by the radar. In this case, the radar backscatter will only carry the information from the top surface, regardless of the penetration depth. On the other hand, if there is a high dielectric contrast deep in the subsurface, the scattering due to that contrast can still have observable impacts on the total backscattering coecient even if a relatively high signal attenuation is experienced in the medium, in which case the radar can therefore potentially \sense" deeper than the nominal penetration depth. Considering a soil dielectric prole ((z) = p (z) for z < z s ) on top of a target dielectric substrate ((z) = s forzz s ), we varyz s and calculate the corresponding backscattering coecient 0 (z s ). To determine whether 0 is still sensitive to the presence of s , we dene a \radar sensing depth" s as s = maxfz s jD + max (z s )> 2 n g; (2.1) where D + max (z s ) = maxD + (z s ) (2.2) 18 0 0.2 0.4 0.6 0.8 z 1 [m] 0 50 100 150 D + max 0.4 0.5 0.6 0.7 0.8 0 0.5 1 1.5 D + max D + max ( s ) n 2 Figure 2.5: Example of calculating radar sensing depth s for a top dielectric layer of 0 1 = 30 with a target dielectric substrate 0 2 = 5 at z 1 . The inset shows a zoom-in view of how s is determined. n represents the radar calibration accuracy and is set to 0.5 dB in this case. and D + (z s ) = 1 2 X pp=hh;vv 0 pp (z s ) 0 pp (z + s ) 2 : (2.3) The radar calibration accuracy, denoted as n , represents the standard deviation of residual radar backscattering cross section measured from known calibration targets such as corner re ectors [24]. n can therefore be considered as the standard deviation of the measured radar backscatter 0 . D + (z s ) is a set of average L 2 -norm type distances between the co-polarized values of 0 hh (z s ) and 0 vv (z s ), which determines the variation of 0 forz + s =fzjz >z s g. If the maximum ofD + (z s ) is larger than 2 n , the presence of s at z =z s is considered detectable because one can distinguish this depth from another larger depth by the dierence of backscattering measurements. In other words, if the variation of 0 forz >z s is smaller than n , then it is dicult for the radar to accurately sense the target substrate because that variation can be due to noise or uncertainties introduced during measurements or post-processing. Fig. 2.5 shows an example of calculating s from the simulated backscattering coecient shown 19 2 10 20 30 40 50 60 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 Depth [m] s (P-band) s (L-band) p (P-band) p (L-band) Figure 2.6: P- and L-band radar sensing depth s and penetration depth p with varying top-layer dielectric constant 0 1 . The calibration accuracy n of the radar systems is assumed to be 0.5 dB. The sharp peak at 0 1 = 5 is because of the perfect match between the two layers and therefore lack of a scattering interface. Note that the sensing depth is considerably larger than the penetration depth for both L-band and P-band. in Fig. 2.4, where the top dielectric prole is considered constant as 0 p (z) = 0 1 = 30 with a target dielectric substrate 0 s = 0 2 = 5. The variations of the sensing depth s and the penetration depth p with respect to 0 1 are shown in Fig. 2.6 for P-band (430 MHz) and L-band (1.26 GHz) radars with an expected n of 0.5 dB, and a two-layer dielectric structure with a dielectric substrate 0 2 = 5. When 0 1 is close to 0 2 , representing a low dielectric contrast, s is much shallower than when 0 1 and 0 2 are well separated. Once 0 1 gets suciently higher than 0 2 , s decreases as the attenuation becomes higher due to the increasing imaginary part 00 1 of the top layer dielectric constant. Since p only accounts for the wave attenuation when traveling within 0 1 , it is a monotonically decreasing function of 0 1 and it signicantly underestimates the depths that radar backscatter can actually detect, as compared to the values of s . As expected, the lower frequency P-band radar has larger s and p than L-band radar, and is capable of detecting a low-dielectric 20 2 10 20 30 40 50 60 1 0.04 0.2 0.4 0.6 0.8 z 1 [m] -30 -25 -20 -15 -10 -5 (a) 0 hh 2 10 20 30 40 50 60 1 0.04 0.2 0.4 0.6 0.8 z 1 [m] -30 -25 -20 -15 -10 -5 (b) 0 vv Figure 2.7: P-band radar backscattering coecients 0 as a function of top-layer dielectric constant 0 1 and depth of the top layer z 1 . Dashed-dot line: corresponding radar sensing depth s . substrate at depths larger than 60 cm, depending on the value of 0 1 . 2.5 Sensitivity of Radar Backscatter to Subsurface Variables Consider a two-layer geometry with a top rough surface (h;l c ) = (0:02 m; 0:2 m); there will be three unknown variables ( 0 1 , 0 2 , z 1 ) to be determined. To show the sensitivity of P-band radar backscatter to these three unknowns, we x either 0 1 or 0 2 at a time, and vary the other two unknowns when calculating 0 hh and 0 vv . The xed values used in the sensitivity studies are 0 1 = 30 (thawed and wet active layer soils) and 0 2 = 5 (frozen soil or permafrost), and the simulated 0 hh and 0 vv are shown in Fig. 2.7 and Fig. 2.8. When the bottom layer is assumed to be a frozen layer ( 0 2 = 5), the sensitivity of 0 to z 1 is high as long as z 1 is varying within the sensing depth s , as shown in Fig. 2.7. The values of 0 become less sensitive to, or saturated as a function of, 0 1 as 0 1 increases, resulting in larger uncertainty when retrieving high dielectric constants. The same challenge is encountered for radar retrievals of soil moisture when high surface soil moisture is present [25]. Fig. 2.8 shows that 0 is relatively insensitive to the bottom dielectric constant 0 2 21 2 10 20 30 40 50 60 2 0.04 0.2 0.4 0.6 0.8 z 1 [m] -25 -20 -15 -10 (a) 0 hh 2 10 20 30 40 50 60 2 0.04 0.2 0.4 0.6 0.8 z 1 [m] -25 -20 -15 -10 (b) 0 vv Figure 2.8: P-band radar backscattering coecients 0 as a function of bottom-layer dielectric constant 0 2 and depth of the top layer z 1 . Dashed-dot line: corresponding radar sensing depth s . compared to the other two variables. This is due to the fact that when electromagnetic waves travel into the bottom layer, there are no other dielectric discontinuities that can scatter the wave back to the upper layers and therefore no multiple scattering takes place in the bottom layer. The main function of 0 2 is to form the dielectric contrast with 0 1 atz 1 ; the variations of 0 are nearly symmetrically opposite when 0 2 is higher vs. lower than 0 1 . Moreover, since the dielectric constant for frozen soils is in the range of 2 to 8, and the 0 variations with respect to 0 2 are negligible within this range, we can assume 0 = 5 for the frozen dielectric substrate (frozen active layer or permafrost) at all times without introducing much error in retrievals. By xing the dielectric constant for the bottom layer, the number of unknowns can be reduced and the retrievals become less ambiguous. 2.6 Inverse Problem To retrieve the unknown parameters of a layered dielectric structure from radar backscattering measurements, we need to form and solve an inverse problem to nd the optimal (or the most probable) solution set for the unknowns. The retrieval 22 performance depends on how well the forward model (including assumptions about the structure and its parameters) can represent the actual scenes imaged by radar, as well as on calibration and processing errors. Ideally, one can retrieve all of the unknown subsurface variables, up to the radar sensing depth, from 0 measurements if there is sucient number of measurements (i.e., equal to, or greater than, the number of unknowns) [34]. In practice, however, the number of measurement points is usually far less than the number of unknowns. Therefore, assumptions based on prior information need to be made in order to limit the complexity of the problem and mitigate the errors from an underdetermined system. The main challenge of retrieving the depth z s at which a dielectric substrate s is located is the ambiguity of 0 with respect to depth variations. Even if the dielectric prole p (z) on top of s is perfectly retrieved, 0 is still an oscillating and decaying function ofz s , as shown in Fig. 2.4. For any given 0 hh and 0 vv , there will be multiple possible solutions forz s . Unlike the techniques applied to ground-penetrating radars (GPR), where the time history of returned signal (both amplitude and phase information) is available, here we use only the amplitude of the polarimetric SAR observations for estimating the dielectric properties and layer depths. The use of co-pol phase dierence could potentially result in more accurate retrievals, however the current phase calibration of AirMOSS image data is unknown and therefore it is not possible to use the co-pol phase. 2.6.1 Cost Function In forward-model-centric inversion schemes, optimization is carried out to search for the solutions that minimize the dierence between observations and forward model 23 predictions. To describe this dierence, we dene the following cost function: L(X) = 1 N X pq;f; inc 0 pq (X;f; inc )d pq (f; inc ) 2 ; (2.4) where N is the number of measurement points and X = ( 0 i ;z i ) represents the un- known variables of interest, including soil layer dielectric constants 0 i and layer depths z i . The values of 0 pq (X;f; inc ) and d pq (f; inc ) are the calculated and measured pq- polarized backscattering cross section at radar frequency f and local incidence angle inc , respectively, where p and q2fv;hg. Since the calibration accuracy of cross- polarized radar backscatter is not yet known for the AirMOSS radar system due to cross-pol contamination [35], only co-polarized radar backscatter, 0 hh and 0 vv , are considered throughout this paper. The number of local minima and the shape of a cost function can greatly impact the performance of retrieval algorithms. For surface soil moisture retrievals, since 0 is a monotonic function with respect to surface roughness h and surface dielectric constant 0 1 , the cost function has a unique global minimum. Fig. 2.9 shows the cost functions of a half-space problem (N L = 1) with dierent surface dielectric constants 0 1 . For this example, the surface roughness is assumed known and xed. Although each cost function in Fig. 2.9 has only one global minimum, the uncertainty of the retrieved 0 1 is higher for larger 0 1 , as represented by the wider range of 0 1 under which the cost function values are below a nominal error criterion set to 2 n . Global or local optimization algorithms can be applied to nd the global minimum of the cost function, which would be the optimum solution for the problems in Fig. 2.9. For the simplest case of ALT retrievals, we consider a two-layer problem with known dielectric constants 0 1 and 0 2 representing the active layer and permafrost, respectively, and one unknown depth z 1 representing ALT. This is just the starting point; we will make the problem more realistic in the next steps. Since 0 is an os- 24 2 10 20 30 40 50 60 1 0 0.2 0.4 0.6 0.8 1 L(X) 1 true = 10 1 true = 20 1 true = 40 n 2 Figure 2.9: Cost function values of one-layer (half-space) problems with dierent rst-layer dielectric constants 0 1 . The dashed line is the error criterion nominally set to 2 n , the variance of calibrated radar backscattering cross section. Vertical axis is truncated at L(X) = 1 for display simplicity. cillating and decaying function (Fig. 2.4), the corresponding cost function will have multiple local minima, as shown in Fig. 2.10 for the noise-free cases. Instead of having possible solutions spreading around the true point as in the case of surface soil moisture retrievals, the possible solutions to z 1 are discrete groups around the ambiguous points. Due to the presence of noise and calibration errors in the mea- sured 0 , the global minimum could occur at one of local minima which have cost function values smaller than 2 n . The solution of the global minimum will therefore not necessarily correspond to the true solution z true 1 , but rather z true 1 could occur at one of the local minima. Regularization is required for the cost function in (2.4) to have a better estimation of z true 1 . 2.6.2 Largest Possible Depths The cost function shown in Fig. 2.10(a) is one of the worst cases for retrieval because the values of 0 (z true 1 ) are close to the ones expected for a half space (N L = 1, 25 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 z 1 [m] 0 0.25 0.5 0.75 1 L(X) z 1 true = 0.40 m n 2 (a) 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 z 1 [m] 0 0.25 0.5 0.75 1 L(X) z 1 true = 0.32 m n 2 (b) Figure 2.10: Cost function values for a two-layer problem with ( 0 1 , 0 2 ) = (30, 5) and one unknown z 1 . (a) z true 1 = 0.4 m. (b) z true 1 = 0.32 m. In (b), red circles and black crosses are the local minima that occur at z 1 z true 1 and at z 1 > z true 1 , respectively. Note that there are two red circles located closely aroundz true 1 = 0.32 m. The dashed line is the error criterion nominally set to 2 n , the variance of calibrated radar backscattering cross section. 0 1 = 30), meaning that there is a large number of local minima with cost function values close to zero. The similarity between the local minima makes it dicult to distinguish which one is more likely to be located around the desired true pointz true 1 . For the cost function shown in Fig. 2.10(b),z true 1 occurs near one of the peaks of 0 so there is noz 1 that is larger thanz true 1 and has 0 values similar to 0 (z true 1 ). Therefore, the local minima that occur at z 1 > z true 1 will result in cost function values larger than L(z true 1 ) and they correspond to the peaks of 0 after 0 (z true 1 ) with distances approximately equivalent to eff;1 =2. The local minima that occur at z 1 < z true 1 still have similar cost function values as L(z true 1 ) due to their 0 values being similar to 0 (z true 1 ). Their distances, however, are all smaller than eff;1 =2 and are either larger or smaller than eff;1 =4. We then can separate the local minima that occur 26 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 z 1 [m] -28 -26 -24 -22 -20 -18 -16 -14 -12 -10 -8 0 [dB] 0 hh (ambiguous) 0 hh (solvable) 0 vv (ambiguous) 0 vv (solvable) (a) 0 5 10 15 20 n 0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 z min [n] eff /4 (b) Figure 2.11: (a) Simulated radar backscattering coecient of a two-layer dielectric structure with varying rst-layer interface z 1 and ( 0 1 , 0 2 ) = (30, 5). The regions where z 1 can be solved by the LPD approach are highlighted. (b) z min [n], the distances between adjacent local minima of L(z 1 ) in Fig. 2.10(b). Red circles and black crosses correspond to the local minima in Fig. 2.10(b) that occur at z 1 z true 1 and at z 1 >z true 1 , respectively. Red star corresponds to n q in (2.5). at z 1 z true 1 from the ones that occur at z 1 > z true 1 , and choose the largest z 1 in the former group as the \largest possible depth (LPD)" to the unknown z 1 for the measured 0 (z true 1 ). If the measured 0 values are located around one of the extrema in the curves of 0 (z 1 ) and are still higher (or lower) than the next maximum (or minimum), the corresponding z 1 's are considered solvable by the LPD approach, as shown in Fig. 2.11(a). For the z 1 's located outside the solvable regions, the LPD solution will result in overestimation of z true 1 due to the ambiguity of the measured 0 values. The LPD approach starts with the cost function curve L(z 1 ), which is computed at a ne step length set to 0.05 cm to prevent numerical errors, as shown in Fig. 2.10(b). We then locate the local minima of L(z 1 ) and record the corresponding z 1 's in a sequence denoted as z min [n]. The distances between adjacent local minima, z min [n] =z min [n+1]z min [n], are plotted in Fig. 2.11(b) as an example. The local- minimum pairs with z min < eff;1 =4 corresponds to the points centered around one 27 0 0.2 0.4 0.6 0.8 z 1 true [m] 0 0.2 0.4 0.6 0.8 1 z 1 est [m] LPD s (a) n = 0.5 dB 0 0.2 0.4 0.6 0.8 z 1 true [m] 0 0.2 0.4 0.6 0.8 1 z 1 est [m] min L(X) s (b) n = 0.5 dB 0 0.2 0.4 0.6 0.8 z 1 true [m] 0 0.2 0.4 0.6 0.8 1 z 1 est [m] LPD s (c) n = 0.3 dB 0 0.2 0.4 0.6 0.8 z 1 true [m] 0 0.2 0.4 0.6 0.8 1 z 1 est [m] LPD s (d) n = 0.1 dB Figure 2.12: Comparison of error performance of z est 1 . (a) LPD approach at n = 0.5 dB. (b) Global minimum approach at n = 0.5 dB. (c) LPD approach at n = 0.3 dB. (d) LPD approach at n = 0.1 dB. The gray shaded regions represent one standard deviation around the average values of z est 1 . of the extrema of 0 (z 1 ), and hence a distance shorter than a quarter wavelength. The sequence index of the last local-minimum pair with z min < eff;1 =4 can be written as n q = maxfnj z min [n]< eff;1 =4g: (2.5) z min [n q ] and z min [n q + 1] are the z 1 's of the pair z min [n q ] and they are centered 28 around the extremum of 0 (z 1 ) at z q = arg max z min [nq ]z 1 z min [nq +1] L(z 1 ): (2.6) The next extremum of 0 (z 1 ) is located at z h = z min [n q + 2] and the LPD can be determined by z LPD = 8 > > < > > : z min [n q + 1]; if L(z q )L(z h ) z min [n q + 2]; if L(z q )>L(z h ): (2.7) The LPD approach is essentially looking for the extremum of the 0 (z 1 ) curve that is the closest to the measured 0 values and determines z LPD based on that. The resolution of the LPD approach will therefore be approximately eff;1 =4, which is the distance between the adjacent extrema of the 0 (z 1 ) curve. We evaluate the LPD approach with Monte-Carlo simulations for dierent levels of radar calibration accuracy n . The simulations consider a two-layer dielectric structure of ( 0 1 , 0 2 ) = (30, 5) with the z true 1 varying from 0.04 m to 0.8 m. The synthesized radar data are generated from SPM simulations and an independent zero- mean Gaussian noise with standard deviation of n (0.5, 0.3, and 0.1 dB) is added to each polarization channel. The results of the LPD approach are shown in Fig. 2.12. When z true 1 is located within the solvable regions above the sensing depth s , the estimated depths z est 1 follow closely the 1:1 line and have smaller standard deviation (or uncertainty) among the simulations. However, when z true 1 is in the ambiguous regions (appeared as spikes in Fig. 2.12), the LPD solutions signicantly overestimate z true 1 . Note that when smaller noise is present, as long as z true 1 has not reached the sensing depth, the solvable regions have narrower spread around the 1:1 line as the probability of synthesized radar data being contaminated into ambiguous regions gets lower. Once z true 1 gets closer to and beyond the sensing depth, z est 1 is distributed around a value close to s as the synthesized 0 values have lost the 29 sensitivity to z 1 , resulting in underestimation of z true 1 . The depths estimated by the global minimum of L(z 1 ) are also shown in Fig. 2.12(b) for the case of n = 0:5 dB, and they have larger biases and uncertainties due to the ambiguity of L(z 1 ) and the presence of noise. 2.6.3 Time-Series Retrieval One way to avoid overestimation when using the LPD approach for the 0 points in the ambiguous regions shown in Fig. 2.11(a) is to use a time-series of radar mea- surements taken when the depth of the dielectric substrate z s is barely moving but the top dielectric prole 0 p (z) is changing. In the context of ALT retrieval, this sce- nario can happen when 0 measurements are taken between the maximum thaw and before the upward freezing front rises signicantly from the permafrost table. For the two-layer problem discussed for the LPD approach, the changes of 0 1 will alter the equivalent wavelength in the rst soil layer, the locations of the extrema of the corresponding 0 (z 1 ) curve will also be shifted. This increases the chances of having z true 1 located outside the ambiguous regions for at least one time snapshot so that the LPD approach can be applied to nd a good estimate of z true 1 . Furthermore, most of the local minima of L(z 1 ) will occur at dierent depths for dierent time snapshots, whereasz true 1 will remain at the same location. Therefore, a joint cost function using all time-series observations will have fewer local minima at which the cost function is lower than the error criterion. The cost function involving N t time-series observations can be dened as L(X(t 1 ;t 2 ;:::;t Nt )) = 1 N X pq;f;tn j 0 pq (X(t n );f; inc (t n ))d pq (f; inc (t n ))j 2 ; (2.8) where X(t n ) = ( 0 i (t n );z i (t n )) are the unknown variables when t n -th time snapshot measurement 0 (t n ) is taken. Fig. 2.13 shows the cost function using a two-element 30 0 0.2 0.4 0.6 0.8 z 1 [m] 0 0.25 0.5 0.75 1 L(X) L(X(t 1 )) LPD (t 1 only) n 2 (a) 0 0.2 0.4 0.6 0.8 z 1 [m] 0 0.25 0.5 0.75 1 L(X) L(X(t 2 )) LPD (t 2 only) n 2 (b) 0 0.2 0.4 0.6 0.8 z 1 [m] 0 0.25 0.5 0.75 1 L(X) L(X(t 1 ; t 2 )) LPD (t 1 & t 2 ) min L(X(t 1 ; t 2 )) (c) Figure 2.13: Example demonstrating cost function values of a two-layer problem with one unknown z 1 , using two snapshots (L(X(t 1 )) and L(X(t 2 ))) and the time-series combined (L(X(t 1 ;t 2 ))). ( 0 1 (t 1 ); 0 2 ) = (30; 5) at t 1 and ( 0 1 (t 2 ); 0 2 ) = (15; 5) at t 2 . The true value of unknown z 1 is set to z true 1 = 0:4 m. The vertical pink dotted lines correspond to the solutions estimated by each time snapshot LPD approach. The vertical black dashed-dot lines are the search space for the LPD-assisted time-series approach. 31 time-series observation ( 0 (t 1 ) and 0 (t 2 )) along with the ones using just the snap- shot at t 1 or at t 2 . The values of 0 1 (t 1 ), 0 1 (t 2 ), and 0 2 are again assumed known to the problem, which leaves z 1 as the only unknown. Due to the shift of the locations of the local minima in the cost functions (L(X(t 1 ) and L(X(t 2 )), the possible set of solutions is signicantly smaller when using the joint cost function L(X(t 1 ;t 2 )). The global minimum of L(X(t 1 ;t 2 )) is still not able to completely resolve the ambiguity between those local minima at which the cost function is lower than 2 n . To regular- ize L(X(t 1 ;t 2 )), we use the LPD solution of each snapshot (z LPD 1 (t 1 ) and z LPD 1 (t 2 )) as prior knowledge of z true 1 , and the z 1 estimate using this time-series approach is determined as z TS = arg min minfz LPD 1 (t 1;2 )grz 1 maxfz LPD 1 (t 1;2 )g+r L(z 1 (t 1 ;t 2 )); (2.9) wherer is a parameter which determines the search space inz 1 and is set to maxf eff;1 (t 1 )=2; eff;1 (t 2 )=2g to account for the resolution of each snapshot LPD solutions. The ex- ample cost function in Fig. 2.13(c) shows how the LPD-assisted time-series approach helps nding a better estimate of z true 1 than the one simply determined by the global minimum of L(X(t 1 ;t 2 )). To show that the time-series approach can remedy the overestimation of the LPD solutions, synthesized data representing snapshots of radar backscattering measure- ments are added one by one in the Monte-Carlo simulations. The value of 0 2 of a two-layer dielectric structure is assumed to be xed at 5, while 0 1 (t n ) is set to 30, 15, and 25, at t 1 , t 2 , and t 3 , respectively. Fig. 2.14 shows that when including more snapshots in the time-series, the overestimation is mitigated as their ambiguous re- gions are mismatched and the LPD solutions can better estimate z true 1 , compared to the LPD solutions using just the snapshot at t 1 in Fig. 2.12(a). The solutions determined by the global minimum of L(X(t 1 ;t 2 )) have also been greatly improved 32 0 0.2 0.4 0.6 0.8 z 1 true [m] 0 0.2 0.4 0.6 0.8 1 z 1 est [m] LPD (t 1 & t 2 ) s (a) 0 0.2 0.4 0.6 0.8 z 1 true [m] 0 0.2 0.4 0.6 0.8 1 z 1 est [m] min L(X(t 1 ; t 2 )) s (b) 0 0.2 0.4 0.6 0.8 z 1 true [m] 0 0.2 0.4 0.6 0.8 1 z 1 est [m] LPD (t 1 & t 2 & t 3 ) s (c) 0 0.2 0.4 0.6 0.8 z 1 true [m] 0 0.2 0.4 0.6 0.8 1 z 1 est [m] min L(X(t 1 ; t 2 ; t 3 )) s (d) Figure 2.14: Comparison of error performance of z est 1 . (a) LPD-assisted time-series approach using t 1 and t 2 . (b) Global minimum approach using t 1 and t 2 . (c) LPD- assisted time-series approach usingt 1 ,t 2 andt 3 . (d) Global minimum approach using t 1 , t 2 , and t 3 . ( 0 1 (t 1 ); 0 2 ) = (30; 5); ( 0 1 (t 2 ); 0 2 ) = (15; 5); ( 0 1 (t 3 ); 0 2 ) = (25; 5). The gray shaded regions represent one standard deviation around the average values of z est 1 . n is set to be 0.5 dB. The vertical dashed lines are the expected sensing depths s for the snapshots t 1 , t 2 , and t 3 . by having time-series observations, and as the number of snapshots increases, the global minimum approach may perform better than the LPD-assisted approach when the regularization is no longer needed. The root-mean-squared error (RMSE) of the z 1 estimates by the LPD-assisted time-series approach is shown in Fig. 2.15. When more snapshot measurements are included in the retrieval and when z true 1 is within the sensing depth s , the RMSE values can be progressively reduced below 10 cm, 33 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 z 1 true [m] 0 0.1 0.2 0.3 0.4 0.5 RMSE [m] LPD (t 1 only) LPD (t 1 & t 2 ) LPD (t 1 & t 2 & t 3 ) s Figure 2.15: Root-mean-squared error (RMSE) of the LPD-assisted time-series ap- proach. as the overestimation of z true 1 by the LPD solutions is mitigated by the time-series observations. However, when z true 1 is beyond the sensing depth s , the RMSE values increase nearly linearly with increasing z true 1 for all cases because the corresponding values of 0 (t n ) have lost the sensitivity to z 1 and the underestimation of z true 1 is inevitable due to the limited sensing depths of radar signals. 2.7 Retrieval of Active Layer Properties The retrieval of active layer properties in this work is performed using a two- element time-series of 0 measurements, one taken in mid-late August (t 1 ) and one taken in early October (t 2 ). Based on the observations made from the in-situ data in Fig. 2.1 and that the measured soil moisture can be seen as a proxy for soil dielectric constant, we assume a three-layer dielectric structure with a top surface roughness of h to represent the surface active layer ( 0 1 ), the bottom saturated active layer ( 0 2 ), and the permafrost ( 0 f = 5), as shown in Fig. 2.2. The top layer dielectric constant 0 1 can vary within the time window to account for the changes in soil moisture or freeze-thaw state in the near-surface layer. The surface may be covered by snow in early October; however, since the dielectric constant of wet snow is also low (3.25 as indicated in [36]) 34 and close to that of frozen soils, the snow cover is lumped into the surface layer and the value of 0 1 represents the degree of frozen state at the surface. Because the water content of saturated soils only depends on soil porosity, 0 2 is considered as time- invariant in the problem. The depths z 1 and z 2 represent the surface layer thickness and ALT, respectively. Althoughz 2 can vary between the two AirMOSS acquisitions due to the upward freezing, we assume that the amount ofz 2 variations is still within the uncertainty of the time-series approach shown in Fig. 2.14(a) and treat z 2 as another unknown that is constant between t 1 and t 2 . The assumption about z 2 may fail if the upward freezing front has risen signicantly att 2 , which could introduce an underestimation of the retrieved ALT. Finally, surface roughness h is also assumed time-invariant and the ratio of correlation length to rms height (l c =h) is set to 10 in this work as the retrieval performance is not expected to be sensitive to errors in knowledge ofl c [26]. There are thus four input 0 measurements ( 0 hh (t 1 ); 0 vv (t 1 ); 0 hh (t 2 ); 0 vv (t 2 )) and seven unknowns consisting of three dielectric constant values ( 0 1 (t 1 ), 0 1 (t 2 ), 0 2 ), three depths (z 1 (t 1 ), z 1 (t 2 ), z 2 ), and surface roughness (h). Before the LPD-assisted time-series approach can be applied to estimate z 2 , we need to estimate the other six unknowns rst. Since the inverse problem we have here is ill-posed, we do not attempt to nd the global minimum of the cost function dened in (2.8); instead, we search for the solution points that have cost function values below an error criterion (L(X(t 1 ;t 2 ))< 10 3 is used in this work). This search is performed by using the simulated annealing (SA) method [34], which is a global optimizer using a randomized search method. The stopping criterion is set to 10 3 so that the SA algorithm can exit a local minimum with cost function value smaller than that when the random search has converged to this point for a certain number of times. For each randomized search output found by the SA algorithm, the vector of unknowns X SA = [h; 0 1 (t 1 ); 0 1 (t 2 ); 0 2r ;z 1 (t 1 );z 1 (t 2 );z 2 ] is estimated. The estimates of all unknowns except z 2 will be considered as known parameters in the LPD-assisted 35 Table 2.1: AirMOSS Flights Used in this Study Flight Line Start Coordinate End Coordinate Acquisition Dates Atqasuk (69.6838°, -156.8095°) (71.0499°, -157.9022°) 2014/08/16, 2014/10/09, 2015/08/29, 2015/10/01 Barrow (71.3803°, -156.7062°) (70.2327°, -154.6705°) Deadhorse (70.4471°, -148.4280°) (68.7929°, -148.8882°) Ivotuk (67.9607°, -156.5759°) (68.8992°, -155.0969°) time-series method, from which the nal solution to z 2 is estimated. The two steps (SA algorithm followed by the LPD-assisted time-series approach) are considered as a retrieval run, and we perform at least 30 runs to obtain the statistics of each unknown variable. For each unknown variable, the mean and the standard deviation of the solutions from all retrieval runs are reported as the retrieved value and the retrieval uncertainty, respectively. 2.7.1 Retrieval Results from Pre-ABoVE IDS Flights The dates and the ight lines of the AirMOSS P-band radar data used in this work to estimate the ALT in 2014 and 2015 are summarized in Table 2.1. The radar pixels are multilook averaged from original 0.5 arcsec resolution to 2 arcsec in order to reduce the speckle eects. The study areas are focused on the ight lines within the Alaska North Slope (the four northern lines in Fig. 1.2). They are located within the continuous permafrost zone and the aboveground vegetation consisting of dwarf shrub and tussock/sedge/moss tundra has minimal impact on P- band radar backscatter and can be ignored, as indicated by the simulations using the scattering model in [37]. To determine the pixels on which the retrieval algorithm will be performed, the land cover information from the National Land Cover Database (NLCD) [38] for the state of Alaska from the year 2011 is used. The pixels classied as open water or developed areas are excluded. We also exclude the forest pixels as the forward scattering model used in the current work only applies to the bare or lightly vegetated surfaces. The forward model can be extended to include taller and denser vegetation, such as forested areas, for retrievals of ALT in lower latitudes, especially 36 south of the Brooks range (the six southern lines in Fig. 1.2). We have already developed this forward model, and have validated its use for the RZSM retrievals in several sites in north America [22, 37, 39]. For application to ALT retrievals in Alaska, the model needs to be parameterized for the corresponding vegetated land covers within the ight lines. This is the subject of our ongoing work but is not within the scope of the present paper. Fig. 2.16 shows the maps of active layer properties (h, 0 1 (t 1 ), 0 1 (t 2 ), 0 2 , z 1 (t 1 ), z 1 (t 2 ),z 2 ) retrieved from the 2015 AirMOSS data over the Deadhorse ight line which covers the northern section of Dalton Highway. The retrieved 0 2 is shown to be much higher than both 0 1 (t 1 ) and 0 1 (t 2 ) as it represents the unfrozen water content of the saturated zone in the active layer. The values of 0 1 (t 2 ) are generally lower than 0 1 (t 1 ) for most of the pixels, as the surface had started to freeze in early October and the unfrozen water content had decreased further. The values of z 1 (t 1 ) and z 1 (t 2 ) are considered as the unfrozen water table depths in August and October, respectively, and the results show that the unfrozen water table had moved deeper in October, as the surface layer became partially frozen and had less unfrozen water content. The retrieval uncertainty of each variable is shown in Fig. 2.17 and is generally larger for deeper subsurface features such as 0 2 andz 2 . Note that the two-layer dielectric model assumed for the active layer is an approximation of the actual soil dielectric prole, which may have more complicated prole dynamics caused by soil moisture gradients and SOM concentration. But also note that the two-layer model does not preclude the presence of an organic layer; the dielectric constant of each layer can represent that of mineral soil, organic soil, or their mixture. The radar retrieval of dielectric constant is agnostic to the type of soil. However, if the values of dielectric constant need to be converted to soil moisture, then the specic soil mixture must be known so that an appropriate dielectric mixture model can be applied for conversion to soil moisture. 37 -149 -148.6 -148.2 Longitude 68.8 69 69.2 69.4 69.6 69.8 70 70.2 70.4 Latitude 0 10 20 30 40 50 60 (a) 0 1 (t 1 ) -149 -148.6 -148.2 Longitude 68.8 69 69.2 69.4 69.6 69.8 70 70.2 70.4 Latitude 0 10 20 30 40 50 60 (b) 0 1 (t 2 ) -149 -148.6 -148.2 Longitude 68.8 69 69.2 69.4 69.6 69.8 70 70.2 70.4 Latitude 0 10 20 30 40 50 60 (c) 0 2 -149 -148.6 -148.2 Longitude 68.8 69 69.2 69.4 69.6 69.8 70 70.2 70.4 Latitude 0 0.05 0.1 0.15 0.2 0.25 0.3 (d) z 1 (t 1 ) [m] -149 -148.6 -148.2 Longitude 68.8 69 69.2 69.4 69.6 69.8 70 70.2 70.4 Latitude 0 0.05 0.1 0.15 0.2 0.25 0.3 (e) z 1 (t 2 ) [m] -149 -148.6 -148.2 Longitude 68.8 69 69.2 69.4 69.6 69.8 70 70.2 70.4 Latitude 0.3 0.35 0.4 0.45 0.5 0.55 0.6 (f) z 2 (ALT) [m] -149 -148.6 -148.2 Longitude 68.8 69 69.2 69.4 69.6 69.8 70 70.2 70.4 Latitude 0 0.01 0.02 0.03 0.04 (g) h [m] -149 -148.6 -148.2 Longitude 68.8 69 69.2 69.4 69.6 69.8 70 70.2 70.4 Latitude (h) NLCD Land Cover (i) NLCD Land Legend Figure 2.16: Retrieved maps of layer dielectric constants i , layer vertical positionsz i , ALT, and surface roughness h using the AirMOSS data acquired for the Deadhorse ight line over the northern section of Dalton Highway in 2015. Similar maps have been generated for 2014 observations but not shown for the purpose of brevity. 38 -149 -148.6 -148.2 Longitude 68.8 69 69.2 69.4 69.6 69.8 70 70.2 70.4 Latitude 0 2 4 6 8 10 12 14 16 18 20 (a) 0 1 (t 1 ) -149 -148.6 -148.2 Longitude 68.8 69 69.2 69.4 69.6 69.8 70 70.2 70.4 Latitude 0 2 4 6 8 10 12 14 16 18 20 (b) 0 1 (t 2 ) -149 -148.6 -148.2 Longitude 68.8 69 69.2 69.4 69.6 69.8 70 70.2 70.4 Latitude 0 2 4 6 8 10 12 14 16 18 20 (c) 0 2 -149 -148.6 -148.2 Longitude 68.8 69 69.2 69.4 69.6 69.8 70 70.2 70.4 Latitude 0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0.18 0.2 (d) z 1 (t 1 ) [m] -149 -148.6 -148.2 Longitude 68.8 69 69.2 69.4 69.6 69.8 70 70.2 70.4 Latitude 0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0.18 0.2 (e) z 1 (t 2 ) [m] -149 -148.6 -148.2 Longitude 68.8 69 69.2 69.4 69.6 69.8 70 70.2 70.4 Latitude 0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0.18 0.2 (f) z 2 (ALT) [m] -149 -148.6 -148.2 Longitude 68.8 69 69.2 69.4 69.6 69.8 70 70.2 70.4 Latitude 0 0.001 0.002 0.003 0.004 0.005 0.006 0.007 0.008 0.009 0.01 (g) h [m] Figure 2.17: Uncertainty maps of retrieved layer dielectric constants i , layer vertical positionsz i , ALT, and surface roughnessh using the AirMOSS data acquired for the Deadhorse ight line over the northern section of Dalton Highway in 2015. The major land cover in the southern and eastern parts of the Deadhorse ight line in Fig. 2.16 is classied as dwarf scrubs, which are usually dominated by tussock tundra and have thicker moss and organic layers, compared to the areas classied as sedges in the north. The surface moss and organic layers have higher porosity and can act as an insulating layer to prevent the active layer from thawing into deeper depths. As a result, the distribution of ALT values along this ight line does not follow the North-South temperature gradient (deep ALT in the south and shallow ALT in the north) and is mainly controlled by land cover types, as shown in the retrieved ALT 39 0.3 0.4 0.5 0.6 ALT [m] Barren Land Herbaceous Wetlands Sedge Dwarf Scrub Shrub (a) ALT - 2014 0.3 0.4 0.5 0.6 ALT [m] Barren Land Herbaceous Wetlands Sedge Dwarf Scrub Shrub (b) ALT - 2015 20 30 40 50 60 2 Barren Land Herbaceous Wetlands Sedge Dwarf Scrub Shrub (c) 0 2 - 2014 20 30 40 50 60 2 Barren Land Herbaceous Wetlands Sedge Dwarf Scrub Shrub (d) 0 2 - 2015 0 0.01 0.02 0.03 0.04 h [m] Barren Land Herbaceous Wetlands Sedge Dwarf Scrub Shrub (e) h - 2014 0 0.01 0.02 0.03 0.04 h [m] Barren Land Herbaceous Wetlands Sedge Dwarf Scrub Shrub (f) h - 2015 Figure 2.18: Distribution of the 2014 and 2015 AirMOSS retrievals (ALT, 0 2 , and h) for the Deadhorse ight line across dierent NLCD land cover types: 31-Barren land, 95-Emergent Herbaceous Wetlands, 72-Sedge/Herbaceous, 51-Dwarf Scrub, 52- Shrub/Scrub. The red central mark in each box indicates the median, and the bottom and top edges of the box indicate the 25th and 75th percentiles, respectively. The whiskers are extended to the 5th and 95th percentiles. in Fig. 2.16(f) and in the box plots in Fig. 2.18(a)(b). The spatial distribution of the retrieved 0 2 can also reveal the SOM content in the active layer. For the areas classied as dwarf scrubs or shrubs, the retrieved 0 2 is mostly higher than 40 as shown in Fig. 2.18(c)(d), which is only possible for organic soils as shown in Fig. 2.16(c), indicating the corresponding areas can potentially have higher SOM content. Moreover, the landscapes covered by tussock tundra are expected to have rougher ground surfaces, which is again consistent with the retrieved surface roughnessh shown in Fig. 2.16(g) and the box plots in Fig. 2.18(e)(f). Fig. 2.19 shows the maps of the retrieved 0 2 and ALT for the other three ight 40 -158 -157.5 -157 Longitude 69.8 70 70.2 70.4 70.6 70.8 71 Latitude 0 10 20 30 40 50 60 (a) Atqasuk 0 2 -156.5 -156 -155.5 -155 -154.5 Longitude 70.2 70.4 70.6 70.8 71 71.2 71.4 Latitude 0 10 20 30 40 50 60 (b) Barrow 0 2 -156.5 -156 -155.5 -155 Longitude 68 68.1 68.2 68.3 68.4 68.5 68.6 68.7 68.8 68.9 Latitude 0 10 20 30 40 50 60 (c) Ivotuk 0 2 -158 -157.5 -157 Longitude 69.8 70 70.2 70.4 70.6 70.8 71 Latitude 0.3 0.35 0.4 0.45 0.5 0.55 0.6 (d) Atqasuk z 2 (ALT) [m] -156.5 -156 -155.5 -155 -154.5 Longitude 70.2 70.4 70.6 70.8 71 71.2 71.4 Latitude 0.3 0.35 0.4 0.45 0.5 0.55 0.6 (e) Barrow z 2 (ALT) [m] -156.5 -156 -155.5 -155 Longitude 68 68.1 68.2 68.3 68.4 68.5 68.6 68.7 68.8 68.9 Latitude 0.3 0.35 0.4 0.45 0.5 0.55 0.6 (f) Ivotuk z 2 (ALT) [m] Figure 2.19: Retrieved maps of bottom layer dielectric constants 0 2 and ALT for the ight lines own in 2015 over Atqasuk, Barrow, and Ivotuk. The East-West gradient in (a) is an artifact related to the wide range of incidence angles within the AirMOSS imagery, which is expected to improve in the future when enhanced calibration of antenna patterns becomes available. 41 0.2 0.3 0.4 0.5 0.6 0.7 0.8 CALM ALT [m] 0.2 0.3 0.4 0.5 0.6 0.7 0.8 Retrieved ALT [m] West Dock (+0.07 m, 0.3) Deadhorse (-0.18 m, 6.2) Franklin Bluffs (-0.24 m, 4.4) Happy Valley 1 ha (+0.03 m, 0.1) Happy Valley 1 km (+0.01 m, 0.02) 56 Mile (-0.15 m, 1.5) Sagwon Hills MNT (-0.10 m, 0.7) Sagwon Hills MAT (+0.05 m, 0.4) Barrow (+0.10 m, 1.4) Atqasuk (-0.06 m, 0.2) Ivotuk (-0.05 m, 0.1) Site Name (bias, 2 ) (a) 2014 0.2 0.3 0.4 0.5 0.6 0.7 0.8 CALM ALT [m] 0.2 0.3 0.4 0.5 0.6 0.7 0.8 Retrieved ALT [m] West Dock (-0.01 m, 0.01) Deadhorse (-0.23 m, 9.2) Franklin Bluffs (-0.20 m, 3.5) Happy Valley 1 ha (+0.02 m, 0.1) Happy Valley 1 km (-0.004 m, 0.002) 56 Mile (-0.16 m, 1.5) Sagwon Hills MNT (-0.10 m, 0.6) Sagwon Hills MAT (+0.04 m, 0.3) Barrow (+0.04 m, 0.2) Atqasuk (-0.07 m, 0.1) Site Name (bias, 2 ) (b) 2015 Figure 2.20: AirMOSS ALT retrievals and CALM in-situ ALT comparison. The red solid (blue open) symbols represent the sites where the in-situ ALT is larger (smaller) than the radar sensing depth s . The error metrics are shown in the parenthesis after the site name in the legend for each site. The rst number in the parenthesis is the retrieval bias; the second number is the 2 statistic introduced in [6]. lines own in 2015 over Atqasuk, Barrow, and Ivotuk. The retrieved ALT (z 2 ) is validated against the probing measurements taken at the CALM sites within the four radar swaths, as shown in Fig. 2.20. The locations, sampling grids, and the dates when ALT measurements were taken for the CALM sites are summarized in Table 2.2. The ALT validation is performed at the site scale, where the retrieved ALT values (from all retrieval runs) of all AirMOSS pixels falling within a CALM grid are taken into account when calculating the average value and standard deviation for the site. The standard deviation of the retrieved ALT at the site scale represents both the retrieval uncertainty and the spatial ALT variability. For example, there are 15 AirMOSS pixels within the 1-hectare CALM grid at Happy Valley in the retrieved maps of Fig. 2.16. The mean values of these pixels range from 0.38 m to 0.5 m and the standard deviations are from 0.1 m to 0.19 m, which represent the retrieval uncertainty of each pixel. The standard deviation of the mean values represents the 42 Table 2.2: CALM in-situ ALT Sites Covered by the AirMOSS Flights in this Study Site Code Site Name Latitude Longitude Measurement Date Sampling Grid U1 Barrow 71.3103° -156.5928° 2014/08/18, 2015/08/21 1-km 2 U3 Atqasuk 70.4514° -157.4008° 2014/08/21, 2015/08/25 1-km 2 U4 West Dock 70.3745° -148.5522° 2014/08/25, 2015/08/25 1-ha U5 West Dock 70.3705° -148.5688° 2014/08/09, 2015/08/11 1-km 2 U6 Deadhorse 70.1613° -148.4653° 2014/08/25, 2015/08/26 1-ha U8 Franklin Blus 69.6739° -148.7219° 2014/08/26, 2015/08/27 1-ha U9A Happy Valley 69.1565° -148.8374° 2014/08/26, 2015/08/25 1-ha U9B Happy Valley 69.1495° -148.8534° 2014/08/13, 2015/08/15 1-km 2 U26 Ivotuk 68.4809° -155.7366° 2014/08/19, N/A 1-km 2 U31 56 Mile 69.6960° -148.6849° 2014/08/13, 2015/08/15 1-ha U32A Sagwon Hills MNT 69.4400° -148.6728° 2014/08/14, 2015/08/16 1-ha U32B Sagwon Hills MAT 69.4001° -148.8084° 2014/08/14, 2015/08/16 1-ha spatial variability at the 2-arcsec scale and is about 0.04 m for the site. The overall uncertainty of all retrieved ALT values within the grid is about 0.14 m. For the sites at West Dock, where both 1-hectare and 1-km 2 grids are sampled and they are adjacent to each other, we use an area-weighted average of their measurements to represent the overall in-situ ALT of the area. We also compare the radar sensing depths s , which are calculated during the retrieval process for each pixel and depend on the retrieved 0 i and z 1 , with the in-situ ALT values. The retrieved ALT values follow the 1:1 line well when the in-situ ALT is within the sensing depth (open symbols in Fig. 2.20) and the retrieval errors are generally less than 0.1 m with an overall bias of +0.05 m and +0.01 m for 2014 and 2015, respectively. When the in-situ ALT is beyond the sensing depth (solid symbols in Fig. 2.20), the retrieved ALT values underestimate the in-situ ALT with an overall bias of -0.13 m and -0.15 m for 2014 and 2015, respectively. Despite the underestimation due to limited sensing depths, the retrieved ALT values still show an increasing trend when the in-situ ALT gets deeper. We also adopt the 2 statistic introduced in [6] to evaluate the retrieved ALT values, which is dened as 2 = ( ALT SAR ALT 0 SD(ALT 0 ) ) 2 ; (2.10) whereALT SAR andALT 0 are the retrieved and in-situ ALT, respectively, andSD(ALT 0 ) 43 is the standard deviation in the in-situ ALT. The 2 value is used to determine whether the agreement between the retrieved and in-situ ALT is within the uncer- tainty of in-situ data (SD(ALT 0 )), with the values less than one indicating a good agreement. For all but three sites, the 2 values are less than one. For the sites at Deadhorse, Franklin Blus, and 56 Mile, the ALT retrievals have 2 values larger than one for both 2014 and 2015 as these sites have the deepest in-situ ALT (0.64 - 0.73 m) among the CALM sites studied. The radar sensing depth poses the physical limitation of using polarimetric SAR of a given system noise oor and a given calibration accuracy to quantify ALT | insucient sensitivity of radar backscatter when thaw depths are beyond sensing depths. Although the proposed method may not be able to fully capture the large gradients in ALT over the vast areas in high-latitude regions, it can still provide valuable information at the regional scale from the spatial patterns of the retrieved active layer properties, which includes not only ALT, but also the dielectric properties and surface roughness. It can also provide a measure of the lower-bound of the absolute value of ALT. Furthermore, the retrieval results do show that the retrieved ALT, 2 , and surface roughness are highly correlated with the land cover types, which are expected to be a key factor controlling the active layer dynamics. The values of the retrieved 1 can also reveal the changes in unfrozen water content during the soil freezing process. 2.7.2 Retrieval Results from 2017 ABoVE Flights The ten AirMOSS ight lines in northern Alaska that were own ve times in 2014 and 2015 during the span of the pre-ABoVE IDS project are referred to as "legacy lines" (Fig. 1.2). As part of the 2017 ABoVE Airborne Campaign (AAC), these legacy lines were re- own three times (early June, mid August, and early October), whereas most other foundational ight lines in Fig. 1.3 were only own two times (May/June 44 Figure 2.21: Retrieved ALT for the legacy ight lines and the ones over Toolik Lake and Teller, Alaska own during the ABoVE Airborne Campaign in 2017. and August/September). Since the time-series for ALT retrievals in this work require at least two radar snapshots after the maximum thaw in mid- or late-August, we only applied the LPD-assisted time-series approach to the AirMOSS ight lines of 2017 AAC that were own both in August and October. The retrieved ALT of these ight lines are shown in the map of Fig. 2.21, which include the ten legacy lines, one line over Toolik Lake (south of the Deadhorse line and north of the Brooks Range), and one line over Teller, Alaska in Seward Peninsula (crossing the Council line). Due to warmer average air temperatures, the retrieved ALT values of the legacy lines south of the Brooks Range and the Teller line are generally larger than the ones in the north, with the exception that the northern part of the Deadhorse line also shows large ALT. As discussed in the previous section for Fig. 2.18, the soil and land cover 45 Figure 2.22: Interannual comparisons of August surface soil moisture and ALT maps retrieved from AirMOSS data acquired over the Deadhorse legacy line in 2014, 2015, and 2017. types are also important factors that can impact the soil thawing processes, which is why the retrieved ALT maps do not necessarily follow the North-South temperature gradient and can show the heterogeneity of ALT distributions on a regional scale. Fig. 2.22 shows the interannual changes in August surface soil moisture and ALT for the Deadhorse legacy line. The August surface soil moisture is converted from the retrieved surface soil dielectric constant 0 1 (t 1 ) using the calibration model developed by Engstrom et al. for near surface soil moisture in the Barrow region [40]. The results show that the surface soil moisture near the Prudhoe Bay area (northern 46 part of the retrieved swath that is close to the shoreline) is much wetter in 2014 compared to the ones in 2015 and 2017, which could be due to heavier precipitation and ooding or more ground ice was thawed in that summer. When surface soil is wet and close to saturation, the heat from the warmer atmosphere is easier to be transfered into subsurface soil because the thermal conductivity of liquid water (0.57 W/m/K) is higher than the one of air (0.023 W/m/K). The heat transfer is even more faster when the soil contains more mineral soil solids as they are better thermal conductors than organic matter [41]. Because the SOM content in Prudhoe Bay area is potentially lower than the southern part of the swath, the high soil moisture can greatly promote larger thaw depths in summer, which could be the reason why the 2014 retrieved ALT around the area is the largest among the three. Similar to the retrieval results of 2014 and 2015 data, the 2017 ALT retrievals are also underestimated when the actual permafrost table is beyond the sensing depth, as shown in Fig. 2.23(a). The underestimation of deep ALT is inevitable for polarimetric SAR-only approach, at least until the radar system calibration accuracy can be fur- ther improved or lower and multiple frequency bands can be used in the future. Fig. 2.23(b) shows the overall retrieval performance of ALT for all the AirMOSS acquisi- tions over the legacy ight lines so far (2014, 2015, 2017). When the ALT is within the sensing depth, the overall RMSE is 0.056 m with a bias of 0.029 m. Although for the sites when ALT is larger than the sensing depth, there is still a positive correlation between the retrieved and in-situ values, which means the retrieved ALT maps can still represent a good spatial distribution about the regional heterogeneity of ALT. All the results have been archived at the Oak Ridge National Laboratory Distributed Active Archive Center (ORNL DAAC): https://doi.org/10.3334/ORNLDAAC/1657. 47 (a) 2017 (b) All acquisition years (2014, 2015, 2017) Figure 2.23: AirMOSS ALT retrievals and CALM in-situ ALT comparison. (a) Val- idation plot of 2017 data. The purple (blue) symbols represent the sites where the in-situ ALT is larger (smaller) than the radar sensing depth s . The error metrics are shown in the parenthesis after the site name in the legend for each site. The rst number in the parenthesis is the retrieval bias; the second number is the 2 statistic introduced in [6]. (b) Validation plot of all acquisition years (2014, 2015, 2017). 48 2.8 Conclusion A permafrost active layer retrieval algorithm using time-series P-band SAR ob- servations is presented. The retrieval process is an inversion of the forward radar backscattering model, where a three-layer dielectric structure is used to model the active layer (surface and saturated layers) and the underlying permafrost. Both layer dielectric constants (representing unfrozen water content) and layer thicknesses (rep- resenting water tables and thaw depths) of the active layer soils are retrieved. The number of dielectric layers assumed in the model is limited by the available time-series measurement points and is aimed at characterizing the most prominent dielectric con- trasts in the active layer. The ability of the P-band radar for active layer retrievals is evaluated by the radar sensing depth, which is dened based on the sensitivity of the radar backscatter to the subsurface dielectric changes and the calibration accuracy of the radar system. Since the radar backscatter is an oscillating and decaying function of the layer interface position, the retrieval may become ambiguous such that multiple layer depths may correspond to similar radar backscattering coecient values. Time-series SAR obser- vations are used to reduce this ambiguity and increase the probability that the true solution for ALT falls within the region where the proposed largest possible depth approach can provide accurate ALT estimates. Mid-late August and early October ights are selected for the time-series to capture the maximum thaw and partially frozen seasons, when the boundary location between thawed active layer and per- mafrost is still observable and can be assumed constant over the time window. The dielectric constant of the lower part of the active layer is also assumed to be con- stant because the unfrozen water content depends only on soil porosity within the saturated zone. Having these time-invariant parameters can reduce the number of unknowns and improve the retrieval accuracy with limited number of observations. The ALT values estimated by the largest possible depth approach are used as prior 49 information of the true ALT values to regularize the search space in the optimization problem for the time-series approach. The combination of these two approaches sig- nicantly improves the ALT retrieval accuracy, as demonstrated by the Monte-Carlo simulations. The time-series retrieval algorithm is applied to the polarimetric SAR data ac- quired by AirMOSS over the ight lines on the Alaska North Slope from August and October of 2014 and 2015. The retrieved ALT is validated against the in-situ ALT measurements at CALM sites and error is quantied based on retrieval biases and 2 values. For the sites where the in-situ ALT is smaller than the radar sensing depth, the ALT retrieval shows a small overall bias of +0.05 m and +0.01 m for 2014 and 2015, respectively. Although the ALT retrievals are in general underestimated for the sites where deep ALT is observed, the retrieved ALT shows a consistent trend when compared to in-situ ALT. The 2 values of ALT retrievals are less than one for most of the study sites, indicating a good match between the retrieved and in-situ ALT. The spatial patterns of the retrieved active layer properties are highly correlated with the land cover types, indicating land cover as one of the main controlling factors at the regional scale. 50 CHAPTER III Harmonized Parametrization of Soil Properties for Multi-Physics Approach 3.1 Introduction Permafrost-aected soils are more complex than the ones in temperate or tropical zones because of seasonal freezing and thawing and the presence of perennially frozen ground underneath [42]. Soils are often saturated from the bottom of active layer because the underlying permafrost impedes subsurface drainage. The conditions of cold temperature and excessive moisture have led to low decomposition rates and the accumulation of organic matter [43, 44], which can strongly impact subsurface thermal dynamics and hydrologic processes in the active layer. The thermal and hydraulic properties of organic soils also vary greatly with the degree of decomposition, which can be quantied by the amount of bers that still remains in organic materials [45]. Fibers are fragments of plant tissue that are large enough to be retained on a sieve of selected size (e.g., a sieve with 0.1-mm open- ings), representing the undecomposed or partially decomposed portions in organic soil materials. Depending on the degree of decomposition, organic/peat soils can be classied as bric (least decomposed), hemic (moderately decomposed), and sapric (highly decomposed) soils. Because the ber content in organic soils is the fundamen- 51 tal parameter that determines bulk density, porosity, and pore size distribution, these organic soil classes can have signicantly dierent water-retention characteristics and hydraulic conductivity [46], which must be taken into account when modeling the permafrost soil dynamics. Radar remote sensing of surface or subsurface soil features, including soil moisture and freeze/thaw state, requires a soil dielectric mixing model that translates the soil-water mixture into dielectric constant (relative permittivity). The amount of unfrozen water content present in the soil is the dominant factor in determining the soil dielectric constant, particularly the ratio of bound water to free water and their dielectric behaviors [47]. Clay and sapric soils can better hold the water molecules around the surface of soil particles due to larger surface areas and micropore volumes, resulting in higher amount of bound water that can exist in a soil medium. To estimate active layer properties, including active layer thickness (ALT), we have developed a baseline algorithm using time-series P-band polarimetric synthetic aperture radar (SAR) observations as detailed in Chapter II [48]. The algorithm considers a three-layer dielectric structure (two in the active layer and one for per- mafrost) and retrieves the corresponding dielectric constants and layer thicknesses using time-series P-band data acquired when the active layer is fully thawed (late August) and partially frozen (early October). The major obstacle to directly retriev- ing soil moisture proles instead of layer dielectric constants was largely due to the lacks of high-quality organic matter prole data, a realistic soil prole parametriza- tion, and a soil dielectric model applicable to a wide range of organic matter content. 52 3.2 Soil Parametrization Scheme 3.2.1 Incorporating Fiber Content into Parametrization Fiber content (FC) is a proxy measure of degree of decomposition of organic soils, which is dened in this work as the dry mass fraction of bers in organic matter [45]: FC = M fiber M om (3.1) where M fiber and M om are the dry masses of bers and organic matter, respectively. Organic matter content (OM) is dened as OM = M om M soil (3.2) whereM soil is the dry mass of total soil. Both FC and OM are highly correlated with bulk density ( b = M soil /V soil , where V soil is the total volume of soil) as more brous or organic materials in the soil generally result in a more porous medium [45] [49]. We thus assume that there is a positive correlation between FC and OM and we use bulk density as the intermediate parameter to relate the two based on existing FC{ b and b {OM relationships reported in literatures. The FC{ b relation in [45] was a quadratic function, which is apparently not valid for high b values. We ret it with an exponential function to cover the full range of b . The tted exponential curve is given as FC = 0:9887e 5:512 b (3.3) The b {OM relation is also commonly expressed as an exponential function [49]: b = bm e rOM (3.4) 53 where bm is the bulk density of purely mineral soil (OM = 0) and r is the decay constant. bm can be calculated from the porosity of mineral soil ( m ) and the specic density of mineral particles ( m = 2.65 g=cm 3 ): bm = (1 m ) m (3.5) wherein we use the regression equation of m on percent sand (%sand) by [50]: m = 0:489 0:00126 %sand (3.6) The decay constant r can be determined by (3.4){(3.6) assuming a bulk density of 0.05 g=cm 3 for purely organic soil (OM = 1). We then combine the curves in (3.3) and (3.4) to express FC as a function of OM: FC = 0:9887 exp(5:512 bm exp(rOM)) (3.7) An example of the FC{ b , b {OM, and FC{OM relations with a sand fraction of 35% are shown in Fig. 3.1, where the data from the U.S. Department of Agriculture- National Cooperative Soil Survey (USDA-NCSS) Soil Characterization Database [51] are also shown for comparisons. Instead of dividing soil solids into just mineral and organic materials as most land surface models do, we incorporate FC into the parametrization scheme and consider three principal solid components: mineral (m), humus (h), and brous (f) materials. The humus and brous materials here represent the portions of soil organic matter that are well-decomposed and slightly decomposed, respectively. Their sub-phase 54 0 0.2 0.4 0.6 0.8 1 1.2 1.4 Bulk Density [g/cm 3 ] 0 20 40 60 80 100 Fiber Content [%] USDA-NCSS Data Boelter 1969 Fitted Exponential Curve Hemic Sapric Fibric (a) 0 20 40 60 80 100 Organic Matter [%] 0 0.5 1 1.5 2 2.5 Bulk Density [g/cm 3 ] USDA-NCSS Data Exponential Curve (b) 0 20 40 60 80 100 Organic Matter [%] 0 20 40 60 80 100 Fiber Content [%] USDA-NCSS Data Model (c) Figure 3.1: (a) FC vs. bulk density. (b) Bulk density vs. OM. (c) FC vs. OM. The sand fraction is 35% for this example. The curves are NOT derived from the USDA-NCSS data, which are shown here merely for comparisons. Note that the ber content in USDA-NCSS data is dened as a percentage by volume, not by mass. volumetric fractions (f vm , f vh , f vf , with f vm + f vh + f vf = 1) can be expressed as f vm = (1OM)= m (1)= b (3.8a) f vh = OM(1FC)= h (1)= b (3.8b) f vf = OMFC= f (1)= b (3.8c) where h = 1.8 g=cm 3 and f = 0.6 g=cm 3 are assumed to be the specic densities of 55 0 20 40 60 80 100 OM [%] 0 0.2 0.4 0.6 0.8 1 Volumetric Fraction, f v [-] Mineral (f vm ) Humus (f vs ) Fibrous (f vf ) (a) 0 20 40 60 80 100 OM [%] 0 0.2 0.4 0.6 0.8 1 [m 3 /m 3 ] Mixed soil Mineral soil Sapric peat Fibric peat (b) Figure 3.2: (a) Sub-phase volumetric fractions (f v ) of mineral, humus, and brous materials and (b) Soil porosity () as a function of organic matter content (OM). The sand fraction is 35% for this example. The porosity at 100% and 45% of OM correspond to the values of bric and sapric peat reported in Letts et al. (2000), respectively. humus and brous materials, respectively, and the porosity is given by = 1 b m (1OM) b h OM(1FC) b f OMFC (3.9) The sub-phase volumetric fractions and porosity are plotted as a function of OM in Fig. 3.2. Although we extended the description of soil by adding FC to characterize the degree of decomposition, the input variables to the parametrization are still min- eral soil texture (%sand and %clay) and OM, as FC being parametrized as a function of OM in (3.7). 3.2.2 Parametrization of Soil Hydraulic Properties The soil hydraulic properties considered in this work include the soil water reten- tion and hydraulic conductivity, for which the equations developed by Campbell [52] are used to relate the soil matric potential ( ) and hydraulic conductivity (k) to the 56 volumetric soil moisture (): = s ( s ) b (3.10a) k =k s ( s ) 2b+3 (3.10b) where s (m 3 =m 3 ) is the soil moisture content at saturation, which is equivalent to the porosity , s (cm) is the air-entry potential, b is an empirical parameter that describes the shape of the water retention curve, and k s (m=s) is the saturated hy- draulic conductivity. To estimate the hydraulic properties of soils with various min- eral texture and OM, we apply a weighted combination approach using the sub-phase fractions in (3.8a{3.8c) and the parameter values (X) that are representative for each principal solid component (X m ,X h ,X f ). The regression equations of hydraulic prop- erties, also known as pedotransfer functions (PTFs), in [50] provide the parameter values corresponding to purely mineral soil (X m ). The literature values reported in [46] are for bric, hemic, and sapric peat, which are mixtures of humus and brous materials. We assume that the bric and sapric peat in [46] correspond to the soils with OM values of 100% and 45%, respectively, and with a sand fraction of 35%. With this assumption, the corresponding porosity values calculated by (5.6) can match the ones reported in [46] (0.93 and 0.83 for bric and sapric peat, respectively), as shown in Fig. 3.2(b). We then use the sub-phase volumetric fractions (f vm;fibric , f vh;fibric , f vf;fibric , f vm;sapric , f vh;sapric , f vf;sapric ) and the hydraulic parameters (X fibric , X sapric ) of bric and sapric peat in [46] to solve for the nominal values of \purely" humus and brous materials (X h ,X f ). Finally, the hydraulic parameters are calculated from the weighted average of the nominal values of principal solid components. For the soil water retention curve in (3.10a), there are two parameters to be determined: the air-entry potential s and the exponent b of the curve. The values of s barely vary among the bric, hemic, sapric peat [46], showing a very weak 57 dependence on degree of decomposition. We thus assume a linear transition from bric peat to sapric peat, followed by a quadratic transition from sapric peat to purely mineral soil. s = 8 > > < > > : s;sapric + OM0:45 10:45 ( s;fibric s;sapric ); if OM > 0:45 s;sapric + sm s;sapric 0:45 2 (OM 0:45) 2 ; if OM 0:45 (3.11) where sm is the air-entry potential of purely mineral soil and is given by the PTF in [50]: sm =10 1:880:0131%sand cm (3.12) and s;fibric = -1.03 cm and s;sapric = -1.01 cm from [46]. For the exponent b, we use the weighted arithmetic mean to calculate b of mixed soils: b =f vm b m +f vh b h +f vf b f (3.13) where b m = 2:91 + 0:159 %clay as given by [50]. To obtain the values of b h and b f , we solve for the following system of equations: 8 > > < > > : f vm;fibric b m +f vh;fibric b h +f vf;fibric b f =b fibric f vm;sapric b m +f vh;sapric b h +f vf;sapric b f =b sapric (3.14) where b fibric = 2.7 and b sapric = 10 [46]. Note that we slightly adjusted the value of b sapric from 12 to 10 to avoid slow drainage in well-decomposed organic soils. The solutions to b f and b h in (3.14) are 0.654 and 21.1668, respectively. For saturated hydraulic conductivity k s in in (3.10b), we adopted the weighted geometric mean to calculate k s of mixed soils: k s =k fvm sm k f vh sh k f vf sf (3.15) 58 which is equivalent to the weighted arithmetic mean of log 10 k s that is used in [50]: log 10 k s =f vm log 10 k sm +f vh log 10 k sh +f vf log 10 k sf (3.16) wherek sm = 7:0556 10 6:884+0:0153%sand as given by [50]. The system of equations to solve for the values of k sh and k sf is 8 > > < > > : k f vm;fibric sm k f vh;fibric sh k f vf;fibric sf =k s;fibric k f vm;sapric sm k f vh;sapric sh k f vf;sapric sf =k s;sapric (3.17) where k s;fibric = 2:8 10 4 m=s and k s;sapric = 1:0 10 6 m=s [46]. Note that we again adjusted the value for sapric peat from 1:0 10 7 to 1:0 10 6 m=s to avoid too slow water ows in well-decomposed organic soils. The solutions tok sf andk sh in (3.17) are 9:412110 4 and 4:951710 9 m=s, respectively. The hydraulic parameters of soils with varying OM values are plotted in Fig. 3.3, showing that the proposed parametrization is able to capture the slow drainage or water ows in well-decomposed organic soils, which is caused by the better water-holding capability of humus materials. 59 0 20 40 60 80 100 OM [%] -30 -25 -20 -15 -10 -5 0 s [cm] Mixed soil Mineral soil Sapric peat Fibric peat (a) 0 20 40 60 80 100 OM [%] 0 5 10 15 b [-] Mixed soil Mineral soil Sapric peat Fibric peat (b) 0 20 40 60 80 100 OM [%] 10 -8 10 -7 10 -6 10 -5 10 -4 10 -3 k s [m/s] Mixed soil Mineral soil Sapric peat Fibric peat (c) Figure 3.3: Soil hydraulic properties as a function of organic matter content (OM). (a) Air-entry potential ( s ). (b) Exponent b of the Campbell water retention curve. (c) Saturated hydraulic conductivity (k s ). The mineral soil texture is to (sand;clay) = (35%; 10%) for this example. The values of bric and sapric peat are from the ones reported in Letts et al. (2000) with minor adjustments forb andk s to avoid too slow drainage in well-decomposed organic soils. 60 3.3 Soil Dielectric Mixing Model A generalized soil dielectric mixing model can be expressed as soil = (1) s + w + () a (3.18) where s , w , a are the dielectric constants of solid soils, water, and air, respectively, and is a constant shape factor [47]. Following the work by Park et al. [53], we in this work simply adopted a linear model ( = 1) to demonstrate the advantages of incorporating the parametrization in the previous sections into a soil dielectric mixing model. The dielectric constant of solid soils can be decomposed as s =f vm m + (f vh +f vf ) om (3.19) where m is the dielectric constant of mineral particles, for which we use the dielectric model of dry rocks in [54] m = 1:96 2:65 +i(0:038 + 0:06=f r ) (3.20) wheref r is the frequency in GHz. As for the dielectric constant of organic matter om , we use om = 1:7 from a dielectric model for a dry vegetation material in [55]. Both m and om have relatively small real parts, hence small contributions to the eective soil dielectric constant soil in (3.18). The dielectric constant of water, however, is the dominant factor in determining soil and can be expressed as w =f v;bw bw +f v;fw fw +f v;ice ice (3.21) where the subscripts bw and fw refer to bound water and free water, respectively. bw is set to 35 +i15 based on the work by Dobson et al. [47]. For fw and ice , a 61 double-Debye dielectric model (D3M) of saline water and a dielectric model of saline ice reported in [56] are used, respectively. To determine the sub-phase volumetric fractions of water in (3.21), we rst con- sider the case when the soil temperature (T ) is above 0 °C, i.e., only liquid water fractions are considered (f v;bw + f v;fw = 1). Following the approach by Park et al. [53], we use the wilting point ( wp ) to assign the bound-to-free water transition point. The wilting point is the minimal amount of water in the soil that the plant requires not to wilt, which is dened in convention as the soil moisture content at −15,000 cm of soil matric potential and can be expressed using the Campbell water retention curve parameters in (3.10a) as wp = s ( wp s ) 1=b (3.22) where wp is set to −15,000 cm. When the soil water content is smaller than the wilting point wp , the water content consists only of bound water, i.e., 8 > < > : f v;bw = 1 f v;fw = 0 if 0 wp (3.23) When the water content exceeds the wilting point, we assume that the bound water content is xed at wp and any additional water content is counted towards the free water, i.e., 8 > < > : f v;bw = wp f v;fw = 1 wp if wp s = (3.24) Fig. 3.4 shows the values of soil dielectric constant with respect to varying soil moisture and OM. The dielectric constant values of well-decomposed organic soils (20% { 60% of OM) are lower than mineral and slightly decomposed organic soils because of the higher amount of bound water content ( bw is much smaller than fw ) 62 0 20 40 60 80 100 OM [%] 0 0.2 0.4 0.6 0.8 1 Soil Moisture [m 3 /m 3 ] 0 10 20 30 40 50 60 70 80 Figure 3.4: Soil dielectric constant as a function of soil moisture and organic matter content (OM) at P-band (430 MHz). hold between micropore structures of humus materials. When the soil temperature drops below 0 °C, not all soil water will be frozen instantly; instead, the free water freezes before the bound water as the latter remains in liquid form due to the adsorptive and capillary forces exerted by soil particles, which is equivalent to a depression of the freezing point. The amount of soil water that is frozen at a certain temperature depends on the exchanges of thermal energy for releasing the latent heat to the surrounding soils, which can only be determined by a soil heat transfer model. To demonstrate the dielectric behaviors of soils at subzero temperatures, we assume that the soil water remains unfrozen for the portion that is smaller than the maximum liquid water content ( liq;max ) allowed at a subfreezing soil temperature (T k <T f ), which is given by liq;max = s ( 10 2 L f (T k Tf) gT s ) 1=b (3.25) where T k and T f are soil temperature and freezing point (K), respectively, L f is the 63 -10 -5 0 5 10 Soil Temperature [ ° C] 0 0.2 0.4 0.6 0.8 1 Volumetric Fraction, v [-] v bw v fw v ice v liquid v water Figure 3.5: Volumetric fractions of soil water phases in the total volume of soil (v = V x =V soil ). -30 -20 -10 0 10 20 30 Soil Temperature [° C] 0 20 40 60 80 100 OM [%] 0 5 10 15 20 25 30 35 40 45 Figure 3.6: Soil dielectric constant as a function of soil temperature and organic matter content (OM) at P-band (430 MHz). The saturation fraction (S w = = s ) is set to 70% in this example. 64 latent heat of fusion (J=kg), and g is the gravitational acceleration (m=s 2 ). Note that this equation also uses the Campbell water retention curve parameters ( s , s , b), for which we can again apply the proposed soil parametrization to obtain the liq;max of soils with various textures and organic matter contents. The volumetric fractions in the total soil volume in Fig. 3.5 shows that the free water freezes quickly when the temperature drops below 0 °C whereas the bound water contributes to the unfrozen water content when temperature is well below the freezing point. Note that the changes in volume from liquid water to ice are also considered in Fig. 3.5. As shown in Fig. 3.6, the dielectric constant values of well-decomposed organic soils can still maintain a value around 10 under subfreezing temperatures due to the existence of higher bound water content. 3.4 Model Validation To validate the soil parametrization scheme and dielectric mixing model proposed in this chapter, we use the eld data that we collected from four sites along the Dalton Highway, Alaska in August 2018. The visited sites include the following: • Sagwon (SW) (69°28' 40.55" N, 148°33' 53.08" W) • Happy Valley (HV) (69°9' 11.00" N, 148°50' 31.00" W) • Ice Cut (IC) (69°2' 33.89" N, 148°49' 34.59" W) • Imnavait Creek (IM) (68°36' 18.74" N, 149°18' 23.98" W) Detailed soil prole measurements were taken at a ne vertical resolution down to the permafrost table, which include dielectric constant, soil moisture, porosity, bulk density, organic matter content, and mineral texture. We use the measured soil moisture , OM, and mineral texture (%sand and %clay) as input parameters to the proposed dielectric mixing model in previous section and calculate the model 65 0 20 40 60 80 0 20 40 60 Depth [cm] HV-1 in-situ model (14.4) 0 20 40 60 80 0 20 40 60 Depth [cm] IC-1 in-situ model (5.3) 0 20 40 60 80 0 20 40 60 Depth [cm] IC-2 in-situ model (15.5) 0 20 40 60 80 0 20 40 60 Depth [cm] IM-1 in-situ model (9.7) 0 20 40 60 80 0 20 40 60 Depth [cm] SW-1 in-situ model (12.7) 0 20 40 60 80 0 20 40 60 Depth [cm] SW-2 in-situ model (9.4) Figure 3.7: Comparison of in-situ and model soil dielectric constants. The numbers in the parentheses in the legend are the Root Mean Square Error (RMSE) of each soil pit prole. dielectric constants for each soil layer samples. The comparisons between the in-situ and model dielectric constants are shown in Fig. 3.7, where the Root Mean Square Error (RMSE) of each soil pit prole is calculated. Without any further ne-tuning for the model parameters (e.g., in (3.18), bw , and the choice of bound-to-free water transition point) against in-situ or in-lab measurements, the model dielectric proles generally agree well with the in-situ dielectric proles. Since the soil dielectric and soil property measurements of a soil pit are not exactly co-located and the soil depths might be slightly biased when compressing the soil for sample compaction, some mismatches or larger errors are observed at certain soil depths. 66 3.5 Conclusion In this chapter, we presented a new soil parametrization scheme that considers both organic matter content and the degree of decomposition, which can capture the dierences in hydraulic properties not only between mineral and organic soils, but also between dierent types of organic materials (humus and bers). A soil di- electric mixing model is also presented to convert soil-water mixture into dielectric constant based on the parametrized water-retention characteristics. These model de- velopments and harmonized soil parameters show the potentials of directly retrieving soil moisture and organic matter content from radar observations, as well as leading towards a synergistic framework of radar remote sensing and permafrost soil mod- eling. Future work includes further model calibrations against extensive in-situ and in-lab measurements of various soil physical properties. In the following chapters, we will apply this soil parametrization scheme and dielectric model to the radar retrieval of soil moisture and organic matter proles in the Arctic-boreal regions. 67 CHAPTER IV Retrieval of Soil Moisture and Organic Matter Proles in the Arctic-Boreal Region 4.1 Introduction The capability of retrieving soil moisture using synthetic aperture radar (SAR) po- larimetry has been demonstrated through missions like Soil Moisture Active Passive (SMAP) and Airborne Microwave Observatory of Subcanopy and Subsurface (Air- MOSS). Radar backscatter is sensitive to soil dielectric properties, which are strongly correlated with soil moisture, texture, and freeze/thaw state. To enable accurate radar retrievals, we often rely on detailed soil maps for soil texture parameters in the forward model to convert dielectric constant to soil moisture. For retrieval of permafrost active layer properties, the spatial and vertical distribution of soil organic matter (SOM) is crucial because it can greatly aect subsurface soil hydraulic prop- erties, hence the dielectric properties of active layer soils. However, the SOM maps available for permafrost regions are still of limited quality and require more detailed soil surveys for better spatial coverage and resolution. Instead of treating organic matter (OM) prole as an ancillary data layer in the retrieval, in this chapter, we propose to retrieve organic matter and soil moisture proles simultaneously using P-band radar backscatter. To this end, eld observations 68 of permafrost-aected soils will be used to characterize and parametrize the OM proles in the active layer. The prole representation needs to be realistic and, at the same time, the number of parameters used in the prole functions must be kept as low as possible to avoid introducing too many unknowns in the inverse problem. The harmonized soil parametrization models described in Chapter III will be used to convert the organic matter and water content into soil dielectric constant. We will demonstrate a retrieval scheme to estimate active layer properties, including active layer thickness (ALT), water table depth, and soil moisture and OM proles. Finally, retrieval results using the AirMOSS data acquired as part of the Arctic-Boreal Vulnerability Experiment (ABoVE) airborne campaign in 2017 will also be presented. 4.2 Prole Representation of Permafrost Soils To observe and characterize the proles of various soil parameters in the per- mafrost active layer, eight soil pits were sampled at a 5- to 7-cm vertical increments at four sites north of the Brooks Range and along the Dalton Highway in Alaska. The site locations are summarized in Section 3.4. The key parameters for soil prole characterization are OM and soil moisture proles. Most soil structural (e.g., bulk density and porosity) and hydraulic properties can be derived from mineral texture and OM content using a pedotransfer function (PTF), for which we will exploit the soil parametrization scheme developed in Section 3.2. Once the porosity is determined from the OM content, we only need the saturation fraction to determine the volumet- ric soil moisture value. The soil pit data are analyzed and tted to a mathematical function so the prole shape can be properly captured while limiting the number of prole parameters. 69 4.2.1 Organic Matter Proles OM content is dened as the weight percentage of OM in a dry soil sample and can be measured by the loss-on-ignition (LOI) method, which combusts the organic carbon at a high temperature and uses the dierences in weights before and after the combustion to estimate how much carbon was in the soil sample. In many land surface models, OM proles is commonly assumed as an exponential function decreasing with soil depths [41] [49]: OM(z) =ae bz +c (4.1) The assumption of this exponential decay from the surface is based on the fact that most OM prole measurements are taken for soil horizons, which have large intervals along the soil prole. The exponential function in (4.1) is only valid when the soil proles considered in the model are sampled at a coarse vertical resolution. In Fig. 4.1, six in-situ OM proles measured from the soil pit samples are shown for compar- isons with tted prole functions. The exponential proles generally do not t well to the in-situ data when there is a surface organic layer thicker than the sampling resolution (5 to 7 cm) and the root-mean-squared error (RMSE) is larger than 10 % in such cases. When multiple points are sampled within the surface organic layer, the OM prole would look like a S-shaped curve and thus should be modeled as a sigmoid function. The generalized logistic function, also known as Richards' curve, is chosen to describe the OM proles in this work: OM(z) =A + KA 1 +e B(zM) (4.2) where A and K are the asymptotes as z approaches1 and1, respectively, B is the decay rate (B > 0), and M indicates the depth where the maximum decay rate occurs. We set K = 0:05 to represent the residual OM content in deep soils, 70 0 50 100 0 20 40 60 Depth [cm] HV-1 in situ exp (1.4) sigmoid (3.6) 0 50 100 0 20 40 60 IC-1 in situ exp (13.3) sigmoid (3.4) 0 50 100 0 20 40 60 IC-2 in situ exp (11.4) sigmoid (6.2) 0 50 100 OM [%] 0 20 40 60 Depth [cm] IM-1 in situ exp (16.8) sigmoid (8.9) 0 50 100 OM [%] 0 20 40 60 SW-1 in situ exp (0.6) sigmoid (1.6) 0 50 100 OM [%] 0 20 40 60 SW-2 in situ exp (18.4) sigmoid (3.9) Figure 4.1: Comparison of in-situ organic matter proles and two prole represen- tations: exponential and sigmoid functions. The numbers in the parentheses in the legend are the RMSE of each tted prole. which is valid for the near-surface soil proles where the cryoturbation eects are not prominent. After tting the generalized logistic function in (4.2) to the in-situ OM data, the tted values ofB parameter are around 50 for all proles, so we adopt this value forB hereafter. The remaining two parameters A andM will be the ones that control the shape of OM proles, with A representing the OM levels at the surface and M as an indicating depth for surface organic horizon. As shown in Fig. 4.1, the RMSE of the sigmoid function t is much smaller than the exponential t when thick surface organic horizon is present and is no larger than 9 % for all proles. With the OM proles properly parametrized, we can calculate other soil parame- ters using the soil parametrization scheme described in Section 3.4. As an example, 71 0 0.5 1 1.5 b [g/cm 3 ] 0 20 40 60 Depth [cm] HV-1 in-situ exp sigmoid 0 0.5 1 1.5 b [g/cm 3 ] 0 20 40 60 Depth [cm] IC-1 in-situ exp sigmoid 0 0.5 1 1.5 b [g/cm 3 ] 0 20 40 60 Depth [cm] IC-2 in-situ exp sigmoid 0 0.5 1 1.5 b [g/cm 3 ] 0 20 40 60 Depth [cm] IM-1 in-situ exp sigmoid 0 0.5 1 1.5 b [g/cm 3 ] 0 20 40 60 Depth [cm] SW-1 in-situ exp sigmoid 0 0.5 1 1.5 b [g/cm 3 ] 0 20 40 60 Depth [cm] SW-2 in-situ exp sigmoid Figure 4.2: Comparison of in-situ bulk density ( b ) proles and two prole represen- tations: exponential and sigmoid functions. The numbers in the parentheses in the legend are the RMSE of each modeled prole. the corresponding bulk density b values for the proles in Fig. 4.1 are shown in Fig. 4.2. Note that these bulk density proles are not obtained from tting the in-situ bulk density data; they are determined by the tted OM proles and the soil parameteri- zation. The results in 4.2 shows that the in-situ bulk density proles also possess the S-shaped characteristics, but with an opposite trend as their values increasing with depths. The bulk density proles calculated from the sigmoid OM proles still out- performs the ones calculated from the exponential OM proles. The errors between the calculated and in-situ proles could come from the mist of the OM proles or from the prediction errors of the soil parametrization model. 72 0 0.5 1 1.5 0 20 40 60 Depth [cm] HV-1 in situ model (0.083) 0 0.5 1 1.5 0 20 40 60 IC-1 in situ model (0.125) 0 0.5 1 1.5 0 20 40 60 IC-2 in situ model (0.007) 0 0.5 1 1.5 S w [m 3 /m 3 ] 0 20 40 60 Depth [cm] IM-1 in situ model (0.057) 0 0.5 1 1.5 S w [m 3 /m 3 ] 0 20 40 60 SW-1 in situ model (0.174) 0 0.5 1 1.5 S w [m 3 /m 3 ] 0 20 40 60 SW-2 in situ model (0.014) Figure 4.3: Comparison of in-situ saturation fraction (S w ) proles and the prole representation in (4.4). The numbers in the parentheses in the legend are the RMSE of each tted prole. 4.2.2 Soil Moisture Proles The volumetric soil moisture prole ((z)) is a product of proles of porosity () and saturation fraction (S w ): (z) =(z)S w (z) (4.3) Because the porosity only depends on the soil mineral texture and OM content, we will focus on the saturation fraction here for describing soil moisture proles. The in- situ saturation proles measured at the soil pits are shown in Fig. 4.3. The bottoms of these proles are mostly close to saturation (S w = 1) due to the underlying permafrost table impeding the water drainage. We then propose to assume the active layer soil prole is fully saturated from the bottom, and has a quadratic transition near the 73 surface. S w (z) = 8 > < > : S w0 1 WTD 2 (zWTD) 2 + 1; if zWTD 1; if z >WTD: (4.4) where S w0 is the saturation fraction at the surface and the boundary between satu- rated (S w = 1) and unsaturated (S w < 1) soils is referred to as the water table depth (WTD). The proposed prole function in (4.4) shows a good t with the in-situ data and the RMSE are generally smaller than 0.09 m 3 /m 3 . Larger errors might occur when the bottom permafrost table is not fully impermeable (e.g., SW-1 soil pit) or the lateral water ows are allowed in the lower mineral soil layers (e.g., IC-1 soil pit). This will be addressed in the future model improvements to better describe the subsurface hydrological processes. 4.3 Radar Retrieval Process 4.3.1 Forward Model To simulate the radar backscatter from permafrost soils, we need to model the active layer and underlying permafrost as a soil dielectric prole. To this end, we have developed the harmonized parametrization as described in Chapter III to estimate the soil dielectric constant for a given OM content and saturation fraction. With the OM and saturation prole functions proposed in (4.2) and (4.4), we now can express the soil dielectric prole using only four input parameters (A, M, S w0 , WTD). Fig. 4.4 shows an example of proles of bulk density, porosity, soil moisture, and soil dielectric constant from the assumed OM andS w proles. The prole parameters for OM prole are set to (A, M) = (0.7, 0.2 m); the prole parameters for saturation prole are set to (S w0 , WTD) = (0.2, 0.2 m). With this combined prole parametrization scheme, we will be able to estimate complex proles of geophysical parameters from radar, rather than just discrete soil dielectric layers as in the baseline retrievals in Chapter 74 0 50 100 OM [%] 0 20 40 60 Depth [cm] 0 0.5 1 [m 3 /m 3 ] 0 20 40 60 Depth [cm] 0 0.5 1 1.5 b [g/cm 3 ] 0 20 40 60 Depth [cm] 0 0.5 1 S w [m 3 /m 3 ] 0 20 40 60 Depth [cm] 0 0.5 1 [m 3 /m 3 ] 0 20 40 60 Depth [cm] 0 20 40 60 80 0 20 40 60 Depth [cm] Figure 4.4: Example of prole representation of soil compositional properties (or- ganic matter content (OM), porosity (), bulk density ( b )), saturation fraction (S w ) and soil moisture (), and the corresponding soil dielectric prole (). The prole parameters are set to (A,M,S w0 ,WTD) = (0.7, 0.2 m, 0.2, 0.2 m) for this example. II, while limiting the number of unknowns in the inverse problem. Figure 4.5: Geometry of a multilayer dielectric structure (left) representing the active layer and underlying permafrost, where (z) are discretized with a spacing of 1 cm and z ALT is the ALT. 75 In the radar scattering model, the soil dielectric prole is discretized as a layered dielectric structure with a layer thickness of 1 cm. The geometry of a multilayer dielectric structure representing the active layer and underlying permafrost is shown in Fig. 4.5. For the dielectric layers within the thawed active layer (z < z ALT ), the above mentioned prole parametrization model is applied, whereas for permafrost (zz ALT ) we assign a dielectric constant of 5 to represent the frozen soils. The top of the dielectric structure is assumed to be a rough surface represented by root-mean- square (rms) height h and correlation length l c , where the ratio of l c =h is assumed to be 10 in this work. A small perturbation method (SPM) up to the rst order developed in [33] is then used to calculate the radar backscatter from the layered dielectric structure with the top rough surface. 4.3.2 Inverse Problem The inverse problem is posed as an optimization problem by minimizing a cost function dened as L(X(t 1 ;t 2 ;:::;t Nt )) = 1 N X pp;tn j 0 pq (X(t n ))d pq j 2 ; (4.5) where X(t n ) = [A;M;S w0 (t n );WTD(t n );z ALT ;h] is the vector of unknowns at time snapshot t n . The values of 0 pp and d pp are the calculated and measured co-polarized (pp = hh, vv) P-band radar backscattering coecients, respectively. A two-element time-series of 0 observations, one taken in mid-late August (t 1 ) and one taken in early October (t 2 ) are considered in the retrieval of soil proles in this work. As shown in Fig. 4.6, we assume that there are no major human or wildre activities between t 1 and t 2 , hence OM prole (A and M) and surface roughness (h) can be considered as time-invariant parameters. Moreover, since the upward freezing front has not risen signicantly from the permafrost table and the active layer is not entirely frozen att 2 , 76 Figure 4.6: Prole models used in the retrieval of active layer soil properties. The organic matter prole is assumed to be time-invariant. For saturation proles, the surface saturation fraction (S w0 ) and water table depth (WTD) can vary between mid-late August and early October, whereas the bottom of the active layer remains fully saturated. ALT can still be observable andz ALT can also be assumed constant betweent 1 andt 2 . There are thus four input 0 observations ( 0 hh (t 1 ), 0 vv (t 1 ), 0 hh (t 2 ), 0 vv (t 2 )) and eight unknowns consisting ofA,M,h,z ALT S w0 (t 1 ),S w0 (t 2 ),WTD(t 1 ), andWTD(t 2 ). The search of the global minimum of the cost function L(X(t 1 ;t 2 )) is performed by using the simulated annealing (SA) method, which is a global optimizer using a randomized search method. 4.4 Results 4.4.1 Synthetic Data We evaluate the OM and soil moisture prole retrievals with Monte Carlo simula- tions for dierent values of M parameter (0.05 m, 0.1 m, 0.2 m), which can indicate the organic layer thickness for each case. The true prole parameters are set to (A, S w0 (t 1 ), WTD(t 1 ), S w0 (t 2 ), WTD(t 2 )) = (0.8, 0.3, 0.2 m, 0.2, 0.3 m) and the ALT is assumed to be 0.5 m. The synthetic radar data are generated from SPM 77 Figure 4.7: Comparison of true and estimated proles of organic matter, porosity, saturation fraction, and soil moisture for the case of M = 0:05m. Other prole parameters are set to (A, S w0 (t 1 ), WTD(t 1 ), S w0 (t 2 ), WTD(t 2 )) = (0.8, 0.3, 0.2 m, 0.2, 0.3 m). The shaded regions represent one standard deviation around the mean estimated values for each parameter. 78 Figure 4.8: Comparison of true and estimated proles of organic matter, porosity, saturation fraction, and soil moisture for the case ofM = 0:1m. Other prole param- eters are set to (A,S w0 (t 1 ),WTD(t 1 ),S w0 (t 2 ),WTD(t 2 )) = (0.8, 0.3, 0.2 m, 0.2, 0.3 m). The shaded regions represent one standard deviation around the mean estimated values for each parameter. 79 Figure 4.9: Comparison of true and estimated proles of organic matter, porosity, saturation fraction, and soil moisture for the case ofM = 0:2m. Other prole param- eters are set to (A,S w0 (t 1 ),WTD(t 1 ),S w0 (t 2 ),WTD(t 2 )) = (0.8, 0.3, 0.2 m, 0.2, 0.3 m). The shaded regions represent one standard deviation around the mean estimated values for each parameter. 80 simulations without adding any noise as the problem is already ill-posed for having more unknowns than observation points and we would like to focus on examining the sensitivity of radar backscatter to the model parameters and the associated retrieval uncertainty. The prole retrieval results are shown in Fig. 4.7-4.9. Due to the am- biguity of soil dielectric constant for dierent OM and soil moisture levels, as shown in Fig. 3.4, the absolute OM values are not well captured by the retrievals. This could also be due to the fact that the sensitivity of radar backscatter to dielectric changes will degrade when large soil dielectric constant value is present. Although the absolute OM levels are not well estimated, the retrievals can still capture the organic layer thickness variations as the parameter M changes. The soil dielectric constant is more sensitive to the changes in saturation fraction (or soil water content) than to the changes in OM content, so the retrieved and true saturation proles have a much better agreement with smaller retrieval uncertainty. Porosity proles are derived from the OM proles so the retrieval accuracy depends on the one of retrieved OM proles. Since soil moisture is the product of porosity and saturation fraction, some details in the prole variations might not be well captured due to the errors contributed from the retrieved porosity prole. Note that only the polarimetric radar backscatter measurements are considered in the retrievals here, and the sensitivity of polarimetric radar data to the desired prole parameters has limits. The retrieval accuracy and uncertainty of OM proles can be improved in the future by introducing more observations by dierent instruments (e.g., passive radiometer does not lose sensitivity when high dielectric values are present) or by integrating other soil physical models that describe dierent physical aspects of soils (e.g., soil hydraulic and thermal properties). 81 4.4.2 AirMOSS Data in 2017 ABoVE Airborne Campaign The proposed forward model and retrieval process are applied to retrieve the soil moisture and OM proles at the Soil moisture Sensing Controller and oPtimal Estimator (SoilSCAPE) site in Happy Valley, Alaska (69°9'11"N, 148°50'39"W). The time-series P-band AirMOSS data used in the retrieval were acquired as part of 2017 NASA ABoVE Airborne Campaign (AAC) on August 13 and October 9, 2017. The retrieved soil dielectric proles using the baseline algorithm (three-layer structure) in Chapter II and the proposed prole representation are shown in Fig. 4.10. The in-situ soil dielectric measurements taken by four METER TEROS-12 soil sensors deployed at the depths of 5, 15, 30, 45 cm are also shown for comparisons. The retrieved soil dielectric proles are generally overestimated but still show a similar prole shape as the in-situ measurements. The overestimation is likely because the soil porosity calculated from the proposed soil parametrization is too large for soils with moderate-low SOM content, resulting in higher values for saturated soil moisture and dielectric constant. Also, the TEROS-12 sensors are operated at 70 MHz, which is very dierent than the P-band frequency (430 MHz) used in this work and can be another factor causing the discrepancy. The retrieved soil moisture and OM proles are shown in Fig. 4.11. The soil moisture considered here represents the unfrozen water content in the soil. The results show that the unfrozen water content near the surface decreases as the top soil is getting frozen when entering the freezing season. The retrieved OM prole shows a thick organic layer (dened as the depth at which the OM content has decreased to 20 %) with a thickness of 20 cm (the retrieved M = 17:8 cm), which is slightly thinner than the in-situ prole sampled around the pixel, and is close to the site average of 24 cm. The parameters used in the soil parametrization and dielectric mixing model will need to be further calibrated against existing in-situ soil measurements and soil databases before being applied to larger areas in the Arctic-boreal regions. 82 0 20 40 60 80 0 0.1 0.2 0.3 0.4 0.5 0.6 Depth [m] 2017-08-13 Baseline Proposed in-situ 0 20 40 60 80 0 0.1 0.2 0.3 0.4 0.5 0.6 Depth [m] 2017-10-09 Baseline Proposed In-Situ Figure 4.10: Comparison of retrieved and in-situ soil dielectric proles at the SoilSCAPE site in Happy Valley, Alaska. 0 0.5 1 [m 3 /m 3 ] 0 0.1 0.2 0.3 0.4 0.5 0.6 Depth [m] 2017-08-13 2017-10-09 0 0.5 1 OM [%] 0 0.1 0.2 0.3 0.4 0.5 0.6 Depth [m] retrieved in-situ Figure 4.11: Retrieved soil moisture and OM proles at the SoilSCAPE site in Happy Valley, Alaska. 83 4.5 Conclusion In this chapter, we advance the baseline algorithm introduced in Chapter II by in- corporating the harmonized soil parametrization scheme described in Chapter III and a realistic prole representation for active layer soils in Section 4.2 into the retrieval. Now we are able to retrieve much more complex active layer soil prole characteristics with only one more unknown introduced in the inverse problem. The retrieved soil moisture and organic matter proles will have a broader audience in the science com- munities as they can directly use the retrieved products in their research activities. Although the approach in this chapter will still suer from the limited sensing depths for ALT retrievals as no additional information or observations related to the deep soil features are provided, the current parametrization scheme and retrieval process can now easily be integrated with other soil physics models or extended wherever new measurements are available. 84 CHAPTER V Joint Retrieval of Active Layer Proles Using PolSAR Backscatter and InSAR Deformation 5.1 Introduction The interferometric synthetic aperture radar (InSAR) is a geodetic technique that uses the dierences in the observed phase of two or more synthetic aperture radar (SAR) images to measure surface deformation or elevation at high spatial resolu- tion [57]. The application of the InSAR technique for the active layer thickness (ALT) estimation, also known as the remotely sensed active layer thickness (ReSALT) algo- rithm, uses the surface deformation to infer ALT [19] [6]. Frost heave happens when the soil water in the active layer freezes in fall and reaches the maximum heave when the entire soil column is frozen. When the active layer thaws in summer the ground surface subsides and the subsidence between the thaw onset and the maximum thaw can be used to estimate ALT. The ReSALT algorithm has been used to monitor the active layer and permafrost dynamics, as well as the thermokarst processes across the Arctic domain [19] [58] [6]. The ReSALT algorithm is well suitable for fast assessment of ALT distributions and thermokarst activities. However, there are some limitations inherit from the assumptions made for active layer soils in its seasonal subsidence model, e.g., the soil 85 prole is assumed to be fully saturated at all times, which might not hold for most areas as lateral water ows, evapotranspiration, and precipitation can all uctuate the soil moisture prole in the active layer. In this chapter, we propose to combine the polarimetric SAR (PolSAR) and InSAR techniques to better estimate the active layer properties, including the soil moisture proles and ALT. There is a great potential of exploiting the complementary strengths in both techniques: PolSAR is sensitive to changes in soil moisture but is limited by the radar sensing depth, whereas InSAR can accurately estimate surface subsidence due to the active layer thaw but large uncertainty is expected without soil moisture prole information. 5.2 Forward Model For a monostatic radar conguration, the incident (E i ) and scattered (E s ) electric elds can related by using a complex scattering matrix (S): 2 6 4 E h E v 3 7 5 s = ( e jkr r ) 2 6 4 S hh S hv S vh S vv 3 7 5 2 6 4 E h E v 3 7 5 i (5.1) where k is the wave number and r is the distance from the radar to the target. The subscripts h and v denote the horizontal and vertical components, respectively. For PolSAR, the radar backscattering cross section of a target is determined as: 0 pq = 4jS pq j 2 (5.2) which is in general a function of the target's geometrical and dielectric properties. For InSAR, the phase term ine jkr is used instead. The two-way phase delay between the radar and the target is given as =2kr = 4 r (5.3) 86 Figure 5.1: Schematic diagrams of active layer conditions throughout a thaw season. The shaded blue regions represent the part of the active layer that is still frozen. is the seasonal subsidence that occurs when the maximum thaw of the active layer is reached. Considering the repeat-pass mode where the motion deviations of each acquisition are assumed zero or decorrelated, the InSAR phase of two SAR images att 1 andt 2 is insar =(t 2 )(t 1 ) =2k [r(t 2 )r(t 1 )] =2k (5.4) where is the surface deformation. Note that is only measured modulo 2, phase unwrapping is needed to add proper multiples of 2 to insar [59,60]. Also, the phase changes between S pq (t 1 ) andS pq (t 2 ) are commonly ignored although it has been well known that the soil moisture changes are observable in InSAR phase measurements [61,62]. 5.2.1 ReSALT Subsidence Model The seasonal subsidence due to the phase change of soil water in the active layer is illustrated in Fig. 5.1. The ground surface continuously subsides from the onset of thawing as ice in the soil is melted to liquid water until reaching the maximum thaw 87 depth in late summer. In the ReSALT subsidence model, the following assumptions are made [19]: 1. The surface subsidence is caused purely by the phase change of soil water in the active layer. The expansion and contraction of soils within the active layer are assumed to be minimal. 2. The active layer is fully saturated for the entire soil column at all times. 3. The active layer is fully frozen in winter and before the thaw onset, no unfrozen water content exists in the soil. With these assumptions, the surface subsidence () at a certain thaw depth (TD) can be expressed as = TD Z 0 (z)S w (z) w i i dz (5.5) where is the soil porosity, S w is the saturation fraction, w and i are the densities of liquid water and ice, respectively. Porosity is dened as the fraction of the volume of voids (V void ) in the total soil volume (V soil ): = V void V soil (5.6) Saturation fraction is dened as the volumetric fraction of soil water occupied in the total porosity, i.e., S w = V water V void (5.7) Multiplying (5.6) and (5.7) will give the soil moisture (): = V water V soil = V void V soil V water V void =S w (5.8) 88 Therefore, the equation in (5.5) is essentially an integral of soil moisture prole mul- tiplied by a factor due to phase change. = TD Z 0 (z) w i i dz (5.9) The typical values for w and i are 1.0 and 0.917 g=cm 3 , respectively, which means that the factor in (5.9) represents a change (decrease) in volume of soil water by about 9.05% when the ground ice melts. Although the active layer is assumed fully saturated at all times (S w (z) = 1), the ReSALT subsidence model still requires the soil porosity prole(z) before calculating the integral in (5.5). The porosity prole is estimated by a mixed soil model that assumes the organic matter content and thus the porosity decrease exponentially with soil depth from the surface pure organic soil to deep pure mineral soil [19]. Typical parameters values for permafrost soils in the Arctic region are used and the same porosity prole is applied to every pixel where the ReSALT algorithm is performed, which could introduce uncertainty in the retrieved ALT values because the soil organic matter (SOC) is known to be heterogeneous in its spatial and vertical prole distribution. This single porosity prole assumption may be overcome by exploiting available SOC maps to properly parametrize the porosity prole at each pixel. The large uncertainty and coarse vertical resolution of these soil maps, however, pose a challenge of applying them to the ReSALT algorithm for a better retrieval performance. To address the potential errors and uncertainties resulted from the full-saturation and single-porosity-prole assumptions, we adopt the permafrost soil prole model developed in Chapter III and Chapter IV to provide a basis of modeling saturation and porosity proles simultaneously. This also allows us to be consistent in terms of soil modeling when integrating the PolSAR forward backscatter model in the subsequent 89 Figure 5.2: Surface deformation in units of meters with respect to varying water saturation proles (S w0 , WTD) and thaw depths (TD). Parameters of organic matter prole are xed to A = 0.8 and M = 0.15 m. sections. With the permafrost soil prole model, the subsidence will now take ve parameters as input: A andM for organic matter prole,S w0 and WTD for saturation prole, and TD for thaw depth. Fig. 5.2 shows the values of with respect to varying S w0 , WTD, and TD when organic matter prole is xed to A = 0.8 and M = 0.15 m. is most sensitive to TD; the typical values for ALT are from 1 cm to 5 cm, approximately corresponding to ALT of 30 cm to 1 m. For a given TD, is generally smaller when the soil prole is dry as fewer water content participate in the phase change. When full-saturation assumption is negated, the ambiguity of with respect to TD becomes prominent as multiple saturation proles could eectively result in similar values. Fig. 5.3 shows the values of with respect to varying A and M when saturation prole (S w0 , WTD) and TD are xed. is generally more sensitive to organic matter prole changes whenA andM are suciently large than 0.2 and 10 cm, respectively, which correspond to a thicker and highly organic surface layer that has larger pore space to host water and can thus contribute more to the surface deformation from the water phase change. As observed from Fig. 5.2 and Fig. 5.3, is still most sensitive to TD and not as much to other parameters. However, the errors in saturation and organic matter proles (or equivalently the soil moisture prole) can easily introduce 90 0 0.2 0.4 0.6 0.8 1 A [-] 0.1 0.2 0.3 0.4 0.5 M [m] S w0 = 0.3, WTD = 0.2 m, TD = 0.5 m 0.02 0.025 0.03 0.035 0.04 (a) 0 0.2 0.4 0.6 0.8 1 A [-] 0.1 0.2 0.3 0.4 0.5 M [m] S w0 = 0.8, WTD = 0.1 m, TD = 0.5 m 0.02 0.025 0.03 0.035 0.04 (b) Figure 5.3: Surface deformation in units of meters with respect to varying organic matter proles (A, M). large uncertainty in TD estimates. Another aspect to look at the weaker sensitivity of to soil moisture proles is that when the measurements are contaminated by noise or processing artifacts, the soil moisture retrievals will not be impacted much. 5.2.2 PolSAR Backscatter Model To estimate the PolSAR radar backscatter from the same target scene observed by the InSAR techniques, we exploit the same soil prole parametrization model as described above for consistency. The PolSAR backscatter model requires modeling the soil prole as a layered dielectric structure, which can be achieved by the methods described in Chapter IV. Compared to , radar backscatter 0 will take one more parameter as input: the roughness h of the top soil surface. Fig. 5.4 shows the simulated values of 0 HH with the same parameter setting as in Fig. 5.2 except the additional parameter h is set to 2 cm. The oscillating behavior due to the multiple scattering within the thawed soils is observed. As expected, radar backscatter is most sensitive to the water saturation near the surface, and is ambiguous with respect to depth variations of WTD and TD. When the soil thaws deeper, radar backscatter also loses the sensitivity to TD as the signal scattered back from the permafrost table 91 Figure 5.4: P-band HH radar backscatter 0 HH in units of dB with respect to varying water saturation proles (S w0 , WTD) and thaw depths (TD). Parameters of organic matter prole are xed to A = 0.8 and M = 0.15 m, and the surface roughness (h) is assumed to be 2 cm in this example. becomes weaker. 0 0.2 0.4 0.6 0.8 1 A [-] 0.1 0.2 0.3 0.4 0.5 M [m] 0 HH -22 -21.5 -21 -20.5 -20 -19.5 -19 (a) 0 0.2 0.4 0.6 0.8 1 A [-] 0.1 0.2 0.3 0.4 0.5 M [m] 0 VV -16.5 -16 -15.5 -15 -14.5 -14 -13.5 (b) Figure 5.5: P-band radar backscatter 0 in units of meters with respect to varying organic matter proles (A,M). (S w0 , WTD, TD, h) is set to (0.5, 0.2 m, 0.5 m, 0.02 m). Fig. 5.5 shows the values of 0 with respect to varyingA andM when saturation prole (S w0 , WTD), TD, and h are xed. When the organic layer thickness M is smaller than WTD (0.2 m in this case), 0 values do not have much sensitivity to either A or M because the organic-to-mineral transitions are within the unsaturated zones, where the porosity changes due to organic matter prole cannot fully translate 92 to changes in overall soil moisture and eective soil dielectric constant. When part of the organic layer is immersed in the saturated zone, the sensitivity of 0 with respect to organic proles becomes more prominent. 5.3 Combined PolSAR-InSAR Retrieval In the ReSALT algorithm, a simple bisection method is used to converge on an estimated TD that matches the modeled model and the observed surface deformation InSAR within an error criterion (e.g., 1 mm). When InSAR represents the seasonal subsidence between the thaw onset and maximum thaw, the estimated TD will be the ALT. The InSAR -TD relation is 1-to-1 if the above mentioned assumptions are met and can be used for a quick, rst-order assessment of the thaw depth distribution. However, those assumptions do not hold all the time and special treatments are often required for a specic site or dataset for optimal retrieval performance. The key to integrate the PolSAR and InSAR approach is to have the same set of soil parameters so that the (z) and S w (z) proles will be identical in both for- ward models and can be estimated simultaneously in the retrieval process. The cost function L involving only the InSAR deformation observations can be dened as L (X) =j model (X) InSAR j (5.10) where model and InSAR are the surface deformation from the forward model in (5.5) and InSAR observations, respectively. X denotes the vector of unknowns. The cost function L involving only the PolSAR backscatter observations is given by L (X) = 1 2 X pp=HH;VV j 0 pp;model (X) pp;PolSAR j 2 (5.11) where 0 pp;model and pp;PolSAR are the pp-polarized backscatter values from the radar 93 (a) (b) (c) Figure 5.6: Cost function values of a active layer soil prole with varying water saturation characteristicsS w0 (-) and WTD (m) and thaw depths TD (m). (a) InSAR only (L ). (b) PolSAR only (L ). (c) Joint PolSAR-InSAR (L joint ). The true values of unknowns are set to S w0 = 0.5, WTD = 0.1 m, TD = 0.5 m, h = 0.02 m. forward backscattering model and PolSAR observations, respectively. Finally, the joint cost function L joint involving both InSAR and PolSAR is given by L joint (X) = 1 L (X) + 1 () 2 L (X) (5.12) where is the uncertainty in measurements and is the PolSAR radiometric calibration accuracy. is a regularization term that controls the weights of the cost function values towards InSAR or PolSAR terms. In this work, is set to 0.5 cm [19] and is set to 0.5 dB [24]. We consider equal contributions from InSAR 94 and PolSAR cost function terms so is set to 1. There are six unknowns to be solved in X, namely X = [A;M;S w0 ;WTD;TD], with the three observation points (, 0 hh , 0 vv ). An example of the cost function values in (5.10), (5.11), and (5.12) is shown in Fig. 5.6 with the same parameter setting as in Fig. 5.2. InSAR only approach (L ) has large uncertainties for estimating TD if inaccurate soil moisture information is provided, as shown in Fig. 5.6(a). For PolSAR only approach (L ) in Fig. 5.6(b), there will still be an ambiguity issue in the depth retrievals (WTD and TD) if no time-series observations are available. With L and L resolving the layer depths and soil moisture respectively, the joint PolSAR-InSAR approach (L joint ) can solve both inverse problems well, as the solutions converging to the points close to the true values in Fig. 5.6(c). 5.4 Results 5.4.1 Synthetic Data We rst evaluate the organic matter and soil moisture prole retrievals with Monte Carlo simulations under a noise-free scenario. The true parameters are set to (A, M S w0 , WTD, TD) = (0.8, 0.2 m, 0.3, 0.2 m, 0.5 m). The retrieved proles are shown in Fig. 5.7. Due to the ambiguity of and 0 with respect toA andM as observed in Fig. 5.5 and Fig. 5.3, the organic matter prole is generally underestimated with large uncertainty. As a result of having InSAR and PolSAR techniques resolving the layer depths and the near-surface dielectric proles, respectively, the retrieved and true saturation fraction proles have a much better agreement and smaller retrieval un- certainty. Although the soil moisture retrievals are degraded by the relatively poorly retrieved organic matter prole, the retrieved soil moisture can still at least capture the near-surface soil moisture level and the overall prole shape in the subsurface. 95 Figure 5.7: Comparison of true and estimated proles of organic matter, porosity, saturation fraction, and soil moisture. True parameters are set to (A,M S w0 ,WTD, TD) = (0.8, 0.2 m, 0.3, 0.2 m, 0.5 m). The shaded regions represent one standard deviation around the mean estimated values for each parameter. The soil moisture prole retrievals can be further improved if more prior information about the true organic matter prole becomes available. The results of the joint PolSAR-InSAR retrieval for thaw depths ranging from 40 cm to 90 cm are shown in Fig. 5.8. As described in Chapter II, the P-band sensing depth for PolSAR approach is about 50 cm; the joint retrieval can apparently resolve a much deeper thaw depth than that. The overestimation of the TD estimates is most likely due to the underestimated organic matter and porosity proles, which 96 Figure 5.8: Comparison of true and estimated thaw depth (TD). Other true param- eters are set to (A, M S w0 , WTD) = (0.8, 0.2 m, 0.3, 0.2 m, 0.5 m). The error bars represent one standard deviation around the mean estimated TD values. makes the InSAR subsidence model to integrate into deeper soil depths to match the observed surface deformation. The retrieval uncertainty is about 10 cm in the TD estimates, which is similar to the one of the baseline retrievals developed in Chapter II. 5.4.2 UAVSAR and AirMOSS Data in 2017 ABoVE Airborne Campaign The joint retrieval of soil moisture proles and ALT is performed using the P-band Airborne Microwave Observatory of Subcanopy and Subsurface (AirMOSS) PolSAR data and L-band Uninhabited Aerial Vehicle Synthetic Aperture Radar (UAVSAR) InSAR interferograms acquired during the 2017 (Arctic-Boreal Vulnerability Exper- iment) ABoVE airborne campaigns. Two Alaskan sites where wild res had oc- curred are selected for demonstrating the joint retrievals here: the Yukon{Kuskokwim (YK) Delta on the west coast of Alaska and Anaktuvuk River on the North Slope in Alaska. The acquisition dates for the L-band UAVSAR interferograms are June 17 and September 17, 2017 for YK Delta, and June 21 and September 16, 2017 for 97 Figure 5.9: Jointly retrieved maps of ALT and surface water saturation fraction S w0 using P-band AirMOSS PolSAR data and L-band UAVSAR InSAR deformation acquired over the Yukon-Kuskokwim (YK) Delta in 2017. Anaktuvuk River. The P-band PolSAR data were acquired on August 14 and August 13 for YK Delta and Anaktuvuk River, respectively. Fig. 5.9 shows the maps of ALT and water saturation at the surface S w0 retrieved jointly from the 2017 AirMOSS and UAVSAR data over the Yukon-Kuskokwim (YK) Delta. When wildre occurs in tundra landscapes in the Arctic region, the re removes a portion of surface organic layer, which has an insulating functionality regulating the thermal uxes exchanged between the atmosphere and subsurface active layer and permafrost. The ALT in the areas within the re scar are therefore deeper than the unburned areas, shown in Fig. 5.9, as the active layer is more prone to the summer thawing. Moreover, after the surface vegetation and organic layer are removed by the wildre, the evapotranspira- tion from the ground surface to the atmosphere becomes weaker, making the surface 98 Figure 5.10: Jointly retrieved maps of ALT and surface water saturation fraction S w0 using P-band AirMOSS PolSAR data and L-band UAVSAR InSAR deformation acquired over the Anaktuvuk River in 2017. soil more saturated compared to the unburned areas. The retrieval results in Fig. 5.9 show that the joint PolSAR-InSAR retrieval is able to capture both ALT and surface soil moisture conditions related to wildre disturbance. Fig. 5.10 shows the maps of the retrieved ALT and water saturation for the ight lines own in 2017 over the Anaktuvuk River. The L-band interferogram was severely contaminated by the atmospheric noise, and as result phase ramp artifacts are clearly observed in the ALT retrievals. As the surface subsidence is the critical and dominant source of information for determining ALT in the joint retrieval, the artifacts related to InSAR deformation will show up prominently. The surface water saturation, however, is still well resolved by the P-band PolSAR backscatter. If more time-series scenes are available (e.g., through satellite missions or more frequent 99 airborne deployments), InSAR data will have higher signal-to-noise ratio (SNR) for avoiding atmospheric noise eects and PolSAR backscatter will have better capability to remove the model ambiguity when retrieving complex subsurface features. 5.5 Conclusion In this chapter, we further advance the PolSAR-based algorithm by incorporating InSAR techniques into the ALT retrievals. Both approaches are harmonized in the sense that the same soil prole parametrization kernel is shared in both models. The InSAR-based approach (the ReSALT algorithm) provides an accurate method to relate the seasonal subsidence due to phase change of soil water to the ALT. However, several assumptions made for the InSAR method might not hold for most areas in permafrost regions, including the ones related to full saturation and fully frozen state in the active layer. The PolSAR-based approach, on the other hand, is limited by the radar sensing depth and the ambiguity of ALT estimates, but is sensitive to the changes in soil moisture proles within the active layer. The complementary nature of both methods enables an accurate joint retrieval algorithm for ALT and soil moisture prole estimations. The results show that the proposed joint retrieval scheme is able to capture the ALT and soil moisture dynamics that respond to the wildre disturbance. 100 CHAPTER VI Conclusion and Future Work 6.1 Conclusion This dissertation investigates radar remote sensing techniques for monitoring per- mafrost landscapes and estimating active layer soil properties, particularly the active layer thickness (ALT) and soil moisture prole distributions in the Arctic-boreal re- gion. To estimate active layer properties from radar, we rst propose a retrieval algo- rithm using time-series P-band polarimetric synthetic aperture radar (PolSAR) ob- servations. The active layer and underlying permafrost are modeled as a three-layer dielectric structure with the layer dielectric constants representing the soil moisture and freeze-thaw state of the layer. To resolve the ambiguity of the retrieved layer thicknesses, an approach of nding the largest possible depths (LPD) is combined with time-series observations, where the ALT is assumed time-invariant between the maximum thaw and before the upward freezing front rises signicantly from the per- mafrost table. The LPD-assisted time-series retrieval algorithm is applied to the radar data acquired by the Airborne Microwave Observatory of Subcanopy and Sub- surface (AirMOSS) P-band radar in August and October 2014 and 2015 over the Alaska North Slope. The retrieval results show that the active layer properties are strongly in uenced by the land cover types at the regional scale, and not as much 101 by the North-South temperature gradient. The LPD-assisted retrieval process with time-series approach is treated as a baseline algorithm in this dissertation, and it provides insights on the advancements that follow. Depending on the organic matter content and degree of decomposition, organic soils can be classied as bric (least decomposed), hemic (moderately decomposed), sapric (highly decomposed), and mineral soils. These soil classes possess dierent hydraulic properties, in particular the capability to adsorb water around the surface of soil particles, which can greatly aect the soil dielectric properties. We propose a method to harmonize the soil parameters that are used in both soil dielectric and hydrologic models so that the bound-to-free water ratio can be derived from the same parameter set. A semi-theoretical dielectric mixing model for active layer soils is developed and serve as a bridge between the physics of electromagnetics and other soil sciences. With this development, the retrieval process will be able to directly retrieve more informative geophysical properties such as soil moisture and organic matter proles, rather than just dielectric proles as in the baseline algorithm. This harmonized parametrization scheme is also a step towards an synergistic framework using the joint-physics of radar scattering and soil process modeling. The outcome of the harmonized soil parametrization enables a more realistic pro- le representation of active layer soils, which in turn allows us to spend the number of unknowns more eciently on the prole variables in the retrieval process, compared to the layer dielectric constants and thicknesses in the baseline retrievals. With the organic matter and surface soil water saturation proles parametrized as a sigmoid and a quadratic function, respectively, we extend the baseline algorithm to estimate ALT and proles of soil moisture and organic matter from time-series P-band PolSAR observations. To further advance the radar retrievals, we incorporate the surface deformation derived from interferometric SAR (InSAR) into the retrieval process. The seasonal 102 subsidence due to the volumetric changes of soil water during the thawing progress is used to infer the thaw depth, or the ALT when the subsidence is measured between the thaw onset and maximum thaw. However, without the knowledge of soil moisture proles at the thaw onset and maximum thaw, the ALT estimates by the InSAR technique can be well biased. To overcome this, we rst use the proposed prole parametrization of active layer soils to integrate the PolSAR-based forward model with the InSAR subsidence model. A joint cost function is formed to exploit the complementary strengths of both methods to improve the ALT and soil moisture retrievals. The retrieval results show that the joint estimation scheme is able to capture the ALT and soil moisture dynamics that respond to the wildre disturbance. 6.2 Future Work 6.2.1 Implementing a joint radar scattering and soil process inversion framework We will utilize the proposed radar retrieval algorithm jointly with a permafrost soil model which considers both subsurface thermal dynamics and hydrologic processes. The integrated framework for joint radar and soil model inversions is illustrated in Fig. 6.1. The radar retrieval and soil process model will be run in an iterative manner explained as follows. For each snapshot of L- or P-band SAR data, we will perform the radar retrieval algorithm to get an initial estimates of soil moisture and organic matter proles, ALT, and freeze/thaw state. Since the existing soil organic carbon (SOC) databases for arctic and boreal regions still have large uncertainties and lim- ited spatial coverage and resolution, the organic matter prole informed by radar retrievals can greatly improve the soil model simulations. The radar-retrieved soil moisture prole will also be used to drive the hydrologic model for the vertical water ows. The soil process model simulations can in turn inform radar retrievals the soil 103 Figure 6.1: Diagram of the joint modeling framework moisture and temperature proles for a regularization scheme to better constrain the underdetermined nature of the radar scattering model. The soil moisture proles com- puted from Richards' equation can help radar retrieval algorithm determine a more physically realistic prole function for retrieving the prole of saturation fraction. The model-simulated soil temperature proles can even greatly reduce the ambiguity of radar backscatter to thaw depths which are deeper than radar sensing depths. We will design and optimize this iterative inversion process to simultaneously improve the radar retrievals and soil model simulations, so that the data products that we deliver will be consistent in terms of the physics of electromagnetics, hydrology, and thermal dynamics. 6.2.2 Estimation of Vegetation Structure Using SAR and Lidar Due to the long wavelength of P-band radar, canopy height, tree radius, and stem density are the most important structural parameters of vegetation for accurate radar backscatter modeling. Waveform lidar can provide accurate canopy height estimates but has a very narrow swath compared to SAR. For the areas imaged by both instruments, we will investigate the relationship between the canopy height 104 Figure 6.2: Estimation of canopy heights for radar-only pixels in Yukon Flats ecore- gion in Alaska proles (CHP) derived from lidar and radar backscattering measurements, and select a subset of them as training data set. Other ancillary data such as land cover types and re disturbance history will also be used in developing a regression model to estimate the canopy heights for radar-only pixels. Moreover, if the major tree species are known for the areas of interest, the diameter at breast height (DBH) can also be obtained through allometric relations. With canopy height and tree radius estimated, we will explore the potentials of estimating stem density from a SAR-Lidar integrated inversion framework, or simply treat stem density as an unknown variable in the radar retrieval process. Outcomes of this future work will allow the radar retrieval algorithms developed in dissertation to be applied to the forested areas in the Arctic- boreal region. As a rst attempt, Random Forests, an extension of Classication and Regression Trees (CART) algorithm, is used to t a model between the lidar-derived canopy heights and a number of data layers including AirMOSS P-band radar backscatter (HH, HV, VV), topography, land cover types (NLCD), and re disturbance history. 105 The preliminary results performed for the Yukon Flats ecoregion in Alaska are shown in Fig. 6.2, with a root-mean-square error (RMSE) of 2.08 m, a bias of 0.03 m, and a R 2 -value of 0.59. We will investigate how to improve the estimation accuracy and integrate this estimation framework with the radar retrieval process. 106 BIBLIOGRAPHY [1] M.-k. Woo, Permafrost Hydrology. Berlin Heidelberg: Springer-Verlag, 2012. [2] B. K. Biskaborn, J.-P. Lanckman, H. Lantuit, K. Elger, D. A. Streletskiy, W. L. Cable, and V. E. Romanovsky, \The new database of the Global Terrestrial Network for Permafrost (GTN-P)," Earth System Science Data, vol. 7, no. 2, pp. 245{259, 2015. [3] G. J. M. De Lannoy, R. D. Koster, R. H. Reichle, S. P. P. Mahanama, and Q. Liu, \An updated treatment of soil texture and associated hydraulic properties in a global land modeling system," J. Adv. Model. Earth Syst., vol. 6, no. 4, pp. 957{979, Dec. 2014. [4] N. R. Peplinski, F. T. Ulaby, and M. C. Dobson, \Dielectric properties of soils in the 0.3-1.3-GHz range," IEEE Trans. Geosci. Remote Sens., vol. 33, no. 3, pp. 803{807, May 1995. [5] S. Bircher, F. Demontoux, S. Razandratsima, E. Zakharova, M. Drusch, J.-P. Wigneron, and Y. H. Kerr, \L-band relative permittivity of organic soil surface layersa new dataset of resonant cavity measurements and model evaluation," Remote Sensing, vol. 8, no. 12, p. 1024, Dec. 2016. [6] K. Schaefer, L. Liu, A. Parsekian, E. Jafarov, A. Chen, T. Zhang, A. Gus- meroli, S. Panda, H. A. Zebker, and T. Schaefer, \Remotely Sensed Active Layer Thickness (ReSALT) at Barrow, Alaska using interferometric synthetic aperture radar," Remote Sensing, vol. 7, no. 4, pp. 3735{3759, 2015. [7] IPCC,Climate Change 2007: Synthesis Report. Contribution of Working Groups I, II and III to the Fourth Assessment Report of the Intergovernmental Panel on Climate Change. Geneva, Switzerland: IPCC, 2017. [8] L. D. H. et al., \Evidence and implications of recent climate change in Northern Alaska and other Arctic regions," Climatic Change, vol. 7, no. 3, p. 251298, 2005. [9] D. A. Walker, G. J. Jia, H. E. Epstein, M. K. Raynolds, F. S. C. III, C. Copass, L. D. Hinzman, J. A. Knudson, H. A. Maier, G. J. Michaelson, F. Nelson, C. L. Ping, V. E. Romanovsky, and N. Shiklomanov, \Vegetation-soil-thaw-depth rela- tionships along a low-arctic bioclimate gradient, Alaska: synthesis of information 107 from the ATLAS studies," Permafrost and Periglacial Processes, vol. 14, no. 2, p. 103123, 2003. [10] S. A. Zimov, E. A. G. Schuur, and F. S. Chapin, \Permafrost and the global carbon budget," Science, vol. 312, no. 5780, p. 16121613, 2006. [11] J. Brown, K. M. Hinkel, and F. E. Nelson, \The Circumpolar Active Layer Moni- toring (CALM) program: Research designs and initial results," Polar Geography, vol. 24, no. 3, pp. 166{258, 2000. [12] S. L. Smith, S. A. Wolfe, D. W. Riseborough, and F. M. Nixon, \Active-layer characteristics and summer climatic indices, Mackenzie Valley, Northwest Terri- tories, Canada," Permafrost and Periglacial Processes, vol. 20, no. 2, pp. 201{ 220, 2009. [13] F. E. Nelson, N. I. Shiklomanov, G. R. Mueller, K. M. Hinkel, D. A. Walker, and J. G. Bockheim, \Estimating active-layer thickness over a large region: Kuparuk River Basin, Alaska, U.S.A." Arctic and Alpine Research, vol. 29, no. 4, pp. 367{378, 1997. [14] C. Gangodagamage, J. C. Rowland, S. S. Hubbard, S. P. Brumby, A. K. Liljedahl, H. Wainwright, C. J. Wilson, G. L. Altmann, B. Daon, J. Peterson, C. Ulrich, C. E. Tweedie, and S. D. Wullschleger, \Extrapolating active layer thickness measurements across Arctic polygonal terrain using LiDAR and NDVI data sets," Water Resources Research, vol. 50, no. 8, pp. 6339{6357, 2014. [15] U. Mishra and W. J. Riley, \Active-layer thickness across Alaska: Comparing observation-based estimates with CMIP5 earth system model predictions," Soil Science Society of America Journal, vol. 78, no. 3, pp. 894{902, 2014. [16] N. J. Pastick, M. T. Jorgenson, B. K. Wylie, S. J. Nield, K. D. Johnson, and A. O. Finley, \Distribution of near-surface permafrost in Alaska: Estimates of present and future conditions," Remote Sensing of Environment, vol. 168, pp. 301{315, 2015. [17] M. A. Rawlins, D. J. Nicolsky, K. C. McDonald, and V. E. Romanovsky, \Sim- ulating soil freeze/thaw dynamics with an improved pan-Arctic water balance model," Journal of Advances in Modeling Earth Systems, vol. 5, no. 4, pp. 659{ 675, 2013. [18] Y. Yi, J. S. Kimball, M. A. Rawlins, M. Moghaddam, and E. S. Euskirchen, \The role of snow cover aecting boreal-arctic soil freeze-thaw and carbon dynamics," Biogeosciences, vol. 12, no. 19, pp. 5811{5829, 2015. [19] L. Liu, K. Schaefer, T. Zhang, and J. Wahr, \Estimating 1992-2000 average active layer thickness on the Alaskan North Slope from remotely sensed surface subsidence," J. Geophys. Res.: Earth Surface, vol. 117, no. F1, p. F01005, 2012. 108 [20] P. C. Dubois, J. van Zyl, and T. Engman, \Measuring soil moisture with imaging radars," IEEE Trans. Geosci. Remote Sens., vol. 33, no. 4, pp. 915{926, Jul. 1995. [21] Y. Oh, K. Sarabandi, and F. T. Ulaby, \An empirical model and an inversion technique for radar scattering from bare soil surfaces," IEEE Trans. Geosci. Remote Sens., vol. 30, no. 2, pp. 370{381, Mar. 1992. [22] A. Tabatabaeenejad, M. Burgin, X. Duan, and M. Moghaddam, \P-band radar retrieval of subsurface soil moisture prole as a second-order polynomial: First AirMOSS results," IEEE Trans. Geosci. Remote Sens., vol. 53, no. 2, pp. 645{ 658, Feb. 2015. [23] V. L. Mironov, R. D. D. Roo, and I. V. Savin, \Temperature-dependable mi- crowave dielectric model for an Arctic soil," IEEE Trans. Geosci. Remote Sens., vol. 48, no. 6, pp. 2544{2556, Jun. 2010. [24] E. Chapin, B. P. Hawkins, L. Harcke, S. Hensley, Y. Lou, T. R. Michel, L. Mor- eira, R. J. Muellerschoen, J. G. Shimada, K. W. Tham, and M. C. Tope, \Im- proved absolute radiometric calibration of a UHF airborne radar," in Proc. IEEE Radar Conference, Arlington, VA, USA, May 2015, pp. 1720{1724. [25] Y. Kim and J. J. van Zyl, \A time-series approach to estimate soil moisture using polarimetric radar data," IEEE Trans. Geosci. Remote Sens., vol. 47, no. 8, pp. 2519{2527, Aug. 2009. [26] S.-B. Kim, L. Tsang, J. T. Johnson, S. Huang, J. J. van Zyl, and E. G. Njoku, \Soil moisture retrieval using time-series radar observations over bare surfaces," IEEE Trans. Geosci. Remote Sens., vol. 50, no. 5, pp. 1853{1863, May 2012. [27] V. E. Romanovsky and T. E. Osterkamp, \Thawing of the active layer on the coastal plain of the Alaskan Arctic," Permafrost and Periglacial Processes, vol. 8, no. 1, pp. 1{22, 1997. [28] T. E. Osterkamp and V. E. Romanovsky, \Freezing of the active layer on the coastal plain of the Alaskan Arctic," Permafrost and Periglacial Processes, vol. 8, no. 1, pp. 23{44, 1997. [29] V. E. Romanovsky and T. E. Osterkamp, \Eects of unfrozen water on heat and mass transport processes in the active layer and permafrost," Permafrost and Periglacial Processes, vol. 11, no. 3, pp. 219{239, 2000. [30] W. L. Cable, V. E. Romanovsky, and M. T. Jorgenson, \Scaling-up permafrost thermal measurements in western Alaska using an ecotype approach," The Cryosphere, vol. 10, no. 5, pp. 2517{2532, 2016. [31] D. Zona, B. Gioli, R. Commane, J. Lindaas, S. C. Wofsy, C. E. Miller, S. J. Dinardo, S. Dengel, C. Sweeney, A. Karion, R. Y.-W. Chang, J. M. Henderson, 109 P. C. Murphy, J. P. Goodrich, V. Moreaux, A. Liljedahl, J. D. Watts, J. S. Kimball, D. A. Lipson, and W. C. Oechel, \Cold season emissions dominate the Arctic tundra methane budget," PNAS, vol. 113, no. 1, pp. 40{45, Jan. 2016. [32] F. T. Ulaby, R. K. Moore, and A. K. Fung, Microwave Remote Sensing: Ac- tive and Passive, Volume II: Radar Remote Sensing and Surface Scattering and Emission Theory. Reading, MA: Addison Wesley, 1982. [33] A. Tabatabaeenejad and M. Moghaddam, \Bistatic scattering from three- dimensional layered rough surfaces," IEEE Trans. Geosci. Remote Sens., vol. 44, no. 8, pp. 2102{2114, Aug. 2006. [34] ||, \Inversion of subsurface properties of layered dielectric structures with random slightly rough interfaces using the method of simulated annealing," IEEE Trans. Geosci. Remote Sens., vol. 47, no. 7, pp. 2035{2046, Jul. 2009. [35] E. Chapin, A. Chou, J. Chen, B. Heavey, S. Hensley, Y. Lou, R. Machuzak, and M. Moghaddam, \AirMOSS: An airborne P-band SAR to measure root-zone soil moisture," in Proc. IEEE Radar Conference, Atlanta, GA, USA, May 2012, pp. 693{698. [36] H. Rott, R. E. Davis, and J. Dozier, \Polarimetric and multifrequency SAR signatures of wet snow," in Proc. IEEE IGARSS, Houston, TX, USA, May 1992, pp. 1658{1660. [37] M. Burgin, D. Clewley, R. M. Lucas, and M. Moghaddam, \A generalized radar backscattering model based on wave theory for multilayer multispecies vegeta- tion," IEEE Trans. Geosci. Remote Sens., vol. 49, no. 12, pp. 4832{4845, Dec. 2011. [38] C. Homer, J. Dewitz, L. Yang, S. Jin, P. Danielson, G. Xian, J. Coulston, N. Herold, J. Wickham, and K. Megown, \Completion of the 2011 National Land Cover Database for the conterminous United States-Representing a decade of land cover change information," Photogrammetric Engineering and Remote Sensing, vol. 81, no. 5, pp. 345{354, 2015. [39] S. L. Durden, J. J. van Zyl, and H. A. Zebker, \Modeling and observation of the radar polarization signature of forested areas," IEEE Trans. Geosci. Remote Sens., vol. 27, no. 3, pp. 290{301, May 1989. [40] R. Engstrom, A. Hope, H. Kwon, D. Stow, and D. Zamolodchikov, \Spatial distribution of near surface soil moisture and its relationship to microtopography in the alaskan arctic coastal plain," Hydrology Research, vol. 36, no. 3, p. 219234, 2005. [41] D. M. Lawrence and A. G. Slater, \Incorporating organic soil into a global climate model," Climate Dynamics, vol. 30, p. 145160, 2008. 110 [42] C. L. Ping, J. D. Jastrow, M. T. Jorgenson, G. J. Michaelson, and Y. L. Shur, \Permafrost soils and carbon cycling," SOIL, vol. 1, pp. 147{171, 2015. [43] C. Kaiser, H. Meyer, C. Biasi, O. Rusalimova, P. Barsukov, , and A. Richter, \Soil moisture retrieval using time-series radar observations over bare surfaces," Journal of Geophysical Research: Biogeosciences, vol. 112, p. G02017, 2007. [44] A. Rodionov, H. Flessa, M. Grabe, O. A. Kazansky, O. Shibistova, , and G. Guggenberger, \Organic carbon and total nitrogen variability in permafrost- aected soils in a forest tundra ecotone," European Journal of Soil Science, vol. 58, no. 6, p. 12601272, 2007. [45] D. H. Boelter, \Physical properties of peats as related to degree of decompo- sition," Soil Science Society of America Journal, vol. 33, no. 4, pp. 606{609, 1969. [46] M. G. Letts, N. T. Roulet, N. T. Comer, M. R. Skarupa, and D. L. Verseghy, \Parametrization of peatland hydraulic properties for the Canadian land surface scheme," Atmosphere-Ocean, vol. 38, no. 1, pp. 141{160, 2000. [47] M. C. Dobson, F. T. Ulaby, M. T. Hallikainen, and M. A. El-Rayes, \Microwave dielectric behavior of wet soil-part II: Dielectric mixing models," IEEE Trans. Geosci. Remote Sens., vol. GE-23, no. 1, pp. 35{46, 1985. [48] R. H. Chen, A. Tabatabaeenejad, and M. Moghaddam, \Retrieval of permafrost active layer properties using time-series P-band radar observations,"IEEETrans. Geosci. Remote Sens., 2019. [49] M. F. Hossain, W. Chen, and Y. Zhang, \Bulk density of mineral and organic soils in the Canadas arctic and sub-arctic," Information Processing in Agriculture, vol. 2, p. 183190, 2015. [50] B. J. Cosby, G. M. Hornberger, R. B. Clapp, and T. R. Ginn, \A statistical exploration of the relationships of soil moisture characteristics to the physical properties of soils," Water Resources Research, vol. 20, no. 6, p. 682690, 1984. [51] National Cooperative Soil Survey, National Cooperative Soil Survey Characteri- zation Database. http://ncsslabdatamart.sc.egov.usda.gov/. [52] G. S. Campbell, \A simple method for determining unsaturated conductivity from moisture retention data," Soil Science, vol. 117, no. 6, pp. 311{314, 1974. [53] C.-H. Park, A. Behrendt, E. LeDrew, and V. Wulfmeyer, \New approach for calculating the eective dielectric constant of the moist soil for microwaves," Remote Sensing, vol. 9, p. 732, 2017. [54] F. T. Ulaby, T. H. Bengal, M. C. Dobson, J. R. East, J. B. Garvin, and D. L. Evans, \Microwave dielectric properties of dry rocks," IEEE Trans. Geosci. Re- mote Sens., vol. 28, no. 3, p. 325336, 1990. 111 [55] F. T. Ulaby and M. A. El-rayes, \Microwave dielectric spectrum of vegetation - Part II: Dual-dispersion model," IEEE Trans. Geosci. Remote Sens., vol. GE-25, no. 5, p. 550557, 1987. [56] C. M atzler, Thermal Microwave Radiation: Applications for Remote Sensing. Stevenage, UK: Institution of Engineering and Technology, 2006. [57] P. A. Rosen, S. Hensley, I. Joughin, F. K.Li, S. Madsen, E. Rodriguez, and R. M. Goldstein, \Synthetic aperture radar interferometry," Proceedings of the IEEE, vol. 88, no. 3, p. 333382, 2000. [58] L. Liu, E. E. Jafarov, K. M. Schaefer, B. M. Jones, H. A. Zebker, C. A. Williams, J. Rogan, and T. Zhang, \InSAR detects increase in surface subsidence caused by an Arctic tundra re," Geophysical Research Letters, vol. 41, no. 11, p. 39063913, 2014. [59] R. M. Goldstein, H. A. Zebker, and C. L. Werner, \Satellite radar interferometry: Two-dimensional phase unwrapping," Radio Science, vol. 23, no. 4, p. 713720, 1988. [60] C. W. Chen and H. A. Zebker, \Phase unwrapping for large sar interferograms: statistical segmentation and generalized network models," IEEE Trans. Geosci. Remote Sens., vol. 40, no. 8, p. 17091719, 2002. [61] B. Rabus, H. Wehn, and M. Nolan, \The importance of soil moisture and soil structure for insar phase and backscatter, as determined by fdtd modeling," IEEE Trans. Geosci. Remote Sens., vol. 48, no. 5, p. 24212429, 2010. [62] S. Zwieback, S. Hensley, and I. Hajnsek, \Assessment of soil moisture eects on l-band radar interferometry," Remote Sensing of Environment, vol. 164, p. 7789, 2015. , 112
Abstract (if available)
Abstract
This dissertation investigates radar remote sensing techniques for monitoring permafrost landscapes and estimating active layer soil properties, particularly the active layer thickness (ALT) and soil moisture profile, in the Arctic-boreal region. Synthetic aperture radar (SAR) polarimetry has been widely used to estimate soil moisture and freeze/thaw state
Linked assets
University of Southern California Dissertations and Theses
Conceptually similar
PDF
Active and passive microwave remote sensing synergy for geophysical parameter estimation
PDF
Electromagnetic scattering models for satellite remote sensing of soil moisture using reflectometry from microwave signals of opportunity
PDF
Physics-based bistatic radar scattering model for vegetated terrains in support of soil moisture retrieval from signals of opportunity
PDF
Physics-based models and microwave sensors for mapping soil carbon in arctic permafrost using long-wavelength radars
PDF
Solution of inverse scattering problems via hybrid global and local optimization
PDF
Ultra-wideband multistatic and MIMO software defined radar sensor networks
PDF
FPGA based L-band pulse Doppler radar design and implementation
PDF
Signal processing for channel sounding: parameter estimation and calibration
PDF
Investigating the role of climate in affecting residential electricity consumption through high spatiotemporal resolution observations
PDF
Design and analysis of reduced complexity transceivers for massive MIMO and UWB systems
Asset Metadata
Creator
Chen, Richard H.
(author)
Core Title
Radar remote sensing of permafrost landscapes and active layer soil properties
School
Viterbi School of Engineering
Degree
Doctor of Philosophy
Degree Program
Electrical Engineering
Publication Date
11/07/2019
Defense Date
07/18/2019
Publisher
University of Southern California
(original),
University of Southern California. Libraries
(digital)
Tag
active layer thickness,Earth Sciences,estimation,OAI-PMH Harvest,permafrost,Radar,remote sensing,soil moisture
Language
English
Contributor
Electronically uploaded by the author
(provenance)
Advisor
Moghaddam, Mahta (
committee chair
), Jenkins, Keith (
committee member
), Molisch, Andreas (
committee member
), Sanders, Kelly (
committee member
), Tabatabaeenejad, Alireza (
committee member
)
Creator Email
chenrh@usc.edu,richard.oy.chen@gmail.com
Permanent Link (DOI)
https://doi.org/10.25549/usctheses-c89-230997
Unique identifier
UC11674873
Identifier
etd-ChenRichar-7901.pdf (filename),usctheses-c89-230997 (legacy record id)
Legacy Identifier
etd-ChenRichar-7901.pdf
Dmrecord
230997
Document Type
Dissertation
Rights
Chen, Richard H.
Type
texts
Source
University of Southern California
(contributing entity),
University of Southern California Dissertations and Theses
(collection)
Access Conditions
The author retains rights to his/her dissertation, thesis or other graduate work according to U.S. copyright law. Electronic access is being provided by the USC Libraries in agreement with the a...
Repository Name
University of Southern California Digital Library
Repository Location
USC Digital Library, University of Southern California, University Park Campus MC 2810, 3434 South Grand Avenue, 2nd Floor, Los Angeles, California 90089-2810, USA
Tags
active layer thickness
estimation
permafrost
remote sensing
soil moisture