A high-resolution (0.05°) global seamless continuity record (2002–2023) of near-surface soil freeze-thaw states via passive microwave and optical satellite data

The global near-surface soil freeze-thaw (FT) states are crucial for understanding complex interactions with hydrological, ecological, and climatic processes. However, current remote sensing of FT states primarily relies on passive microwave remote sensing, which, despite its all-weather monitoring capabilities, suffers from low spatial resolution. This limitation restricts its application to hydroclimatological scales, precluding its use in finer-scale studies such as soil erosion and hydrometeorological applications. To address this issue, this study introduces a novel downscaling approach that integrates passive microwave and optical satellite data to generate a long-term (2002–2023), high-resolution (0.05°) dataset of global near-surface FT states, ensuring daily seamless continuity. The dataset was validated against in situ measurements, demonstrating that the high-resolution product maintains an overall accuracy of 83.78 %, consistent with the coarse-resolution microwave-based dataset, while offering enhanced spatial detail. Comprehensive global trend analyses provided new insights into the dynamics of FT cycles, revealing that the average annual number of frost days in regions north of 45° N is 187.8 ± 12.7 d, with 14.35 % of the area showing a decreasing trend in frozen persistence. Additionally, the average annual number of freeze onset dates is 240.3 ± 7.2, and 9.10 % of the area exhibits a trend of delayed freeze onset. The high-resolution record enables accurately monitoring FT states and providing detailed information for a refined understanding of hydrological and ecological effects globally. The global 0.05° near-surface soil FT state dataset (FT-HiDFA) is freely available at https://doi.org/10.11888/Cryos.tpdc.301551 (Zhao et al., 2024b).

1 Introduction

Frozen ground, including permafrost and seasonally frozen ground, covers nearly 66×10⁶ km², accounting for 52.5 % of the total global land surface area (Kim et al., 2011; McDonald and Kimball, 2006). Most of these regions experience seasonal freeze-thaw (FT) state cycles, which refer to the phase transition between water and ice in the soil pores of surface layers (Zhang et al., 2010). These FT state cycles significantly impact surface runoff, energy balance, and carbon cycling (Kimball et al., 2004; Wang et al., 2024), thereby influencing climate (Peng et al., 2016; Poutou et al., 2004), hydrological (Gouttevin et al., 2012; Gray et al., 1985), ecological (Black et al., 2000), and biogeochemical (Panneer Selvam et al., 2016; Schaefer et al., 2011; Xu et al., 2013) processes. Therefore, the dynamics of surface soil FT states are recognized as a significant indicator of global climate warming, with broad implications for the Earth's environment and human society (Beer et al., 2018; Yang et al., 2013).

Over the past few decades, numerous studies have focused on detecting soil FT status. Although many traditional methods based on in situ observations offer advantages in terms of accuracy and reliability (Wei et al., 2011; Cary et al., 1979; Zhang et al., 2007), their labour-intensive nature, coupled with sparse spatial coverage and limited representativeness, restricts their ability to meet the comprehensive requirements of frozen ground research. In contrast, satellite remote sensing offers an effective and rapid alternative, providing extensive spatial and temporal coverage for monitoring FT status.

Due to the critical importance of monitoring global surface FT states in hydroclimatology research, FT has been identified as an essential climate variable (ECV) by the Global Climate Observation System (GCOS) (Bojinski et al., 2014). FT cycles occur frequently and are highly spatially heterogeneous due to variations in topography, vegetation, soil properties, and snow characteristics (Chang et al., 2015; Jiang et al., 2020). As a result, the GCOS has specified explicit requirements for FT estimation. Although an ideal spatial resolution of 1 km and a temporal resolution of 6 h are recommended, the application of FT monitoring can be improved if FT records achieve a spatial resolution of 10 km and daily temporal resolution. However, current remote sensing techniques are unable to directly meet these requirements for monitoring FT dynamics.

Satellite microwave remote sensing methods are less susceptible to potential degradation from solar illumination effects and atmospheric cloud/aerosol contamination while remaining sensitive to changes in landscape dielectric properties due to liquid water content variations in soils. This characteristic makes them highly applicable in consistent FT detection (England, 1990; Kou et al., 2017; McDonald et al., 2004b; Wu et al., 2022; Zhao et al., 2014). Passive microwave remote sensing is an effective technique for monitoring global surface FT processes owing to radiometers' shorter revisit intervals and larger coverage areas (Kim et al., 2011; Kou et al., 2017; Zuerndorfer and England, 1992). Additionally, it can detect microwave radiation from specific depths within the surface soil and is sensitive to changes in soil dielectric properties, thereby enhancing its applicability in FT detection (McDonald et al., 2004b). In recent years, with the development of new-generation passive microwave radiometers, various corresponding FT discrimination algorithms have been constructed, including the dual-index algorithm (Han et al., 2015; Judge et al., 1996; Zuerndorfer and England, 1992; Zuerndorfer et al., 1990), the decision tree algorithm (Jin et al., 2009), the discriminant function algorithm (DFA) (Hu et al., 2019; Kou et al., 2017, 2018; Zhao et al., 2011), the seasonal threshold method (Kim et al., 2011), and the polarization ratio (PR)-based algorithm (Rautiainen et al., 2016; Roy et al., 2015). However, the relatively coarse spatial resolution of current passive microwave radiometers limits the retrieval of high-resolution information (Chai et al., 2014; Han et al., 2015; Zhao et al., 2011; Zhou et al., 2016). Moreover, spatial heterogeneity, often caused by mixed pixels resulting from the coarse spatial resolution, introduces uncertainty into FT data derived from passive microwave remote sensing.

Optical remote sensing techniques provide high-spatial-resolution information, such as land surface temperature (LST), which can be used to infer surface FT states. For example, the Moderate Resolution Imaging Spectroradiometer (MODIS) on the Aqua and Terra satellites provides LST products with high accuracy, potentially allowing for monitoring global FT states. However, LST primarily reflects the linear changes in temperature during the FT process and cannot capture the abrupt changes in dielectric properties that occur between the frozen and thawed soils. In addition, these products are significantly affected by discontinuities in coverage due to cloud contamination, vegetation, and snow cover (Cary et al., 1979; Langer et al., 2013; Chen et al., 2021; Running, 1998). Recent research focus is exemplified by efforts to generate seamless datasets (Li et al., 2018; Yao et al., 2023; Yu et al., 2022; Zhang et al., 2022, 2020; Zhao et al., 2024a).

Active microwave sensors, divided into radars and scatterometers, receive backscattering coefficients as the echo signals, which are primarily related to soil structure and dielectric properties. Active microwave remote sensing, particularly through synthetic aperture radars (SARs), provides observations of landscape FT states at resolutions on the kilometre scale or finer. Consequently, many active microwave FT discrimination algorithms have been developed, including the FT threshold value method (Du et al., 2015; Kim et al., 2012; Kimball et al., 2004; Way et al., 1997), a change detection algorithm (Frolking et al., 1999), and an edge detection algorithm (Canny, 1986; McDonald et al., 2004a). However, most satellite-based SARs have longer revisit periods than passive microwave radiometers, preventing them from meeting the temporal resolution requirements of FT monitoring.

Existing methods fall short in fulfilling the scientific requirements for high-resolution detection of FT states and are predominantly dependent on microwave observations. To overcome these limitations, downscaling techniques based on multi-source data fusion have been proposed (Giorgi and Mearns, 1991; Hanssen-Bauer et al., 2005; Xu et al., 2019). Among these, statistical downscaling methods (Bierkens et al., 2000; Vaittinada Ayar et al., 2016) are commonly applied in high-resolution monitoring studies (Fan et al., 2005; Hertig and Jacobeit, 2008), which quickly establish statistical or empirical relationships between downscaling predictors and predictands, thereby enabling the reconstruction of coarse-resolution data at finer scales. This has inspired researchers to integrate multi-source remote sensing data, leveraging the high sensitivity of passive microwave signals to FT dynamics and the finer spatial details of optical datasets (Zhao et al., 2017).

The detection of frozen and thawed soil can rely on two critical characteristics: temperatures below 0 °C and the presence of ice. Therefore, remote sensing observations can be utilized to infer these conditions by monitoring the physical temperature and liquid water content. Microwave and optical remote sensing each offer distinct advantages in this context. Microwave observations are extensively employed to detect FT dynamics through various discrimination algorithms. In particular, the DFA has proven to be highly effective in generating consistent, long-term, daily global FT state products. This algorithm, which utilizes brightness temperature (TB) data from the Advanced Microwave Scanning Radiometer for EOS (AMSR-E) and its successor, the Advanced Microwave Scanning Radiometer 2 (AMSR2), has demonstrated both high validation accuracy and reasonable consistency (Hu et al., 2019; Wang et al., 2019a). One of the key discriminant indicators in this algorithm is TB at 36.5 GHz due to its strong correlation with LST. The other important indicator is the quasi-emissivity (Qe), i.e. the ratio between TB at 18.7 GHz and TB at 36.5 GHz, which captures the variations in landscape dielectric properties (ice $/$ water content) across different FT conditions (Zhang et al., 2010; Zhao et al., 2011). Despite its strengths in physical mechanisms, the coarse spatial resolution of passive microwave sensors constrains their application in detailed monitoring. In comparison, optical remote sensing provides data with high spatial resolution, offering an alternative for detecting detailed FT states. A previous study attempted to generate high-resolution FT maps, but the optical data used were restricted to LST and did not incorporate parameters characterizing soil ice $/$ water content (Hu et al., 2017). The apparent thermal inertia (ATI), derived from optical data, has been shown to effectively monitor soil moisture conditions (Qin et al., 2013; Song and Jia, 2016; Van Doninck et al., 2011; Veroustraete et al., 2012; Verstraeten et al., 2006). Thus, integrating ATI into data fusion might further enhance the ability to generate high-resolution FT records (Liang et al., 2024; Yao et al., 2023; Yu et al., 2022; Zhang et al., 2022).

The objective of this study is to enhance the spatial resolution of FT detection products without compromising the accuracy of passive microwave-based FT products, as derived from AMSR-E/2 TB data through the DFA. To achieve this enhancement, the study employs downscaling indicators, specifically the MODIS-based LST and ATI, which serve to encapsulate soil moisture information. Both the original coarse-resolution and the resultant high-resolution FT records are validated by in situ soil temperatures, thereby facilitating an assessment of accuracy variations post-downscaling. Subsequent trend analyses of the high-resolution FT records are conducted to reflect the detailed dynamics of FT states. The resultant downscaled, high-resolution FT states align with the expanding spatial and temporal resolution requirements of GCOS for FT monitoring, thereby providing a valuable tool in cryospheric and ecological studies.

2 Data

2.1 AMSR-E and AMSR2 brightness temperature

The near-surface FT downscaling method integrates data from passive microwave and optical remote sensing. TB observations from passive microwave sensors were obtained from AMSR-E and its successor, AMSR2. AMSR-E, which was carried on NASA's Aqua satellite, operated from May 2002 to October 2011 and provided six microwave bands (6.9, 10.65, 18.7, 23.8, 36.5, and 89 GHz), each available in both horizontal (H) and vertical (V) polarization (Cho et al., 2017; Knowles et al., 2006). AMSR2, launched in May 2012 and mounted on the Japan Aerospace Exploration Agency's (JAXA's) Global Change Observation Mission-Water 1 (GCOM-W1) satellite, retains most of AMSR-E's physical properties, except for a larger antenna reflector and an extra C-band channel (7.3 GHz) (Cho et al., 2017). Both sensors share the same spatial resolution of 0.25° and provide measurements at 13:30 (ascending) and 01:30 (descending) local time at the Equator. This study utilized TB data from the 18.7 H and 36.5 V channels of the AMSR-E/2 Level-3 TB standard product to discriminate near-surface FT states at a 0.25° grid resolution. All available TB data from 2002 to 2023 were utilized, excluding the missing observations between AMSR-E and AMSR-2.

2.2 Global spatiotemporally continuous MODIS LST dataset

The spatiotemporally continuous MODIS LST dataset was utilized in this study, derived from two products: MODIS/Terra LST Daily L3 Global 0.05° CMG (MOD11C1) and MODIS/Aqua LST Daily L3 Global 0.05° CMG (MYD11C1). These products are separately obtained from the Terra and Aqua satellites, both of which cross the Equator at the local time of 13:30 and 01:30 during ascending and descending orbits, respectively. However, cloud contamination introduces data gaps in these products. To address this issue, a data interpolation and reconstruction method was applied, enabling the generation of spatiotemporally continuous LST records. Furthermore, the clear-sky LSTs were corrected to all-weather LSTs, enabling the retrieval of more realistic information (Yu et al., 2022; Zhao and Yu, 2021). These datasets under clear-sky and all-weather conditions exhibit satisfactory accuracy and are therefore suitable for high-resolution optical remote sensing inputs in the development of the downscaling algorithm.

2.3 GLASS Albedo

Spatiotemporally continuous land surface albedo products are essential for estimating global ATI, as outlined in the introduction. In this study, we used the GLASS02B06 dataset, which is part of the Global Land Surface Satellite (GLASS) project. The GLASS albedo series demonstrates accuracy comparable to that of the MODIS MCD43 albedo product, validated against FLUXNET site observations (Liu et al., 2013). Furthermore, GLASS albedo products address spatial gaps arising from cloud contamination and snow cover, thereby ensuring seamless and continuous products. This study employed the long-term clear-sky albedo dataset from GLASS02B06, which captures surface albedo at a 0.05° spatial resolution and provides new data every 8 d.

2.4 Land cover maps

In this study, the ancillary land cover maps were derived from the MODIS Land Cover Type CMG Yearly L3 Global 0.05° (MCD12C1) dataset (Friedl and Sulla-Menashe, 2022). With a spatial resolution of 0.05°, this dataset aligns with the resolution of the optical datasets mentioned earlier. It follows the 17-class International Geosphere-Biosphere Programme (IGBP) classification framework, as summarized in Table 1. The land cover dataset, illustrated in Fig. 1, was specifically utilized to mask out pixels of three IGBP land cover classes: water bodies, urban and built-up lands, and snow and ice. Furthermore, the corresponding type percentage dataset was used to filter out pixels dominated by large water bodies, which were then explicitly marked in the FT data record. Notably, no permafrost distribution data were used as input or constraints in generating the FT dataset.

Figure 1 IGBP global land cover map for 2019, depicting the spatial distribution of ground measurement sites within the ISMN and Tibet-Obs networks.

Table 1 IGBP land cover classification framework.

2.5 In situ soil temperature

Validation of the downscaled FT dataset was conducted using soil temperature measurements obtained from 41 dense observation networks and three sparse networks (SCAN, SNOTEL, USCRN). Among these, 42 networks were provided by the International Soil Moisture Network (ISMN), a global collaboration that delivers in situ measurements of soil moisture and related variables (Dorigo et al., 2013, 2021). Additionally, two networks (Naqu and Pali) were obtained from the Tibetan Plateau Observatory (Tibet-Obs), which focuses on plateau-scale monitoring of soil moisture and temperature (Su et al., 2011; Zhang et al., 2021a, b). Although the sampling time points vary across networks, the measurements are consistently conducted at hourly intervals. This study selected long-term in situ soil temperature data from 1027 stations within 44 global networks, all measured at a depth of 0–5 cm to match the penetration depth of the passive microwave observations. Figure 1 illustrates the spatial distribution of these stations.

2.6 Data pre-processing

The inter-calibration of different satellite instruments is for establishing consistent data records of Earth's environment (Chander et al., 2013). Although AMSR-E and AMSR2 share numerous similarities, calibrating these two sensors is necessary due to a 5 K difference in observed TB caused by different calibration procedures and other issues (Okuyama and Imaoka, 2015). Hu et al. (2019) introduced an inter-calibration linear model utilizing overlapping TB observations and the least squares method. This inter-calibration model was applied to TB data at the 18.7 and 36.5 GHz bands, generating a long-term passive microwave TB dataset. The inter-calibration model equations are as follows:

$$\begin{array}{}\text{(1)}& {\displaystyle }{\mathrm{TB}}_{\mathrm{AMSRE}\mathit{\_}\mathrm{18.7}H}=\mathrm{1.0189}\times {\mathrm{TB}}_{\mathrm{AMSR}\mathrm{2}\mathit{\_}\mathrm{18.7}H}-\mathrm{5.2717},\text{(2)}& {\displaystyle }{\mathrm{TB}}_{\mathrm{AMSRE}\mathit{\_}\mathrm{18.7}V}=\mathrm{1.0577}\times {\mathrm{TB}}_{\mathrm{AMSR}\mathrm{2}\mathit{\_}\mathrm{18.7}V}-\mathrm{16.2042},\text{(3)}& {\displaystyle }{\mathrm{TB}}_{\mathrm{AMSRE}\mathit{\_}\mathrm{36.5}H}=\mathrm{1.0073}\times {\mathrm{TB}}_{\mathrm{AMSR}\mathrm{2}\mathit{\_}\mathrm{36.5}H}-\mathrm{4.7723},\text{(4)}& {\displaystyle }{\mathrm{TB}}_{\mathrm{AMSRE}\mathit{\_}\mathrm{36.5}V}=\mathrm{1.0135}\times {\mathrm{TB}}_{\mathrm{AMSR}\mathrm{2}\mathit{\_}\mathrm{36.5}V}-\mathrm{6.3914},\end{array}$$

where TB_{AMSRE_18.7 H} and TB_{AMSR2_18.7 H} denote TB at 18.7 GHz in horizontal polarization obtained from AMSR-E and AMSR2, respectively. TB_{AMSRE_18.7 V} and TB_{AMSR2_18.7 V} correspond to vertical polarization at the same frequency. Other terms in the equations follow the same conventions.

To meet the requirements of daily ATI estimation, a linear interpolation method was employed to transform the GLASS 8 d shortwave clear-sky albedo datasets into daily data (Sohrabinia et al., 2014).

3 Methodology

3.1 FT discrimination from passive microwave observations

The DFA is a surface FT discrimination method developed using AMSR-E and AMSR2 data, demonstrating high accuracy compared to existing FT products. Near-surface FT variations are closely associated with soil temperature and moisture, which are reflected by the TB at 36.5 GHz in vertical polarization (TB_36.5 V) and the quasi-emissivity (Qe).

Qe is defined as the ratio of the TB at 18.7 GHz in horizontal polarization (TB_18.7 H) to TB_36.5 V. This ratio serves as an indicator of soil moisture, as microwave TB at these frequencies is highly sensitive to changes in water content. The near-surface soil FT process involves dynamic changes in soil and vegetation moisture states, which lead to significant variations in soil dielectric properties. Due to the substantial difference in dielectric constants between liquid water and ice, these phase transitions cause pronounced dielectric changes, a phenomenon further complicated by the distinct biogeochemical mechanisms of different vegetation types. The Qe index introduced in this study comprehensively characterizes the combined dielectric effects of both soil and vegetation during the FT process. Thus, Qe not only reflects soil moisture variations but also indirectly captures the biophysical status of vegetation and soil conditions, enabling consideration of the unique biogeochemical processes associated with various vegetation types during the FT process.

Therefore, TB_36.5 V and Qe were selected as key parameters for FT discrimination (Zhao et al., 2011).

The DFA is further parameterized to separately detect FT status during ascending and descending orbits (Wang et al., 2019b), as expressed by the following equations:

$$\begin{array}{}\text{(5)}& {\displaystyle }{\mathrm{FTI}}_{\mathrm{A}}=-\mathrm{0.123}\times {\mathrm{TB}}_{\mathrm{36.5}V}+\mathrm{11.842}\times \mathrm{Qe}+\mathrm{20.650},\text{(6)}& {\displaystyle }{\mathrm{FTI}}_{\mathrm{D}}=-\mathrm{0.209}\times {\mathrm{TB}}_{\mathrm{36.5}V}+\mathrm{9.384}\times \mathrm{Qe}+\mathrm{43.697},\text{(7)}& {\displaystyle }\mathrm{Qe}={\displaystyle \frac{{\mathrm{TB}}_{\mathrm{18.7}H}}{{\mathrm{TB}}_{\mathrm{36.5}V}}},\end{array}$$

where FTI_A and FTI_D represent the FT status for ascending and descending orbits, respectively. Through these equations, a long-term 0.25° FT record was generated using AMSR-E and AMSR2 TB data.

3.2 Estimation of apparent thermal inertia

ATI, as a substitute indicator for thermal inertia (TI), has proven effective in monitoring soil moisture conditions, thereby facilitating the downscaling of coarse-resolution FT states (Qin et al., 2013; Song and Jia, 2016; Van Doninck et al., 2011; Veroustraete et al., 2012; Verstraeten et al., 2006). ATI quantifies the temperature rise due to the Earth's surface absorbing radiant energy. Given that water has a relatively large heat capacity, higher soil water content provides greater resistance to external thermal fluctuations (Qin et al., 2013).

In this study, ATI is introduced as an indicator of soil moisture and is calculated as follows:

$$\begin{array}{}\text{(8)}& \mathrm{ATI}=C{\displaystyle \frac{\mathrm{1}-a_{\mathrm{0}}}{\mathrm{DTA}}},\end{array}$$

where a₀ is the actual surface albedo, derived from interpolated daily surface albedo data. C denotes the solar correction factor, which accounts for spatial and temporal variations in solar flux due to latitude and solar declination. The diurnal LST cycle's amplitude, denoted as DTA, represents the largest variation in LST observed within a single day.

C is calculated as follows:

$$\begin{array}{}\text{(9)}& \begin{array}{rl}C& =\sin \mathit{\phi }\sin \mathit{\delta }{\left(\mathrm{1}-{\tan }^{\mathrm{2}}\mathit{\phi }{\tan }^{\mathrm{2}}\mathit{\delta }\right)}^{\mathrm{1}/\mathrm{2}}\\ & +\cos \mathit{\phi }\cos \mathit{\delta }\arccos \left(-\tan \mathit{\phi }\tan \mathit{\delta }\right),\end{array}\end{array}$$

where φ is the latitude and δ represents the solar declination, calculated by

$$\begin{array}{}\text{(10)}& \begin{array}{rl}\mathit{\delta }& =\mathrm{0.006918}-\mathrm{0.399912}\cos \left(\mathrm{\Gamma }\right)\\ & +\mathrm{0.070257}\sin \left(\mathrm{\Gamma }\right)-\mathrm{0.006758}\cos \left(\mathrm{2}\mathrm{\Gamma }\right)\\ & +\mathrm{0.000907}\sin \left(\mathrm{2}\mathrm{\Gamma }\right)-\mathrm{0.002697}\cos \left(\mathrm{3}\mathrm{\Gamma }\right)\\ & +\mathrm{0.00148}\sin \left(\mathrm{3}\mathrm{\Gamma }\right),\end{array}\end{array}$$

where Γ represents the day angle, given by

$$\begin{array}{}\text{(11)}& \mathrm{\Gamma }={\displaystyle \frac{\mathrm{2}\mathit{\pi }\left(n_{\mathrm{d}}-\mathrm{1}\right)}{\mathrm{365.25}}},\end{array}$$

and n_d is the day number of the year.

The daily LST amplitude DTA in Eq. (8) can be calculated as

$$\begin{array}{}\text{(12)}& {\displaystyle \frac{\mathrm{DTA}}{\mathrm{2}}}={\displaystyle \frac{\begin{array}{c}n{\sum }_{i=\mathrm{1}}^n\cos \left(\mathit{\omega }{t}_i-\mathit{\psi }\right)T_i\\ -n{\sum }_{i=\mathrm{1}}^n\cos \left(\mathit{\omega }{t}_i-\mathit{\psi }\right){\sum }_{i=\mathrm{1}}^n{T}_i\end{array}}{\begin{array}{c}n{\sum }_{i=\mathrm{1}}^n{\cos }^{\mathrm{2}}\left(\mathit{\omega }{t}_i-\mathit{\psi }\right)\\ -{\left({\sum }_{i=\mathrm{1}}^n\cos \left(\mathit{\omega }{t}_i-\mathit{\psi }\right)\right)}^{\mathrm{2}}\end{array}}},\end{array}$$

where T_i is the LST recorded at time t_i and ω is the angular velocity of the Earth's rotation. The parameter n represents the number of LST observations for one pixel in a day (n≡4 in this study), and ψ is the phase angle, calculated as

$$\begin{array}{}\text{(13)}& {\displaystyle }\mathit{\psi }=\arctan \left(\mathit{\xi }\right)+\mathit{\pi },\text{(14)}& {\displaystyle }\mathit{\xi }={\displaystyle \frac{\begin{array}{c}\left(T_{\mathrm{1}}-T_{\mathrm{3}}\right)\left(\cos \left(\mathit{\omega }{t}_{\mathrm{2}}\right)-\cos \left(\mathit{\omega }{t}_{\mathrm{4}}\right)\right)\\ -\left(T_{\mathrm{2}}-T_{\mathrm{4}}\right)\left(\cos \left(\mathit{\omega }{t}_{\mathrm{1}}\right)-\cos \left(\mathit{\omega }{t}_{\mathrm{3}}\right)\right)\end{array}}{\begin{array}{c}\left(T_{\mathrm{2}}-T_{\mathrm{4}}\right)\left(\sin \left(\mathit{\omega }{t}_{\mathrm{1}}\right)-\sin \left(\mathit{\omega }{t}_{\mathrm{3}}\right)\right)\\ -\left(T_{\mathrm{1}}-T_{\mathrm{3}}\right)\left(\sin \left(\mathit{\omega }{t}_{\mathrm{2}}\right)-\sin \left(\mathit{\omega }{t}_{\mathrm{4}}\right)\right)\end{array}}},\end{array}$$

where Eq. (14) is used to calculate the diurnal LST cycle and requires clear-sky LST observations at four specific times during the day: Aqua/night (01:30), Terra/day (10:30), Aqua/day (13:30), and Terra/night (22:30). These observations are derived from the spatiotemporally continuous MODIS clear-sky LST products. The method of ATI calculation is adapted from Van Doninck et al. (2011). Consequently, a long-term 0.05° daily ATI dataset has been prepared for downscaling FT status.

3.3 Spatial downscaling of FT status

The FTI is characterized by its quantitative form, expressed in decimal values instead of binary values (Zhao et al., 2017). This feature enables its integration with other satellite-derived parameters, thereby enhancing its utility in remote sensing applications. As defined by soil temperature and emissivity from microwave data in Eqs. (5) and (6), the FTI exhibits a strong correlation with LST, thereby demonstrating its relevance in thermal analysis (Zhao et al., 2017). Moreover, the emissivity properties of soil are susceptible to changes in soil water content, which can be captured through ATI estimation.

Consequently, the assumed linear relationship between microwave and optical datasets is mathematically expressed as follows:

$$\begin{array}{}\text{(15)}& \mathrm{FTI}=a\cdot \mathrm{LST}+b\cdot \mathrm{ATI}+c,\end{array}$$

where the coefficients a, b, and c are determined through linear regression fitting. The MODIS/Aqua LST dataset was selected to ensure temporal consistency with the AMSR-E TB product.

Given the annual FT cycles, the downscaling approach was applied to the data on a yearly basis, following four steps as illustrated in Fig. 2.

Figure 2 Schematic diagram of spatial downscaling methodology.

The 0.05° LST and ATI data were resampled to align with the lower 0.25° resolution of the microwave-derived FTI. This resampling process involved averaging the data across 5×5 grid cells.

The FTI, LST, and ATI datasets were divided into yearly vectors, and the data within each vector were arranged in chronological order.

Linear regression fitting was performed on these three data vectors of each pixel, resulting in six coefficient matrices for a, b, and c at the ascending and descending times:

$$\begin{array}{}\text{(16)}& {\mathrm{FTI}}_{\mathrm{LR}}=a_{\mathrm{LR}}\cdot {\mathrm{LST}}_{\mathrm{LR}}+b_{\mathrm{LR}}\cdot {\mathrm{ATI}}_{\mathrm{LR}}+c_{\mathrm{LR}},\end{array}$$

where the abbreviations “HR” and “LR” refer to high and low resolution, respectively. Regression fitting was performed using only those dates with valid 0.25° FT data, together with the corresponding 0.05° optical data on those dates.

Once the coefficients a_LR, b_LR, and c_LR were expanded into 5 × 5 grids at a spatial resolution of 0.05°, they were applied in a regression model with the high-resolution (0.05°) LST and ATI data to generate the downscaled FTI values:

$$\begin{array}{}\text{(17)}& {\mathrm{FTI}}_{\mathrm{HR}}=a_{\mathrm{LR}}\cdot {\mathrm{LST}}_{\mathrm{HR}}+b_{\mathrm{LR}}\cdot {\mathrm{ATI}}_{\mathrm{HR}}+c_{\mathrm{LR}}.\end{array}$$

This process produced a finer spatial resolution compared to the original 0.25° FT maps derived from microwave data. Additionally, the land cover product was utilized to identify areas characterized by permanent snow and ice, substantial water cover, and urban or built-up lands, ensuring an accurate representation in the downscaled dataset.

3.4 Validation

Validating the downscaled FT product is essential to ensure its accuracy and reliability. In this study, soil temperature data from ISMN and Tibet-Obs ground stations were used to validate the downscaled FT results. Because satellite observations are divided into ascending and descending orbit data, the FT records were validated independently for each orbit. Frozen and thawed states in in situ soil temperature data are classified as follows:

$$\begin{array}{}\text{(18)}& \left\{\begin{array}{ll}T\le \mathrm{0},& \text{frozen}\\ T>\mathrm{0},& \text{thawed}\end{array}\right..\end{array}$$

To validate satellite data against ground measurements, the in situ observations must first be selected and pre-processed. In situ observations temporally closest to the satellite observation times are selected. ISMN observation data are recorded at hourly intervals, while the MODIS LST product provides observations at 20 min increments. Thus, in situ data recorded at the hour closest to the satellite observation time are used for validating FT classification. Missing data are excluded from the validation process. Due to the properties of microwave sensors, data gaps may occur in the 0.25° FT product, where the corresponding in situ observations are not included in the accuracy validation. Furthermore, only ISMN data labelled as “Good” quality is selected for validation.

Accuracy was calculated separately for the 0.25 and 0.05° FT products. The validation process, illustrated in Fig. 3, was conducted using a single, larger 0.25° grid cell as an example.

Figure 3 Schematic diagram of the validation method for comparing the FT product before and after downscaling.

FT conditions derived from in situ observations were considered the actual results. For a grid cell with a resolution of 0.25°, the average of all in situ soil temperature observations within the cell was calculated. The grid cell was subsequently classified according to the criteria delineated in Eq. (18), distinguishing between frozen and thawed states.

The FT classification results for the larger grid cell from both products were validated using in situ data. For the 0.25° product, the FT classification of the grid cell was directly compared to the in situ observation results. For the 0.05° product, each 0.25° grid cell contains 25 pixels. If the majority of these pixels (i.e. more than 13) were classified as frozen, the entire grid cell was labelled as frozen in the 0.05° product; otherwise, it was labelled as thawed.

The classification accuracies of the two products were calculated. Within each grid cell, the numbers of true and false FT classifications over the entire time series were counted to determine the discrimination accuracy using the following formula:

$$\begin{array}{}\text{(19)}& \mathrm{Accuracy}={\displaystyle \frac{\mathrm{FF}+\mathrm{TT}}{\mathrm{FF}+\mathrm{FT}+\mathrm{TF}+\mathrm{TT}}},\end{array}$$

where FF denotes the number of classifications where both the in situ observation and the satellite observation indicate a frozen state and FT denotes the number of classifications where the soil is frozen but the satellite observation indicates thawed conditions. TF and TT are defined similarly.

The same validation procedures were applied to all other grid cells containing in situ observations. The discrimination accuracy before and after downscaling was calculated and compared to assess any changes in accuracy resulting from the downscaling process.

3.5 Trend analysis

Trend analysis for global FT cycles is crucial for understanding temporal changes across different regions, which may indicate impacts associated with climate change. To analyse the annual global FT trend based on the downscaled product, this study employs two key parameters: the frost day and the frozen onset date.

A frost day is defined as the number of days within a year in which the lowest recorded temperature is less than 0 °C. The minimum temperature can be inferred to be below 0 °C if freezing is detected at the satellite's descending time for the same pixel. Therefore, in the satellite FT product, a frost day is identified as a day when the pixel is classified as frozen. In this study, the calculation and analysis of frost days rely on the descending FT product. The freeze onset date refers to the first day a pixel begins freezing and remains in a frozen state for more than 2 weeks. The initial day of this freezing period is designated as the freeze onset date. Both the frost day and the freeze onset date offer significant insights into the global distribution of near-surface FT variations. These two parameters are essential for comprehending the global temporal and spatial dynamics of freezing events.

To further analyse the trends of frost days and freeze onset dates, time series data for both parameters were generated for global-scale analysis using Sen's Slope and the Mann–Kendall (MK) test.

The equation for Sen's Slope is defined as

$$\begin{array}{}\text{(20)}& \mathit{\beta }=\mathrm{Median}\left({\displaystyle \frac{X_j-X_i}{j-i}}\right),\forall j>i,\end{array}$$

where X_j and X_j represent the data values at times i and j, respectively, and Median(⋅) represents the median calculation. This equation computes the median of the slopes between all pairs of data points. A positive β reflects an increasing trend, while a negative β indicates a decreasing trend.

The MK test is a statistical procedure that is employed to evaluate the null hypothesis that the data are independent and identically distributed (Mann, 1945; Kendall, 1948). It is employed to detect trends in time series data without assuming a linear relationship. Given a series X_t of length n, the null hypothesis posits that the values in X_t are independent. The MK test statistic S is calculated as

$$\begin{array}{}\text{(21)}& S=\sum _{i=\mathrm{1}}^{n-\mathrm{1}}\sum _{j=i+\mathrm{1}}^n\mathrm{sgn}\left(X_j-X_i\right),\end{array}$$

where n is the sample size, defined as the number of data points, and sgn(X_j−X_i) is the sign function, defined as

$$\begin{array}{}\text{(22)}& \mathrm{sgn}\left(X_j-X_i\right)=\left\{\begin{array}{ll}\mathrm{1},& \mathrm{if}{X}_j>X_i\\ \mathrm{0},& \mathrm{if}{X}_j=X_i\\ -\mathrm{1},& \mathrm{if}{X}_j<X_i\end{array}\right..\end{array}$$

The test statistic S evaluates the cumulative number of positive differences that exceed negative differences. In this study, because the frost day and freeze onset date of each year in the same pixel are unique, the variance of S is computed as

$$\begin{array}{}\text{(23)}& \mathrm{Var}\left(S\right)={\displaystyle \frac{n(n-\mathrm{1})(\mathrm{2}n+\mathrm{5})}{\mathrm{18}}}.\end{array}$$

If the sample size n>10, the standard normal test statistic Z is calculated as

$$\begin{array}{}\text{(24)}& Z=\left\{\begin{array}{ll}\frac{S-\mathrm{1}}{\sqrt{\mathrm{Var}\left(S\right)}},& \mathrm{if}S>\mathrm{0}\\ \mathrm{0},& \mathrm{if}S=\mathrm{0}\\ \frac{s+\mathrm{1}}{\sqrt{\mathrm{Var}\left(S\right)}},& \mathrm{if}S<\mathrm{0}\end{array}\right..\end{array}$$

Trend testing is conducted at a specified significance level, α. If the absolute value of Z exceeds the critical value $Z_{\mathrm{1}-a/\mathrm{2}}$, the null hypothesis of no trend is rejected. In this study, the test was performed at α=0.05, corresponding to a 95 % confidence interval. Under these conditions, if |Z|>1.96, the null hypothesis is rejected. The results of the test are categorized into five trend types, as summarized in Table 2. Sen's Slope quantifies the direction of temporal variation as positive or negative, while the MK test evaluates the statistical significance of these trends. This combined approach allows for the assessment of both the magnitude and significance of trends in frost days and freeze onset dates, providing insights into changes in global FT cycles that may be linked to climate change.

Table 2 Classification of trends based on the MK test.

4 Results and discussions

4.1 Result of downscaled FT product

According to the methodology described in Sect. 3.3, a global near-surface FT states dataset with a 0.05° resolution was generated, spanning the years 2002 to 2023. This dataset provides daily FT information for both daytime and nighttime, corresponding to the satellite's ascending (13:30) and descending (01:30) orbits, respectively. The left panel of Fig. 4 presents a comparison between the original 0.25° FT discrimination product and the downscaled FT product for 1 January 2019 at 01:30 (descending orbit). The dataset classifies land surface conditions into the following categories: thawed soil, frozen soil, snow and ice, urban and built-up lands, water bodies, regions with missing data, and water-influenced areas (labelled as “water” in Fig. 4).

Figure 4 Maps of the FT product for (a) the global scale and (b) the Qinghai-Tibet Plateau at 01:30 on 1 April 2019, derived from microwave data with a 0.25° spatial resolution. Additionally, panels (c) and (d) display the same conditions in the downscaled FT product at a higher resolution of 0.05°.

An important point to consider is that areas near water bodies may be misclassified as frozen due to the high soil moisture content detected by microwave TB data. However, these areas are unlikely to freeze, particularly in coastal regions where temperatures consistently remain above 0 °C. To improve the reliability of the FT product, these water-influenced areas should be carefully delineated and flagged.

In addition to the classification of surface conditions, another important feature of the dataset is its treatment of missing data. In the original 0.25° FT product, data gaps occur due to the inherent swath gaps of the AMSR-E/2 satellites, resulting in missing values for some grid cells on certain days. By contrast, our downscaling approach for the 0.05° product effectively overcomes this limitation. For each 0.25° pixel, regression fitting is performed using only those dates with valid 0.25° FT data, together with the corresponding 0.05° optical data on those dates. As the 0.05° optical data are spatiotemporally continuous, regression can always be performed on any date with valid 0.25° FT observations. The fitted regression model is then applied to the complete year of 0.05° optical data, resulting in a seamless, gap-free 0.05° FT dataset. This missing-data treatment in the downscaling process ensures that the high-resolution FT product contains no data gaps, as demonstrated in Fig. 4. The right panel of Fig. 4 presents a comparison of the FT product for the Qinghai-Tibet Plateau, highlighting the enhanced resolution of the 0.05° product in capturing FT states. The downscaled product more precisely classifies soil freezing and thawing states, as well as specific surface features that are more difficult to detect in the 0.25° product. This demonstrates the enhanced capability of the 0.05° product to capture FT cycles at the local scale, compensating for the resolution limitations of the original 0.25° product.

This downscaled product, with its enhanced spatial resolution, allows for a more detailed analysis of global FT dynamics, providing a more reliable interpretation of global FT patterns. This capability is crucial for understanding FT cycles and their associated environmental impacts.

Figure 5 Maps of frost days in April 2019 for (a) the global scale and (b) the Qinghai-Tibet Plateau, derived from microwave data with a 0.25° spatial resolution during the descending orbit time. Additionally, panels (c) and (d) show the same conditions in the downscaled descending FT product at a resolution of 0.05°.

Figure 5 illustrates frost days for April 2019 at two spatial resolutions, 0.25 and 0.05°, for both the global scale and the Qinghai-Tibet Plateau.

The left panel compares global frost days for a single month, showing that the 0.25° product exhibits a more gradual variation in frost days across the Northern Hemisphere. In contrast, the 0.05° product reveals more detailed and intricate patterns of frost day variation. For instance, in southern Russia and parts of North America, the 0.05° product captures more nuanced frost day variations, highlighting significant differences in frost day counts between neighbouring grid cells. Meanwhile, the 0.25° product displays relatively uniform changes on frost days, reflecting less variation due to its coarser resolution.

The right panel of Fig. 5 compares frost days for the Qinghai-Tibet Plateau, illustrating finer variations at the higher spatial resolution. This demonstrates that the downscaled product provides a more detailed distribution of frost days while also capturing the trend of frost day variation at a finer spatial scale.

In summary, the downscaled FT product offers significantly more detailed information, especially in regions with complex FT dynamics. This improvement in spatial resolution is crucial for our research, as it enables a more refined understanding of global FT processes and their environmental impacts. The ability to capture finer details and regional variations is a key advantage of the downscaling approach, which effectively meets the primary objectives of this study.

4.2 Validation with in situ soil temperature

In this section, the coarse- and high-resolution FT products are validated against in situ soil temperature data from 44 networks. The scatter plots of accuracy across all ground stations are shown in Fig. 6, with detailed validation results provided in Table 3. The downscaled products achieved overall accuracies of 87.63 % and 83.78 % at ascending and descending orbit times, respectively. These accuracies are comparable to those of the original products, which had overall accuracies of 87.72 % and 84.08 %. In addition, Table 4 summarizes the classification accuracies of both products across various land cover types for both ascending and descending passes. The land cover types are defined by the ESA CCI Land Cover 2010 classification values provided by the ISMN sites, enabling a more comprehensive assessment of performance under different physiographic conditions.

Figure 6 Scatter plots comparing the downscaled FT product at a 0.05° resolution with the FT product at a 0.25° resolution. The plots show data for (a) ascending orbits and (b) descending orbits. The colour bar represents the density of data points.

Table 3 Validation results for the coarse- and high-resolution FT products across various networks. Networks 1 to 27 typically experience FT cycles, whereas networks 28 to 44 do not exhibit FT phenomena. “Num” represents the number of FT product samples that matched in situ measurements at the corresponding sampling times and were used for validation.

Table 4 Validation results for the coarse- and high-resolution FT products across various land cover types. Land cover types are sourced from the ESA CCI Land Cover 2010 classification values provided within the ISMN site dataset.

For networks experiencing FT cycles, the accuracies of the high-resolution products are 85.91 % and 81.55 % at ascending and descending orbit times, respectively. For networks without FT phenomena, the accuracies are 97.43 % and 97.97 %. The colour gradient in the plots of Fig. 6 reflects the data density, with points shifting toward yellow as the density increases. Additionally, in 65.91 % and 56.92 % of networks, accuracy exceeds 90 % at both ascending and descending orbit times for both FT products.

These findings suggest that the data points are most densely concentrated in the range where accuracy values lie between 0.9 and 1. Both ascending and descending orbit results closely align with the one-to-one line, indicating that the accuracy of the FT classification remains largely unchanged after the downscaling process. Despite the enhanced spatial resolution, the downscaled product preserves the accuracy and consistency of the original resolution.

These results are consistent with the conventional understanding of downscaling, which aims to enhance spatial resolution without compromising the accuracy of the original product. As a result, the downscaled FT product provides higher spatial resolution and more detailed information while maintaining the integrity and quality of the original observations.

4.3 FT dynamics and trends

4.3.1 Trend analysis of frost days

After validating the FT dataset's accuracy with in situ observations, we calculated the average number of frost days and their trends over the period 2003–2023. Figure 7a illustrates the spatial distribution of the annual average frost days in the Northern Hemisphere. In high-latitude regions (north of 45° N), the average number of frost days is approximately 187.8 ± 12.7 d (spatial standard deviation). In contrast, frost days in lower-latitude regions vary due to seasonal shifts.

Figure 7 Northern Hemisphere maps showing (a) average annual number of frost days and (b) trend analysis results for the 0.05° downscaled FT product during the descending satellite pass from 2003 to 2023.

While frost days generally increase with latitude, this pattern is not globally consistent. For example, the Qinghai-Tibet Plateau, despite its relatively low latitude, experiences a higher number of frost days due to its high elevation. In the Southern Hemisphere, soil freezing is rare, except in localized regions along the Andes Mountains, where freezing occasionally occurs under specific climatic conditions.

The spatial trends of annual frost days are shown in Fig. 7b. MK trend analysis identifies a decreasing trend in approximately 14.35 % of global land areas, with 2.67 % exhibiting statistically significant declines. This trend reflects the impact of climate warming, especially across much of the Eurasian continent and Alaska, with particularly obvious effects in Russia and the Qinghai-Tibet Plateau.

In contrast, about 11.17 % of the global land area shows an increasing trend in frost days, of which 1.55 % are statistically significant. Regions with increasing frost days include North America and West Asia. These opposing trends underscore the complexity of regional climate dynamics, revealing that while many areas are warming with fewer frost days, localized cooling in some areas results in more frequent freezing events.

The analysis of frost day trends demonstrates the utility of the FT dataset while providing valuable insights into regional variations in climate change impacts.

To further explore the climatic and geocryological significance of this metric, we compared the spatial distribution of annual frost days in 2017 with independently derived permafrost maps over the Qinghai-Tibet Plateau (Zhao, 2017; Zou et al., 2017), as illustrated in Fig. 8. This qualitative comparison reveals a notable spatial agreement, demonstrating that regions with high annual numbers of frost days, as detected by remote sensing, are generally consistent with areas classified as permafrost in existing reference datasets. Specifically, the average annual number of frost days within permafrost-classified pixels is approximately 278.85, indicating a potential spatial correspondence between these two metrics. This finding highlights the potential of the annual number of frost days as a valuable proxy for assessing permafrost extent and its spatial variability.

Figure 8 Maps of (a) permafrost distribution (Zhao, 2017; Zou et al., 2017) and (b) annual frost days derived from the downscaled descending FT product in 2017 over the Qinghai-Tibet Plateau.

4.3.2 Trend analysis of freeze onset date

A statistical analysis was performed to examine the average freeze onset date and its trend from 2003 to 2023. Figure 9a shows the spatial distribution of the annual average freeze onset date across the Northern Hemisphere. Soil freezing predominantly occurs in the Northern Hemisphere, where the average freeze onset date in high-latitude regions is approximately 240.3 ± 7.2 d (spatial standard deviation). A general trend is observed, with freezing beginning earlier in higher-latitude regions compared to lower-latitude areas. However, this pattern is not uniform. On the Qinghai-Tibet Plateau, the high elevation causes freezing to occur earlier, leading to lower freeze onset dates. In the Southern Hemisphere, soil freezing is limited to small areas, primarily along the Andes Mountains, where freezing also begins earlier in the year.

Figure 9 Northern Hemisphere maps showing (a) annual average freeze onset date and (b) trend analysis results for the 0.05° downscaled FT product from 2003 to 2023.

Figure 9b illustrates the spatial patterns of trends in freeze onset dates. The MK trend analysis identifies a decreasing trend in freeze onset dates for approximately 9.10 % of the global land area, of which 1.22 % exhibits statistically significant reductions. A decrease in freeze onset date indicates an earlier start to freezing in these regions, with the most pronounced trends observed in eastern Russia. Conversely, approximately 7.57 % of the global land area demonstrates an increasing trend in freeze onset dates, with 0.91 % exhibiting statistically significant increases. An increase in freeze onset date reflects a delayed start to freezing. These regions are more geographically dispersed, reflecting the widespread influence of global warming on soil FT dynamics.

Investigating frost day and freeze onset trends demonstrates the applicability of the FT dataset and provides deeper insight into regional variations in climate change effects on global FT patterns.

4.4 Discussion

4.4.1 Enhanced capabilities of high-resolution FT products

This study utilizes passive microwave and optical observations to achieve high-resolution detection of FT states. Initially, passive microwave sensors were used to monitor global surface FT states, generating an original FT dataset at a 0.25° resolution. Due to their minimal susceptibility to atmospheric conditions, such as cloud cover and aerosols, these sensors enable continuous tracking of FT transitions across diverse landscapes and climatic zones globally. Furthermore, the associated discrimination algorithm demonstrates high classification accuracy, rendering it suitable for global-scale applications. However, the original dataset's coarse spatial resolution of 0.25°, combined with inherent seams in microwave sensor data, limits its capacity to capture detailed surface variations. In contrast, optical observations, including LST and albedo, offer higher spatial resolution and provide more detailed surface information. Additionally, the development of seamless, long-term data products not only addresses the data seams in the original dataset but also maintains a short revisit interval, enabling daily surface monitoring.

Therefore, to integrate the advantages of both data sources, a bivariate regression model was developed, effectively downscaling the original 0.25° resolution dataset to a finer 0.05° resolution. This improved resolution allows for a more detailed representation of FT states, facilitating analyses of regional FT dynamics, particularly in areas with complex terrain, diverse vegetation, or significant seasonal variations. As shown in Fig. 10, the 0.05° resolution FT product reveals finer details of FT cycles that were previously undetectable in the 0.25° dataset. These enhancements are critical for research requiring precise identification of localized surface condition changes, such as those affecting soil moisture and ecosystem carbon dynamics.

Figure 10 Time series plots comparing the 0.25° FT product values and the 0.05° frozen ratios in the downscaled product at the SQ21 station of the NGARI network, distinguished by (a) daytime and (b) nighttime observations.

The methodology employed in this study not only eliminates seams in the original FT dataset but also integrates data from multiple sources, including MODIS optical products and microwave TB observations. This fusion provides a robust and continuous long-term FT record, essential for understanding global FT cycles and their environmental implications. Moreover, the enhanced spatial resolution enables detailed analyses of changes in local FT conditions, offering valuable insights for hydrological modelling, climate research, and ecosystem management. These advancements expand the applicability of the FT record, contributing significantly to diverse scientific fields.

4.4.2 Uncertainties associated with ground validation

Ground temperature data serve as a reliable basis for validating satellite-derived FT products. In this study, soil temperature measurements at 0–5 cm depth were used for validation rather than surface-level (0 cm) measurements due to their superior representation of near-surface soil FT states. Additionally, the utilization of 0–5 cm depth data ensures comparability with satellite-based FT discrimination, as the typical penetration depths of passive microwave observations at 18.7 and 36.5 GHz are approximately 3–5 and 0–2 cm, respectively, both under 5 cm. While 0 °C is commonly employed as the threshold for distinguishing frozen and thawed states in validation using in situ temperature data, this criterion may not always accurately reflect actual ground conditions. Various factors, including wind, vegetation cover, and snow cover, can influence the true ground state even when surface temperatures are at or below 0 °C.

Currently, there is no universally accepted standard for defining absolute temperature boundaries for frozen and thawed ground. However, the use of 0–5 cm soil temperature remains the most physically meaningful for passive microwave FT validation, as it provides a more direct and reliable indicator of the near-surface state than either air temperature or deeper soil temperature measurements. Despite this, under certain conditions, such as areas with strong winds or heavy snow cover, this approach may lead to classification errors, particularly when the ground thawing process is disrupted at temperatures near 0 °C. To address these limitations, future studies should concentrate on identifying optimal soil depths and establishing more precise temperature thresholds to enhance validation accuracy.

4.4.3 Limitations and challenges

Although this study has effectively generated high-resolution FT data by integrating passive microwave observations with optical datasets, further refinements are required to address the remaining challenges. First, despite the resolution enhancement to 0.05°, it is still relatively coarse for applications demanding finer-scale detail. In complex surface environments, this resolution may not adequately capture localized variations in FT states. To reduce the influence of spatial heterogeneity, higher-resolution optical observations are necessary. For instance, utilizing a 1 km LST and albedo dataset for downscaling would align the resulting FT record with this finer spatial resolution. While 1 km LST datasets are available, the development of more stable and higher-resolution LST products remains critical. The quality of the downscaled FT dataset is closely reliant on the quality of the optical observations employed. Future research should prioritize the use of seamless, high-quality optical datasets with short revisit intervals and higher spatial resolution to produce more precise and reliable FT classification products.

Second, the accuracy of the downscaled FT dataset is inherently determined by FT discrimination algorithms, which rely on lower-resolution microwave remote sensing. As outlined in Sect. 4.2, the downscaling process preserves the accuracy characteristics of the original FT data. Therefore, improvements in the FT discrimination algorithms are necessary for achieving further gains in classification accuracy.

Additionally, during the data fusion process, combining optical observations with microwave data, various factors may affect the correlation coefficients and the performance of the linear regression model. These factors include land cover types, vegetation conditions, terrain elevation, and seasonal variations, all of which can introduce instability in model performance and impact FT classification accuracy. To address these limitations, future studies should incorporate additional surface parameters, including the normalized difference vegetation index (NDVI), to improve the model's adaptability and accuracy across diverse surface conditions.

5 Data availability

The global 0.05° near-surface soil FT state dataset (FT-HiDFA, 2002–2023) is freely available at the National Tibetan Plateau/Third Pole Environment Data Center. The dataset can be accessed via https://doi.org/10.11888/Cryos.tpdc.301551 (Zhao et al., 2024b) or https://cstr.cn/18406.11.Cryos.tpdc.301551 (last access: 28 January 2025).

6 Conclusion

This study developed a comprehensive, high-resolution global dataset of surface FT states by integrating multi-source remote sensing data. The primary objective was to enhance spatial resolution and provide more detailed global monitoring of FT states. Coarse-resolution FT data were initially derived from passive microwave TBs, while long-term ATI data were calculated using optical observations. By leveraging the complementary strengths of passive microwave and optical remote sensing products, a high-resolution daily near-surface FT dataset was produced, offering a finer representation of surface FT conditions.

Subsequently, the coarse- and high-resolution FT products were validated using ground-based in situ observations. This validation facilitated a thorough evaluation of accuracy changes associated with the downscaling process and enabled a comprehensive trend analysis on a global scale using the enhanced-resolution FT records.

The findings have significant implications for understanding and monitoring ecological and hydrological responses to climate dynamics. The analysis revealed intricate patterns of frost days and freeze onset dates, demonstrating the varying regional influences of climate change. The downscaled product, with its enhanced spatial resolution, provided detailed insights into FT dynamics, which are crucial for studies requiring precise identification of local surface condition changes, such as those impacting soil moisture or ecosystem carbon dynamics.

This integration of multi-source remote sensing data marks a significant advancement in monitoring Earth's surface processes. It demonstrates the potential to improve climate models and other environmental assessments. Furthermore, the proposed methodological framework is adaptable for similar applications globally, enhancing the predictive capabilities for climate-related phenomena and supporting environmental decision-making. By providing a robust and continuous long-term FT record, this study contributes to the development of hydrological models, climate studies, and ecosystem management, thereby expanding the applicability of the record in diverse scientific fields.