Mapping Freezing and Thawing Surface State Periods With the CYGNSS Based F/T Seasonal Threshold Algorithm

Freeze/Thaw (F/T) surface state detection and imaging has been demonstrated as a capability of NASA's Cyclone Global Navigation Satellite System mission. The objective of the research presented here is the application of an existing F/T retrieval algorithm to multiple target areas over several years in order to demonstrate the ability to provide long term trending of F/T behavior. Metrics are developed and evaluated to characterize F/T trends related to global warming. They include the annual fraction of time spent frozen and thawed as a function of location and year, and the earliest fall freeze and earliest spring thaw as a function of location and year. The results demonstrate the capability of Global Navigation Satellite Systems Reflectometry to provide long term trending of F/T behavior relevant to studies of climate change.


I. INTRODUCTION
C LIMATE change is affecting many socio-economic activities, including agriculture production because of the warmer temperatures and the reduced water availability [1], [2], [3]. Climate change action requires the accurate determination of many parameters that characterize atmosphere, ocean, and land processes, their interaction; as well as the ability to predict future changes. Freeze/Thaw (F/T) surface state detection and monitoring is important for quantifying carbon, energy, and water fluxes, and their impact on land cover change. At present, there is a lack of operational space-borne measurements of active layer thickness and extension. In this study, the Cyclone Global Navigation Satellite System (CYGNSS) mission Global Navigation Satellite Systems Reflectometry (GNSS-R) data [4], [5], [6], [7], [8], [9], [10], [11], [12], [13] are used to characterize F/T surface state changes due to water-ice phase conversion. These changes are an indirect measurement that can be used to characterize the active layer. The main objective is to demonstrate the potential of GNSS-R based F/T surface state monitoring by future high-inclination orbit GNSS-R missions.
CYGNSS was selected by NASA as a low-cost and highscience Earth Venture Mission in 2012. In 2016, it was successfully launched into space [14], [15]. The original mission Manuscript  objective was wind speed estimation over tropical cyclones (TCs). Traditional remote sensing satellites are not able to access TCs inner core when there are strong precipitations i.e., heavy rain, while CYGNSS does. The CYGNSS orbital configuration was designed to accomplish this goal. Consequently, CYGNSS operates from a low Earth orbit with an inclination angle of ∼35 • . The mission coverage is within the latitude band ∼ [−37. 5 37.5] • . At present, the science team has significantly expanded the number of scientific applications over land and atmosphere.
CYGNSS offers an unprecedented spatio-temporal sampling of the Earth surface, helping to further understand F/T surface state dynamics. Consequently, in 2020, we initiated a research project to explode this strong potential for F/T studies. A seasonal-threshold algorithm (STA) [19] was developed and validated using surface temperature data as provided by the European Centre for Medium-Range Weather Forecast (ECMWF) ERA5-Land numerical reanalysis model. The agreement between CYGNSS and ERA5-Land F/T maps was found to be quite good in winter and summer (see Fig. 10 and Table I in [19]). During the transitional spring and fall months, the F/T missed detection increases, possibly as a result of uncertainties in the reference models during these periods (see Fig. 10 and Table I in [19]).
In addition, CYGNSS-derived F/T surface state was evaluated and an inter-comparison with the SMAP-radiometer F/T data product was performed (see Fig. 11 in [19]). The SMAP F/T product was derived from the normalized polarization ratio of radiometer measurements. A good agreement was found with both CYGNSS and ERA5-Land maps. Finally, a time-series analysis was performed (see Fig. 12 in [19]) over two key areas with high and low altitude terrain. The scale factor Δ(t) was selected, which is the fundamental input observable used by the retrievals. Over the high-altitude target area, a significant increment was found from fall Δ(t) ∼0.2 to winter Δ(t) ∼1.2. Over the low-altitude terrain, Δ(t) remained low throughout the year. This manuscript corresponds to the second part of our research project. Here, we apply our previously validated F/T CYGNSSbased retrieval algorithm [19] over different target areas within the CYGNSS coverage (the Tibet Plateau, the Andes Mountains, and the Rocky Mountains). The ability to detect and monitor long term trends related to global warming are demonstrated by examining the year-to-year changes in the following F/T related parameters: 1) The yearly number of frozen and thawed months at each pixel.
2) The first month after the warm season in which each pixel is frozen (earliest fall freeze).
3) The first month after the cold season in which each pixel is thawed (earliest spring thaw). First, we introduce the selected datasets in Section II. Then, the different target areas are described in detail in Section III. A summary of the STA and the specific considerations for this study are presented in Section IV. The results are presented in Section V, including the interpretation of the maps over the freezing and thawing periods. Final discussions are included in Section VI. Finally, Section VII concludes this article.

A. CYGNSS GNSS-R Data: Properties
CYGNSS Earth's surface sampling properties are determined by the GNSS-R multibistatic radar geometry, which is encouraging for the purpose of remote sensing. In the CYGNSS scenario, this geometry is defined by 32 Global Positioning System (GPS) transmitters and 8 GNSS-R receivers, with an orbital height of ∼20 000 km and ∼520 km, respectively. CYGNSS is an 8microsatellites single-plane GNSS-R constellation. Each single on-board GNSS-R instrument is connected to two down-looking antennas, which point to opposite sides of the ground track with an incidence angle of ∼28°. Each of these receivers is configured to collect scattered GPS signals along four specular directions simultaneously. As such, CYGNSS enables to study the Earth system along 32 tracks simultaneously with a wide range of satellite's incidence angles ∼[0 70]°.
The CYGNSS v3.0 product is selected for this study. This version of the CYGNSS family products incorporates dynamic monitoring of transmitted GPS power [22]. This strategy is required to afford the calibration challenge given by the variable power mode of the GPS Block IIF and IIR-M satellites. As such, this product represents an important improvement as compared to the CYGNSS v2.1 version. The selected temporal window for this study is 31 months, which includes two complete full F/T cycles along 2019 and 2020, and the first 7 months of 2021. This window is long-enough so as to enable the evaluation of the capability of our retrieval algorithm in determining temporal changes of the F/T surface state.

B. Topographic Features: Elevation and Roughness
Given the CYGNSS latitudinal coverage, this F/T study is performed over high-elevation areas because of the higher probability to find F/T transitions as compared to low-elevation terrain. The reference surface elevation is derived from the Global Multiresolution Terrain Elevation Data (GMTED 2010)  In this study, the 250 m version is selected. Surface elevation and topographic roughness index (TRI) are derived from this product, and results are up-scaled to 5 km spatial resolution.
While surface elevation (see Fig. 1) is an important parameter for the selection of the target areas, TRI [24] is also relevant for the analysis over these areas. TRI is zero over flat surface, while positive TRI values appear over mountainous areas with steep slopes. TRI is estimated as the mean of the absolute differences of the surface elevation of a center-cell and the neighboring 8 cells as follows: where H is the surface elevation given by the GMTED 2010 DEM, N = (n − 1)/2, and n is any odd integer smaller than the number of cells in the shortest edge of the raster.

C. ECMWF ERA5-Land Surface Temperature
In this study, the reference surface temperature in the first layer of the soil (0-7 cm) is given by the ECMWF ERA5-Land [see Fig. 2(a) and (b)] [25]. The first layer of the soil is selected because GNSS L-band signals can penetrate the top ∼5 cm of the soil. ERA5-Land is a reanalysis dataset that enables a reliable understanding of several land variables from the last 7 decades. The land surface model used in ERA5-Land is the Carbon Hydrology-Tiled ECMWF Scheme for Surface changes over Land (CHTESSEL) [26]. ERA5-Land corresponds to the improved version of the ERA5 products. It uses atmospheric variables (air temperature and air humidity) to control the modelbased estimates, but these data are not directly assimilated to generate the products. ERA5-Land products have been successfully cross-compared with a significant number of in situ measurements (2001-2018), models, and space-borne reference datasets [27].

A. Tibet Plateau
The Tibet Plateau is also known as the "Third Pole" of the Earth. It is the largest high-altitude (∼4000 m) area on the planet [28], [29]. It has a great impact in water storage and supply in Asia. Both permafrost and seasonal frozen ground can be found over the Tibet Plateau. The dominant International Geosphere Biosphere Program (IGBP) land cover types are grassland, barren, and shrublands. The TRI is high (>70) over the Himalayan Mountains and low (<20) over the central part of the Tibet Plateau [see There is a strong temperature gradient between both seasons, with extended permafrost areas over high-altitude terrain.
During the last decades, the extent of permafrost is being reduced due to global warming induced by climate change [30]. On the other hand, recent studies indicate that, because of global warming, shrublands, open forests, grasslands, and water bodies are growing in the area, while the extension of evergreen and deciduous broadleaf forests, and barrens is reducing [31].

B. Andes Mountains
The length of the Andes Mountains is ∼8000 km along the west coast of South America. This area is quite heterogeneous in terms of topography and climate. It is mostly covered by   ]. There is a strong surface temperature gradient between both seasons, without permafrost areas over the temporal window of this study i.e., from January 2019 to July 2021. A major part of studies for the detection of potential permafrost over this region were based on models because of the difficult access to these remote and high-altitude mountains. In 2012, a non-negligible permafrost probability in the latitude range Lat = [−36 −23]°was reported [33].
The impact of climate change is expected to be stronger on the arid Andes, significantly reducing the amount of available water for millions of people [34], [35]. This is a clear example of how climate change will have negative socio-economic impacts on some countries. On the other hand, over different regions such as the Tibet Plateau (China), Canada, and Siberia (Russia), climate change could transform cold deserts into extended areas suitable from agriculture, which is currently one of the major needs, e.g., in China.

C. Rocky Mountains
The length of the Rocky Mountains is ∼4800 km from western Canada to the southwestern United States of America (USA). There is a wide range of climates in the Rocky Mountains, from semiarid shortgrass prairies to alpine tundra. The dominant IGBP land cover types are grassland, evergreen broadleaf forest, and deciduous needleleaf forest. The TRI is from high-tomoderate (>50) over the Rocky Mountains and low (<20) over the prairies [see Fig. 5(a)]. The ECMWF ERA5-Land surface temperatures are also shown in winter [see Fig. 5(b)] and summer [see Fig. 5(c)]. There is a strong temperature gradient between both seasons.
Climate change is reducing snowpack levels in the Rocky Mountains, which is the main water storage for the central part of USA, including highly populated areas such as Colorado, with its important industrial/technological hub.

A. Fundamental Observable
The fundamental observable used in this study is the CYGNSS v3.0 power DDM after radiometric calibration. Power DDMs are the true power as measured by an ideal analog sensor without additional postprocessing. This selection enables a more accurate and focused estimation of the surface reflectivity than using data products such as the signal-to-noise ratio or processed observables such as the normalized bistatic radar cross section, which depend on external algorithms and scattering assumptions.

B. Metadata and Quality Flags
CYGNSS products are dynamically improving in performance and are increasing the number of available metadata and quality flags. In this study, the following metadata are used: 1) transmitter equivalent isotropically radiated power (EIRP); 2) transmitter and receiver antenna gains (G t and G r ); 3) distance between each CYGNSS spacecraft and the nominal specular points (D t ); 4) distance between each GPS satellite and the nominal specular points (D r ); 5) geographical coordinates of the nominal specular point.
The equivalent overall quality flag over land surfaces is applied to select the highest-quality dataset: large attitude control error, black body DDM, GNSS-R instrument reconfigured, invalid cyclic redundancy check, test antenna pattern DDM, channel idle, confidence in DDM noise floor is low, big gap in noise floor, big gap in low noise amplifier temperature, reflected DDM containing direct signal, low confidence in EIRP estimation, radio frequency interference is detected, bistatic radar cross-section (BRCS) DDM specular point bin delay error, BRCS DDM specular point bin Doppler error, GPS positioning velocity time SP3 propagator error, specular point nonexistent error, BRCS look up table (LUT), antenna data LUT range error, baseband framing error, flight software compensation shift error, spacecraft attitude out of nominal mode, sampling period error, invalid roll state, GNSS-R instrument not valid antenna selection, specular point is in antenna sidelobe, and low gain in zenith antenna.

C. Processed Observable
Final DDMs are obtained after the interpolation of the improved 1700 × 1100 bins DDMs. In so doing, a spline method is used. This improves the quality of the original 17×11 bins fundamental DDMs. This improvement is implemented to minimize the potential impact of topography in the tracking of the peak of the DDMs, which is used for the estimation of the surface reflectivity.

D. Reflectivity Estimation
The reflectivity is estimated as it follows: where λ is the signal wavelength, P is the peak value of the power DDMs, and N is the noise floor, which are calculated from the power DDMs. Several previous F/T studies [16], [17], [18], [19] have applied this reflectivity estimation, which is based on the Friis transmission formula [36], [37]. Pioneering investigations discussed the influence of topography and up-welling biomass cover [38], [39], [40], [41]. The main parameter is the vegetation optical depth (VOD), which causes the attenuation of the reflected signal. Over the three selected case-studies in this study, VOD levels are low (<0.2) and can be neglected, especially in the winter season. Over the highest-interest target area in this study, i.e., the central part of the Tibet Plateau [see Fig. 3(a)], the topographic roughness is low (<20). All of these considerations triggered our decision to use the reflectivity estimation method shown in (2). There is a relationship between the reflectivity and the Fresnel reflection coefficient, which depends on the soil permittivity * s . When the soil surface changes from frozen to thawed, the permittivity increases. The ability of CYGNSS to measure surface reflectivity is high because of the strong sensitivity of L-band GNSS signals to the permittivity, which depends on the phase of the water, and because of the higher biomass penetration depth as compared to higher frequencies, starting at C-band.

E. Surface Gridding Strategy
CYGNSS surface sampling properties are pseudorandom. Surface gridding and spatio-temporal averaging strategies are needed to generate continuous land surface spatial information [42], [43]. The size of the grid is a tradeoff between the required spatial resolution for the study of a specific surface variable (soil moisture, biomass, inland waters) and the temporal resolution associated with the dynamics of the specific variable. The temporal resolution depends on the CYGNSS revisit time.
The size of the grid should account for the available number of CYGNSS measurements within the required time interval. Additionally, a smoothing approach is needed to account for the pseudorandom nature of the surface sampling and the impact of the incidence angle on the spatial resolution of each individual CYGNSS measurement. This is especially important in areas with rough topography, so as to minimize the impact of individual multiangle topographic samples by averaging across the range of angles seen during the averaging time.
In this study, a 0.1°× 0.1°grid is used and a moving-averaging filter of 0.1°× 0.1°is applied at steps of 0.1°. The rationale behind this selection is to increase the number of measurements per pixel so as to generate a statistically significant population. Additionally, regarding this respect, the temporal averaging is set to 1 month to enable the generation of the required number of measurements, and because this is a long term study focused on the yearly evolution of the F/T cycle. Because of the pseudorandom sampling properties, the grid size is larger than in [19] to reduce the differences between the three different target areas.

F. Seasonal Threshold Algorithm
The retrievals are based on our previously validated F/T STA [19]. This algorithm evaluates the relationship between time series of CYGNSS-derived reflectivity Γ and seasonal reference frozen and thawed states. For a measurement at time t, the seasonal scale factor Δ(t) is defined as follows [19]: where Γ(t) is the reflectivity measurement estimated at time t, and Γ fr and Γ th are reflectivity measurements corresponding to frozen and thawed reference states, respectively. The frozen and thawed references states are calculated as the 10% and the 90% points in the cumulative distribution function of the reflectivity at each pixel, during winter (Γ fr ) and summer (Γ th ). Different F/T surface states correspond to different observation times. Two threshold levels T fr and T th are used for the F/T classification. The surface state is frozen when Δ(t) > T fr . The surface state is thawed when Δ(t) < T th . Both threshold levels T fr and T th are selected at the optimum operating points of the receiver operating characteristic (ROC) freeze and thaw curves, using ERA5-Land surface temperatures as proxy data. In this study, the optimum operating points of the ROC freeze and thaw curves are defined at the points of inflection when the slope of the curve crosses from above to below unity. This point marks the transition from more true than false detections to more false than true.
The performance of the STA is directly related to the estimation of the surface reflectivity Γ. This justifies the selection of the fundamental observable and the methodology used to generate the improved observable. The reflectivity Γ depends on the Fresnel reflection coefficient, and thus, on the surface permittivity * s . As such, the scale factor Δ(t) is directly related to transitions from frozen to thawed soil. Δ(t) is computed at each pixel over the target areas. Thus, the STA accounts for the vegetation and land cover type heterogeneity with a spatial resolution given by the grid size. Since Δ(t) is a relative factor based on differential measurements, the impact of the per pixel land cover type is mitigated. Consequently, the two threshold levels T fr and T th can be used independently of the characteristics of each pixel.

V. RESULTS
In this section, the results over the selected target areas, i.e., the Tibet Plateau (see Fig. 3), the Andes Mountains (see Fig. 4), and the Rocky Mountains (see Fig. 5), are presented. Fig. 6 maps the total number of months in 2019 for which the surface was frozen and thawed in each of the three target areas. The assessment is done pixel-by-pixel in order to allow the effects of geographic features to be resolved. Seasonal frozen ground is located within the Tibet Plateau, while permanent thawed soil is mainly located at the south of the Himalayas. Additional maps for the year 2020 have also been generated and difference between from year to year can used to examine longer term climatic trends.
For each of the three target areas, Fig. 7 maps the first month of the year in which each pixel is thawed, which marks the beginning of the spring season, together with the first month after summer in which each pixel is frozen, marking the beginning of the winter season. The freezing and thawing periods strongly determine the flow of soil water over the surface. This movement has an impact on vegetation growth, as well as on changes in land cover type, which in turn has an influence on the frequency of the F/T cycle [44]. Monitoring these changes could influence agriculture development plans, which is a fundamental economic concern.
In the case of the Tibet Plateau, it is found that the beginning of the winter season appears sooner over areas with a higher altitude (see Fig. 1) than in the central part of the Plateau. On the other hand, the beginning of the spring season appears later over the Himalayas than over the central part of the Tibet Plateau, because of the decreasing temperatures with increasing altitude.
The thawing of soils during the thawing period triggers the delivery of organic carbon and acceleration of soil respiration. These results indirectly suggest that CO2 and CH4 emissions are dynamically being geo-distributed [45], [46]. Future long-term sustainable GNSS-R datasets could enable to further understand the effects of climate change on a global scale [47].
A similar methodology for the analysis is performed over the Andes and the Rocky Mountains. Seasonal frozen regions are found in these target areas because of the high temperatures in summer [see Fig. 4(b) and Fig. 5(c)], which prevent the existence of permafrost. Both target areas have a lower altitude than the Tibet Plateau and a smaller surface extension, but the geographic variability of the altitude is stronger. The central part of the Tibet Plateau has a roughly constant altitude of ∼4000 m, while the highest altitudes are clustered over the Himalayas ∼8000 m.
Over the Andes Mountains, the number of frozen months is relatively low ∼3-4. It is worth commenting that the complete Andes region within the CYGNSS coverage was evaluated, but F/T cycling was only found in the selected target area Lat = [−37 .5 −25]°. It appears that in this target area, the F/T cycling also depends on the latitude, being higher at lower latitudes because of the longer duration of the cold season, although the highest elevations are towards to northern part of the target area (see Fig. 1). These results, and the overall methodology presented in this study could help to complement and revisit the few model-based studies [32] that have been performed in this region.
Over the Rocky Mountains, the number of frozen months is also relatively low ∼3-4. The surface elevation is a parameter with an important influence in the number of frozen months. The freezing period indicates that the first month when a frozen pixel can be found is closer to the end of the year than over the Tibet Plateau and the Andes Mountains.
Finally, Fig. 8 is included for a comparison among the three target areas. This figure shows the temporal evolution of the extension of the total frozen area over each target region. A repeatable seasonal variability is found with periodic changes from summer to winter. Over the Tibet Plateau, the extension of the frozen area does not decrease down to 0 km 2 , which indicates the existence of permafrost [see Fig. 8 Additionally, a clear trend is found that shows how the extension of the frozen area is gradually reducing from 2019 to 2021, except for the peak in 2021 over the Rocky Mountains. This is due to a quite cold winter season in the region [48], and the greater susceptibility to annual fluctuations because of the lower extension as compared to the Tibet Plateau and the Andes. Although the moderate length of the time-period under study 2019-2021 prevents making conclusive remarks about the impact of climate change, these results demonstrate the potential of our F/T retrieval algorithm for future long term studies using GNSS-R.

VI. FINAL DISCUSSIONS
Passive microwave multichannel and high-frequency sensors such as, e.g., Special Sensor Microwave/Imager (SSM/I) [49], Advanced Microwave Scanning Radiometer-Earth (AMSR-E) [50], and AMSR-2 [51] have been used to study F/T surface state. The fundamental operating principle is the different emission properties observed by the different channels. L-band microwave missions such as Soil Moisture Ocean Salinity (SMOS) [52] and SMAP [53] have also been considered for F/T mapping because of the enhanced sensitivity to the soil surface permittivity at L-band. Additionally, synthetic aperture radar (SAR) missions such as, e.g., C-band ERS-1 [54], L-band Phase Array type L-band SAR (PALSAR) [55], dual-pol C-band Sentinel-1 [56], and airborne P-band [57], as well as C-band Meteorological Operational (MetOp) Advanced Scatterometer (ASCAT) [58] have been operated for F/T detection using backscattered signals properties.
The spatial and temporal resolution of SAR observational systems is ∼100 m and ∼7-14 days, respectively. On the other hand, the spatial and the temporal resolution of microwave radiometers is ∼25 km and ∼3 days, respectively. Both types of sensors are well suited remote sensing options for F/T studies. However, at present, there are still limitations generating F/T products with adequate spatio-temporal sampling and spatial resolution simultaneously to understand the dynamics of F/T transition periods. GNSS-R offers a promising addition to the existing estimates of the F/T state because of the expected significant increase in number of missions in the coming years, including the current constellation of private satellites operated by Spire [59] and the upcoming Hydrology using GNSS reflections (HydroGNSS) mission by European Space Agency [60]. Future activities could include the study of F/T daily-fluctuations to explore the dynamics of the F/T cycle. GNSS-R, SAR, and microwave radiometry mission F/T products could be used in a synergistic manner to enhance these studies.
The expected impact of this scientific-technological development is high because of the influence on climate change actions and socio-economic aspects, including crop production and geopolitics. China has ∼8% of the Earth' land surface suitable for agriculture, but the ∼25% of the world's population. USA has ∼15% of the Earth' land surface suitable for agriculture, but the ∼4% of the world's population. China is the major customer of USA agriculture exports. The ∼60% of China's trade is by sea. Thus, the economic security of China strongly depends on the control of South China Sea. In the near future, warmer temperatures will significantly increase the available land extension suitable for agriculture in the Tibet Plateau. This could have an important influence on the world's international relationships.

VII. CONCLUSION
Soil F/T transition monitoring is essential for quantifying climate change and hydrologic dynamics over cold regions. F/T cycling is a key-parameter for the evolution of ecosystems. The F/T changes of seasonally frozen ground are an important indicator of climate change. This article is performed using 31 months of CYGNSS data over three different high-altitude target areas: the Tibet Plateau, the Andes, and the Rocky Mountains. These regions are selected based on global surface elevation data from the GMTED 2010 DEM. TRI and ERA5-Land surface temperature are also used to improve the interpretation of the results. Our F/T STA is applied to evaluate the yearly number of frozen and thawed months at each pixel, the first month after the warm season (earliest fall freeze) in which each pixel is frozen as well as the first month after the cold season (earliest spring thaw) in which each pixel is thawed. Results show the capability to capture interannual changes in these parameters and to generate maps of their spatial variation. Extensive permafrost areas are found over the Tibet Plateau. Over the Andes and the Rocky Mountains, only seasonal frozen ground is detected. From 2019 to 2021, the loss of the maximum extension of frozen area in the Tibet Plateau and the Andes Mountains is ∼100 000 km 2 and ∼15 000 km 2 , respectively.
We demonstrate the potential of GNSS-R to further understand the geo-distribution and long term temporal variability of CO2 and CH4 at a global scale. The higher spatio-temporal sampling of GNSS-R as compared to more traditional remote sensing techniques could open new insights in monitoring highly dynamic F/T surfaces processes. Future high-inclination orbit GNSS-R missions, e.g., [59], [60] will enable the application of our retrieval algorithm over higher latitude regions. Within the new space era, the number of SmallSats and constellations of SmallSats will significantly increase the number of data sources as well as both the surface temporal sampling and the latency time. This will enable real climate change action, and it will increase the relevance of GNSS-R and remote sensing in future geopolitics.