Development of a Pixel-Wise Forest Transmissivity Model at Frequencies of 19 GHz and 37 GHz for Snow Depth Inversion in Northeast China

Abstract

Satellite passive microwave remote sensing has been extensively used to estimate snow depth (SD) and snow water equivalent (SWE) across both regional and continental scales. However, the presence of forests causes significant uncertainties in the estimations of snow parameters. Forest transmissivity is one of the most important parameters for describing the microwave radiation and scattering characteristics of forest canopies. Although many researchers have constructed models for the functional relationship between forest transmissivity and forest vegetation parameters (e.g., stand growth and accumulation), such relationships are strongly limited by the inversion accuracy of vegetation parameters, forest distribution types, and scale-transformation effects in terms of regional or global scale applications. In this research, we propose a pixel-wise forest transmissivity estimation model (Pixel-wise γ Model) based on long-term series satellite brightness temperature (TB) data for the satellite remote sensing inversion of snow parameters. The model performance is evaluated and applied in SD inversion. The results show that the SD inversion errors RMSE and Bias are 9.8 cm and −1.5 cm, respectively; the SD inversion results are improved by 41% and 84% after using the Pixel-wise γ Model, compared with the forest transmissivity model applied in the GlobSnow v3.0 product. The proposed forest transmissivity model does not depend on forest cover parameters and other ground measurement parameters, which greatly improves its application scope and simplicity.

1. Introduction

The high-precision inversion of snow depth (SD) and snow water equivalent (SWE) in forest areas is a difficult and open problem in the field of snow remote sensing. Optical remote sensing has advantages for high-resolution mapping of snow covered areas, but it is powerless in the estimation of SD/SWE, and even more so, for the detection of snow under forests, and is also susceptible to cloud influence [1,2]. Passive MicroWave (PMW) remote sensing has unique advantages for SD and SWE inversion and mapping, but mainstream PMW snow remote sensing products, such as the Advanced Microwave Scanning Radiometer 2 (AMSR-2) SD product from the Global Change Observation Mission 1st-Water [3] and the Global Snow Monitoring for Climate Research (GlobSnow) SWE product from the European Space Agency [4,5,6], have not achieved satisfactory accuracy in forest areas [7,8], and cannot meet the requirements of meteorological, hydrological and climate models. Forest snow has a very complex composition; in addition to snow, it also includes dielectric layers such as surface litter, tundra vegetation, shrubs (clumps), trunk (canopy), branches and leaves. The interaction between layers, as well as their radiation and scattering, are factors of interference modeling and signal decoupling. Therefore, in order to realize the high-precision inversion of snow parameters in forest areas, it is necessary to have a clear understanding of the radiation transmission process and the emission and scattering mechanisms of forest snow coupled systems. Quantitatively exploring the relationship between the geometric parameters (such as tree height, diameter at breast height (DBH), canopy density, gap area, etc.,) and physical parameters (such as temperature, etc.) of the forest canopy medium layer and the characteristics of microwave radiation (such as polarization, frequency, etc.), is helpful to further improve and develop the relevant Radiative Transfer (RT) model and provide theoretical support for the inversion of snow parameters in forest areas.

Snow cover is an important scatterer in PMW remote sensing observation, which attenuates the underlying surface microwave signal received by the sensor. Various other components (SD, snow density, frozen/non-frozen soil, etc.), especially in forest parameters, also affect the brightness temperature measured by satellite sensors, which gives the remote sensing inversion of snow parameters under the forest great uncertainty, which in turn seriously limits the accuracy of traditional algorithm inversions of SD and SWE. Therefore, the decoupling of the forest canopy signal and the underlying surface signal (snow, frozen/non-frozen etc.) is particularly important in microwave radiation modeling. In earlier research, the influence of forest terrain on the inversion accuracy of snow parameters has generally been reduced by simply introducing forest coverage parameters into the algorithm as a correction factor [9,10]. However, due to the simplification of the physical process, the universality of the algorithm is unsatisfactory. In recent years, scholars have made some progress in developing and applying empirical or semi-empirical models (including look-up table method, Bayesian assimilation method, etc.) and complex models considering canopy and trunk scattering/absorption effects [11,12,13], but subject to the acquisition of accurate information parameters of vegetation, this is still a significant challenge in the application of satellite remote sensing inversion [14]. Therefore, simplified RT models, such as τ-ω models, are often used to characterize forest radiative transfer processes [15,16,17]. In PMW satellite signals of the forest snow observation scene, forest radiation is a function of microwave frequency, polarization mode, satellite observation angle and forest effective physical temperature, etc. Among these parameters, forest canopy microwave transmissivity (hereinafter referred to as forest transmissivity) is one of the important parameters for describing the absorbing and scattering characteristics of forests. The question of how to parameterize the forest transmissivity accurately is one of the key points to correct the forest effect, and to realize the high-precision inversion of snow parameters under forests.

Among the SD/SWE inversion algorithms based on RTM in forest areas, the previous transmissivity parametric models are often obtained by establishing empirical or semi-empirical functional relationships based on point scale ground measurement data. For example, microwave transmissivity of forest vegetation has been explored through observation [18,19,20,21,22,23,24] and modeling perspectives [25,26,27], respectively. However, scale conversion is rarely considered; that is, when the point scale is extended to the PMW pixel scale, the applicability of the model is uncertain. For example, although the transmissivity model expressed as a function of Growing Stock Volume (GSV) is applied to GlobSnow products, its uncertainty still increases exponentially with the increase of forest coverage [14]. In view of the above limitations, we need to consider how to parameterize or solve the forest transmissivity directly from the pixel scale.

Zeroth order τ–ω models [15,16,17] aim to simplify complex microwave – tree interactions into robust, tractable equations, requiring only a small number of parameters, and thereby making them suitable for large-scale applications, especially when there is insufficient auxiliary data, such as canopy structure, that can support higher-order models. Based on the τ–ω model, this study develops a Pixel-wise model to estimate the winter forest transmissivity of horizontal (H) and vertical (V) polarizations at frequencies of 19 GHz and 37 GHz for each year of 2014–2019 in forest area of Northeast China. The model performance is evaluated, and forest transmissivity data sets are applied in SD inversion, and satisfactory results are obtained.

2. Study Area and Data Set

2.1. Study Area

Northeast China (38°72′–53°55′ N, 115°09′–135°52′ E) is located in the mid-latitude of the northern hemisphere, covering a total area of approximately 1.2 million km² and with an elevation range of 0–2667 m. The climate of Northeast China is controlled by the high latitude East Asian monsoon, which includes warm temperate and cold temperate zones latitudinally, as well as semi-arid, semi-humid and humid zones from west to east. The geographical location and climatic conditions of this region creates a relatively complex surface cover system. According to the GLC_FCS30-2015 global land-cover products [28] (free access at https://doi.org/10.5281/zenodo.3986872, accessed on 8 July 2020), Northeast China is dominated by forests, croplands and grasslands, as shown in Figure 1a, accounting for about 40.4%, 36.9% and 13.1% of the total area, respectively. Forest cover is widely distributed in the mountainous areas, including the Greater Khingan Mountains, Lesser Khingan Mountains and Changbai Mountains areas [29]. The forest cover fraction data within each passive microwave satellite pixel (25 km × 25 km) are calculated (Figure 1b). The areas of deciduous broad-leaved forests and deciduous coniferous forests of forestry land in Northeast China are superior, which are 73.1% and 22.4% of the total, respectively. In addition, Northeast China is one of the three primary snow-covered regions (Xinjiang, Northeast China and the Qinghai-Tibet plateau) in China [30,31]. In winter, snowfall is frequent and evenly distributed, and snow cover lasts from November until April. The maximum snow depth (MSD) in Northeast China forest area generally shows a spatial distribution characteristic, gradually deepening from southwest to northeast, with variations ranging from 5–9 cm incrementally to 21–32 cm [32]. Rich snow resources and forest diversity make Northeast China an ideal region for forest snow research.

2.2. Input Data for Forest Microwave Transmissivity Modeling

The investigation period is from September 2014 to March 2019, covering five winter freezing seasons across Northeast China.

2.2.1. SSMIS Brightness Temperature

Time series TB data from the SSMIS sensor currently operating on board the Defense Meteorological Satellite Program (DMSP) spacecraft F17 were used for solving the transmissivity of 19 GHz and 37 GHz in Northeast China’s forested areas. The TBs at 19, 37 and 91 GHz were obtained from the National Snow and Ice Data Center (NSIDC, https://nsidc.org/data/NSIDC-0032, accessed on 31 July 2019) at the EASE Grid format, with a nominal resolution of 25 km × 25 km [33]. The dynamic range of TBs is 55–320 K, and the ground observation angle is 53.1°. To match the spatial scale of different remote sensing products, the original TB data set was resampled into a standard latitude and longitude gridding (0.25° × 0.25°) format. Only descending (approximately 05:00–07:00, local time) orbital data was applied for reducing the uncertainty caused by wet snow and atmospheric turbulence. The TB values with sudden unrealistic temporal changes had been filtered out according to the criterion $d\mathrm{TB}>3\cdot \mathrm{Std}(d\mathrm{TB})$, meaning TBs are chosen as normal values within three standard deviations.

2.2.2. ERA5-Land Daily Air Temperature Data

ERA5-Land is an updated version repeated after a series of improvements based on the reanalysis data of ERA5, with higher spatial resolution (9 km in ERA5-Land vs 31 km in ERA5), and the time span has been from January 1950 to now. The ERA5-land product was generated on the basis of the ERA5 climate dataset by the European Centre for Medium-range Weather Forecasts (https://cds.climate.copernicus.eu/, accessed on 10 June 2021). In this study, the daily air temperature elements (Unit: K) at 2 m from ERA5-Land data set were used, which refers to the air temperature at 2 m above the surface of land and is obtained by interpolating between the lowest layer of the model and the Earth’s surface. Due to the different spatial resolutions of SSMIS (standard latitude and longitude gridding 0.25° × 0.25°) and ERA5-Land (0.1° × 0.1°), the ERA5-Land air temperature product was resampled to the same spatial size as SSMIS. The processing operations of re-griding were necessary for aligning all relevant information sources, including forest cover fraction, land cover type, ERA5-Land derived air temperature, TB observations and Pixel-wise transmissivity model outputs on a coincident grid.

2.3. Validation Data

2.3.1. Microwave Transmissivity Measurement Data

In mid-to-late December 2016, we conducted a series of ground-based microwave radiometric observation experiments in the forest. Downwelling TBs at 19- and 37-GHz for horizontal and vertical polarizations from the forest canopy were measured at 14 tree sample plots using microwave radiometers. The geographical location of the radiometric sites are shown in the red triangle in Figure 1b. Based on the measured data of TBs and air temperature, combined with the sky background TBs data, the in-situ forest transmissivity for each frequency and polarization was calculated, following the simplified approach provided by [19]. A further detailed description of this dataset can be found in [24]. In this study, this dataset is considered as a reference for the actual value of ground measurements of forest transmissivity.

2.3.2. China’s Snow Survey Data

China’s snow survey (also known as Investigation on Characteristics and Distribution of Snow Cover in China) data sets were collected from a total of six snow measuring lines, as well as several quadrats and stations, spanning across Northeast China, Xinjiang, and Qinghai-Tibet plateaus, during the winter period of 2017–2019 [34,35]. In this study, only the line measurement dataset in Northeast China was used, and the experimental area was mainly distributed around the Greater Khingan Mountains, Lesser Khingan Mountains, Sanjiang Plain, Songneng Plain and Changbai Mountains. The entire snow season was divided into three periods, with reference to snowfall frequency and the surface snow accumulation state: accumulation period, stabilization period and melting period, which lasted for two years; the observations were repeated three times a year. The distribution of snow pit sampling sites is illustrated via the blue ink dots in Figure 1b. As of 18 March 2019, 443 snow pit sample sites were collected in Northeast China, of which 329 were located in SSMIS pixels, meeting the condition of forest cover fraction ≥10%. Snow pit profile measurement parameters mainly include SD, snow grain size, snow layer dielectric constant, snow liquid water content, snow layer hardness, snow layer temperature, stratified snow density and overall snowpack density.

3. Models and Methodology

3.1. Zeroth Order τ–ω Radiative Transfer Model

The TB above the forest canopy at certain nadir angles (e.g., 53.1° for SSMIS observations), and with different frequencies and polarizations, can be modeled using a simple microwave vegetation radiative transfer model called the τ–ω model, which calculates the following main contributions (ignoring the subtle atmospheric emission reflected by the underlying surface) as:

$$T_{B,ac}^{f,p}=T_{B,\mathit{underlying}}^{f,p}+T_{B,\mathit{forest}}^{f,p}+T_{B,forest-r}^{f,p}$$

(1)

where $T_{B,ac}^{f,p}$ is the upwelling microwave emission (unit: K) above forest canopy measured by satellite sensor; $T_{B,\mathit{underlying}}^{f,p}$ is the upwelling emission (unit: K) of the underlying surface (containing soil and snow); $T_{B,\mathit{forest}}^{f,p}$ is the upwelling emission (unit: K) of the forest canopy; $T_{B,forest-r}^{f,p}$ is the downwelling forest vegetation emission (unit: K) reflected by the underlying surface; superscript f represents frequencies 19 and 37 GHz in this study; superscript p is the vertical (V) or horizontal (H) polarization. Following the theory and method proposed by Mo et al. (1982) [15], with vegetation scattering and transmission coefficient (τ–ω) and ground emissivity, Equation (1) can be further expressed as:

$$T_{B,ac}^{f,p}=\gamma \cdot E_{}^{f,p}\cdot T_{underlying}+\left(1-{\omega }^{f,p}\right)\left(1-{\gamma }^{f,p}\right)\cdot T_{\mathit{forest}}+\left(1-{\omega }^{f,p}\right)\left(1-{\gamma }^{f,p}\right)\cdot {\gamma }^{f,p}\cdot \left(1-E_{}^{f,p}\right)\cdot T_{\mathit{forest}}$$

(2)

$${\gamma }^{f,p}=\exp \left(-{\tau }^{f,p}\cdot {\cos }^{-1}\theta \right)$$

(3)

where $E$ is the emissivity of the underlying surface (unitless), $T_{underlying}$ is the effective underlying surface temperature (unit: K), and $T_{\mathit{forest}}$ is the forest vegetation temperature (unit: K). ω is the forest single-scattering albedo (unitless); γ is the effective forest transmissivity (unitless) —which can be represented by the effective forest vegetation optical depth τ, as shown in Equation (3) — and $\theta$ is the incidence viewing angle toward nadir of the radiometer (unit: Rad).

The parameter $E$ is linked to the dielectric properties of the underlying surface, which are significantly influenced by properties of the ensemble medium, such as soil moisture and various features of snow medium covering the ground, including snow depth, volume fraction, wetness, snow cover fraction [36,37,38], etc. As $E$ is generally unknown, the parameters τ and ω are difficult to retrieve directly from the above equations [15,34]. Nevertheless, the parametric solution for τ (or γ) is achievable under some reasonable assumptions. Due to the relatively small magnitude of ω, we assumed ω to be negligible for this study according to literature practices [39,40,41]. In addition, most previous studies [11,19,42,43,44,45] assumed that $T_{air}\approx T_{\mathit{forest}}\approx T_{Gr}$, thereby Equation (2) yields to a simplified relationship:

$$T_{B,ac}^{f,p}=\left[1-{\left({\gamma }^{f,p}\right)}^2\cdot \left(1-E_{}^{f,p}\right)\right]\cdot T_{forest}$$

(4)

Further, the solution of γ is:

$${\gamma }^{f,p}=\sqrt{\frac{T_{forest}-T_{B,ac}^{f,p}}{\left(1-E_{}^{f,p}\right)\cdot T_{forest}}}$$

(5)

In the following, we solved Equation (5) to derive γ at the satellite pixel scale.

3.2. Pixel-Wise Forest Microwave Transmissivity Modeling

Owing to the large heterogeneity in land surface types within the coverage range of satellite pixels, the retrieved values of forest transmissivity are considered as the “effective transmissivity” at the pixel scale, rather than that at the forest plot scale, which also implies that the outputs of the Pixel-wise forest microwave transmissivity model (hereafter referred to as the “Pixel-wise γ Model”) are not dependent on the forest cover fraction parameter. To build the Pixel-wise γ Model, the following assumptions are given based on the foundation of previous research, in addition to those mentioned in Section 3.1:

(1)

Snow characteristics, effective permafrost emissivity and other environmental parameters are considered to be approximately the same for forest and non-forest areas within the same pixel and are set to equal values in the modeling [43];

(2)

The value of Pixel-wise γ remains almost constant during each winter season (September to March of the following year). This assumption is reasonable since previous studies [40,43,46] usually assumed that the winter γ at 19 and 37 GHz remains fixed over a more extended time period (generally the period of growth cycles in forest growing stock or aboveground biomass).

A flowchart of the modeling Pixel-wise γ for all four frequency and polarization combinations (19 GHz and 37 GHz, H and V polarizations) in freezing conditions is presented in Figure 2.

Three primary steps are included in the retrieval scheme.

Step 1: Determine the dates of a subset of stable dry snow conditions

Considering the feature matching accuracy error of land cover type data products and the improvement of computational efficiency, we took the lead in filtering out pixels with <10% forest cover (i.e., ff_forest <10%); meanwhile, to avoid the retrieval error caused by water bodies and wetlands, we excluded pixels with a total coverage of more than 40% water bodies (i.e., ff_water >40%).

From Equation (5), it can be deduced that the γ retrieval necessarily requires the determination of three parameters, $T_{B,ac}^{f,p}$, $E_{}^{f,p}$ and $T_{forest}$. The parameters $T_{B,ac}^{f,p}$ and $T_{forest}$ were obtained from the time-series SSMIS TB data and the ERA5-Land air temperature products, respectively, while parameter $E_{}^{f,p}$ needs to be qualified by setting a threshold condition. As aforementioned, $E_{}^{f,p}$ is related to the underlying surface dielectric properties, and therefore the fluctuation range of $E_{}^{f,p}$ for a relatively stable soil/snow state should theoretically be relatively insignificant. Consequently, our primary effort is to determine a relatively stable period or date of the underlying surface state, with a view to providing a constraint on $E_{}^{f,p}$. Thus, we proposed a rule for determining the dates of soil freezing with dry snow cover by referring to the identification algorithm mentioned in the snow product production specification of AMSRE [47]:

$$255 K\le T_{B\_91V}\le 261 K\hspace{0.33em}\&\hspace{0.33em}{T}_{forest}\le 273.15 K,$$

(6)

where $T_{B\_91V}$ is the TB values (unit: K) at 91 GHz V polarization channel, and $T_{forest}$ is the air temperature (unit: K). The range of values for $E_{}^{f,p}$ will be given in Step 2.

After searching pixel by pixel, and day by day, the datasets ({TB_19H, TB_19V, TB_37H, TB_37V, T_foerst}) of multiple dates that satisfy the above conditions were acquired. To circumvent the influence of microstructural properties of deep or wet snow (which usually occurs in late winter and early spring) on TB —for example, deep snow weakens the emissivity ($E_{}^{f,p}$) while wet snow enhances it, leaving it in a highly unstable state—the “autumn priority principle” as the screening condition was proposed. Specifically, for those that meet the above discriminatory conditions of Equation (6), we give priority to the data set ({TB_19H, TB_19V, TB_37H, TB_37V, T_foerst}) in autumn, and the eligible data set before December 31 of each year will be chosen and input to the model.

Since the influence of the atmosphere on the measured PMW TB is particularly pronounced when observations are made at high frequency microwaves (i.e., greater than 10 GHz) [48], the atmospheric decoupling of TBs at 19 GHz and 37 GHz is required. Following the simplified approach [44,49], the corrected TB ($T_{B,noatm}^{\uparrow }$) from the satellite measurements can be derived:

$$T_{B,noatm}^{\uparrow }=\frac{T_{B,sat}^{}-T_{B,atm\uparrow }}{t_{atm}}$$

(7)

where $T_{B,noatm}^{\uparrow }$ (unit: K) corresponds to the parameter $T_{B,ac}^{f,p}$ in Equation (5), $T_{B,sat}^{}$ represents the TB values (unit: K) acquired from SSMIS, $T_{B,atm\uparrow }$ is the upward atmospheric radiations (unit: K) and $t_{atm}$ is the atmospheric transmissivity (unitless). The parameters $T_{B,atm\uparrow }$ and $t_{atm}$ can be obtained from the atmospheric statistical model [50,51], which will be given in Section 3.3.2.

Step 2: Determine the effective emissivity of dry snow-covered frozen ground

In step 1, we have screened out the state period/date that satisfies the condition (Equation (6)) that the soil is frozen and covered with dry snow, pixel by pixel. Next, we need to determine the range of values for parameter $E_{}^{f,p}$ of the underlying surface in this state. Therefore, it needs to be particularly explicit that $E_{}^{f,p}$ here refers to the “effective emissivity”, meaning the underlying surface containing both snow medium and frozen soil medium is regarded as an integrated medium. Previous studies on the emissivity of snow-free frozen ground or snow surface have been conducted to show that: (1) Based on radiometric experiments in the farmland area of Northeast China, Wu et al. (2017) [52] measured the emissivity of snow-free frozen ground of 19 GHz and 37 GHz horizontal polarization to be 0.92 and 0.93, respectively; (2) Mätzler (1994) [53] pointed out that the polarization difference of permafrost emissivity is very small, and the polarization difference of 19 GHz and 37 GHz is 0.0068 and 0.0009, respectively. (3) Snow emissivity can vary from 0.4 to 1.0 due to complicated snow properties [53,54,55]; (4) Yan and Weng (2003) [56] retrieved the time series of globally mean land emissivity at all of seven SSM/I channels, during 1993–2002, over snow surface, giving the monthly mean emissivity with the value of 0.81, 0.91, 0.75 and 0.84 for 19 GHz (H), 19 GHz (V), 37 GHz (H) and 37GHz (V), respectively. Since the difference in emissivity polarization at the pixel level is not completely clear, this study directly selects the average value at dual polarization of each frequency band as the upper and lower limits of the threshold condition to participate in the model run. In addition, taking into account the complexity of the radiation characteristics of the forest floor, and referring to the results of preliminary trial calculation of the model, we set $E_{}^{f,p}$ as (increment is set to 0.01):

$$E_{\_19 GHz\_(H\hspace{0.33em}and\hspace{0.33em}V)}=0.83: 0.01: 0.93$$

(8)

$$E_{\_37 GHz\_(H\hspace{0.33em}and\hspace{0.33em}V)}=0.76: 0.01: 0.86$$

(9)

So far, three input parameters ($T_{B,ac}^{f,p}$, $T_{forest}$ and $E_{}^{f,p}$) of the Pixel-wise γ Model have been determined. Following this, substitute Equations (8) and (9) and the filtered data set ({TB_19H, TB_19V, TB_37H, TB_37V, T_forest}) from Step 1 into Equation (5) for iterative operations.

Step 3: Iterative optimization

By processing the first two steps, a matrix of γ retrieval values was generated pixel by pixel. Using hypothesis (2) in Section 3.2 that the Pixel-wise γ retrieved value remains almost unchanged in the winter snow period of each year, the optimization criterion, namely, the standard deviation (Std) of the changing sequence γ, should be minimal. In addition, it is necessary to add the restriction: 0 < γ < 1, to ensure the validity of γ in the iterative process. Ultimately, a multi-day γ series with minimum Std was picked in the matrix and averaged as the model output, i.e., the γ retrieval value.

3.3. Evaluation Approaches

To validate the output of the Pixel-wise γ Model at the satellite pixel scale (~25 km × 25 km), we proposed two validation schemes: one was to directly compare the “model output” with the “true value” using a limited number of measured γ data, and the other was to retrieve snow parameters (including SD and SWE) by embedding the Pixel-wise γ Model into the radiative transfer model system in order to complete indirect validation. The commonly used root mean-square error (RMSE), Bias and correlation coefficient (r) were used as evaluation metrics.

3.3.1. Direct Validation with Measured Data of Transmissivity

Due to spatial scale effects, the difficulty of verifying the retrieved values of γ at the pixel scale using the measured sample point data is undoubtable. Nevertheless, the measured data can visualize the reliability of the retrieved dataset in some ways. The forest microwave transmissivity measurement dataset contained 14 sample sites cumulatively, corresponding to being encapsulated in 8 SSMIS pixels, i.e., some pixels cover 2–3 sample sites. For the calculation of γ in multiple-site pixels, the average value was used. The reliability of the model output dataset is evaluated by directly comparing the difference between the measured and retrieved values.

3.3.2. Indirect Validation with Retrieval of SD and SWE

It is known that the removal of forest effects in the manufacturing framework of the Globsnow SD/SWE product [5,6,57] is based on the τ–ω model, with integrating the forest γ model into the overall RT modeling system. Following the modeling idea of Globsnow, we note that the Pixel-wise γ Model can also be easily incorporated into a snow emission model, such as the HUT model [58], to account for microwave attenuation in a forest canopy. Thus, the indirect validation can be achieved by replacing the parameter γ in the module of its RT modeling system and comparing the SD/SWE retrieval accuracy. Both direct forward modeling of SSMIS TB using in situ data for properties of snow cover and the inversion of forward models to retrieve SD/SWE were applied in this study. Note that, whilst we applied the forward model provided by the GlobSnow SWE product (v3.0), we used a different inversion algorithm (different cost function) to GlobSnow. The forward model information for the satellite measurement of TB simulations was summarized in

Table 1. Sub-module description and key parameter information of the forward model.

Model	Purpose	Key Parameters	Unit	Data Source or Value Taking
Atmospheric Statistical Model [46,47,50,51]	To calculate the atmospheric upwelling and downwelling emission and transmissivity	Air temperature	K	ERA5-Land
		Statistical coefficient	Unitless	Recalibration based on atmospheric physical model (MPM, [54]) and sounding data
Zeroth-order τ–ω Radiative Transfer Model [15]	To simulate the interaction between canopy and electromagnetic wave	Forest transmissivity	Unitless	Newly developed Pixel-wise forest transmissivity model in this study
		Single scattering albedo	Unitless	0
		Forest canopy temperature	K	ERA5-Land, assuming the forest canopy temperature and air temperature are the same
HUT Snow Emission Model [58]	To simulate snow emission of 19 GHz and 37 GHz	Snow temperature	K	China’s snow survey data
		Snow volumetric humidity	%	China’s snow survey data
		Snow density	kg/m3	China’s snow survey data
		SWE	cm	China’s snow survey data
		Snow grain size	mm	China’s snow survey data
Rough Bare Soil Reflectance Model [61]	To simulate the ground surface microwave emissivity	Soil surface roughness	mm	3
		Relative permittivity of permafrost	Unitless	$6-j$
		Q	Unitless	0.01-parameter calibrated
		H	Unitless	0.09-parameter calibrated
		NP	Unitless	0.92-parameter calibrated
		NV	Unitless	0.92-parameter calibrated

.

The atmospheric models used in this study included two types: a physical model with sounding data as input, i.e., the Millimeter Propagation Model (MPM) developed by Liebe (1989) [59], and the other model proposed by Pulliainen et al. (1993) [51], which was an empirical model based on statistical data. However, considering that actual atmospheric measurement data are difficult to obtain simultaneously over a large area, we modified the Atmospheric Statistical Model (ASM) using the MPM model for the practical needs of reducing the computational effort and simplicity of application. The specific expression of the ASM is as follows:

$$T_{B-\mathrm{atm},\uparrow (\downarrow )}={\alpha }_{\uparrow (\downarrow )}\cdot T_{air}\cdot \left(1-t_{atm}\right)$$

(10)

$$t_{atm}={\tau }_{atm}(1)+\zeta \cdot {\tau }_{atm}(2)$$

(11)

$${\tau }_{atm}(1)=[0.9851;0.9795;0.9211;0.8326;0.8624;0.6656]$$

(12)

$${\tau }_{atm}(2)=[0.0088;0.0275;0.2069;0.4642;0.2746;0.8163]$$

(13)

where $T_{air}$ is the surface air temperature (unit: K) which can be acquired from ERA5-Land product, ${\alpha }_{\uparrow (\downarrow )}$ is the approximate atmospheric profile factor (unitless) to determine the effective temperature [50]; $t_{atm}$ (unitless) is from statistical results [51,60]; Equations (12) and (13) are given statistical parameters, corresponding to frequencies of 6.9-, 10.65-, 18.7-, 23.8-, 37- and 89 GHz, respectively. $\zeta$ is a statistical coefficient that equals 0.0368 when applied in Finland and needs to be recalibrated when used in other regions. Using the MPM model output data set as a reference, we corrected the coefficient $\zeta$ from 0.0368 for Finland to 0.2406 for Northeast China by least squares fitting.

Next, it is necessary to parameterize the HUT snow emission model and rough bare soil reflectance model separately, using the pixel observations in the snow survey dataset, which can be referred to [57], where the optimized parameters of the soil model have been listed in Table 1. As for the HUT model, the focus of optimization is on its extinction coefficient, i.e.,

$${\kappa }_e=8\cdot f\cdot {\left[1.5\cdot \left(1-{\mathrm{e}}^{-1.5\cdot D_m}\right)\right]}^{3.7}$$

(14)

where f is the frequency (unit: GHz), and $D_m$ is the measured snow grain size (unit: mm).

After completing the initialization of the forward model parameters, an iterative process using an increment of 0.1 cm on SD (range: 0–40 cm) and of 0.1 mm on SWE (range: 0–80 mm) was used to retrieve the SD/SWE values, respectively, that minimized the RMSE between the measured (SSMIS) and simulated (simu) TB for incidence angle (θ) of 53° such that:

$$Minimize\left\{\frac{1}{2{\sigma }^2}\left\{\begin{array}{l}{\left[(T_{B,simu}^{19V}-T_{B,simu}^{37V})-(T_{B,SSMIS}^{19V}-T_{B,SSMIS}^{37V})\right]}^2+\\ {\left[(T_{B,simu}^{19H}-T_{B,simu}^{37H})-(T_{B,SSMIS}^{19H}-T_{B,SSMIS}^{37H})\right]}^2+\\ {\left[(T_{B,simu}^{19V}-T_{B,simu}^{19H})-(T_{B,SSMIS}^{19V}-T_{B,SSMIS}^{19H})\right]}^2+\\ {\left[(T_{B,simu}^{37V}-T_{B,simu}^{37H})-(T_{B,SSMIS}^{37V}-T_{B,SSMIS}^{37H})\right]}^2\end{array}\right\}\right\}$$

(15)

where $\sigma$ is the standard deviation of SSMIS measurements, which is assumed to be 2 K, according to the sensitivity characteristics of the SSM/I instrument [62].

Further, we chose the γ retrieval algorithm provided by the GlobSnow 3.0 product (hereafter referred to as the GlobSnow algorithm) as the object of comparison, with all other sub-module input parameters held constant, and only the parameter γ was retrieved by the Pixel-wise γ Model (hereafter referred to as the Pixel-wise algorithm) replaced in the forward model. The dual-frequency dual-polarization γ retrieval formula provided by the GlobSnow algorithm [46] is as follows:

$$tran{s}_{veg}=e^{-ke\cdot GSV}$$

(16)

where $ke$ is the forest vegetation extinction coefficient, which is 0.01, 0.007, 0.012 and 0.011 at 19 GHz and 37 GHz for H and V polarization, respectively; and GSV is the forest stem volume (unit: m³ ha⁻¹), which can be obtained from ESA BIOMASAR [63]. Finally, the Pixel-wise algorithm and the GlobSnow algorithm were applied to solve SD and SWE in reversion, respectively.

4. Results

4.1. Validation Results

4.1.1. Accuracy Validation Based on Measured Transmissivity

The results of the comparison between the retrieved and measured values of γ are shown in

Table 2. Comparison results of measured transmissivity and transmissivity estimated by Pixel-wise transmissivity model.

Point	Row_Id	Col_ Id	γ19H_Obs	γ19V_ Obs	γ37H_Obs	γ37V _Obs	γ19H_ Est	γ19V_ Est	γ37H_ Est	γ37V_ Est
Mohe01	5	31	0.63	0.57	0.60	0.58	0.62	0.61	0.60	0.59
Huma01	10	37	0.68	0.61	0.58	0.59	0.65	0.61	0.54	0.57
Huma02			0.68	0.61	0.58	0.59	0.65	0.61	0.54	0.57
Huma03	10	38	0.69	0.70	0.59	0.60	0.65	0.60	0.54	0.52
Huma04			0.69	0.70	0.59	0.60	0.65	0.60	0.54	0.52
Huma05			0.69	0.70	0.59	0.60	0.65	0.60	0.54	0.52
Xunke01	21	55	0.66	0.60	0.57	0.60	0.71	0.62	0.64	0.55
Xunke02			0.66	0.60	0.57	0.60	0.71	0.62	0.64	0.55
Xunke03			0.66	0.60	0.57	0.60	0.71	0.62	0.64	0.55
Xunke04	21	56	0.58	0.50	0.53	0.55	0.69	0.64	0.59	0.56
Yichun01	25	56	0.54	0.28	0.47	0.49	0.63	0.66	0.51	0.47
Yichun02			0.54	0.28	0.47	0.49	0.63	0.66	0.51	0.47
Yichun03	24	57	0.36	0.20	0.32	0.34	0.65	0.62	0.53	0.49
Yichun04	24	58	0.34	0.64	0.54	0.56	0.65	0.66	0.53	0.50

. The examples of γ presented in bold refer to the data points in which the difference between the estimated γ and the measured value with the Pixel-wise γ Model is beyond ±0.1, demonstrating that the estimation accuracy is relatively poor. It can be found that the γ at the individual frequencies/polarizations of the pixels where the sampling points Yichun02, Yichun04, Xunke04, and Yichun01 are located, and the γ at the full frequencies/polarizations of the pixels where the sampling point Yichun03 is located are overestimated, and the estimation accuracy is extremely low. One of the reasons for this phenomenon is that the overall forest density or accumulation on the SSMIS footprint is not as high as the actual measured sample sites, especially for the sample sites in Yichun, where the forest type is evergreen coniferous forest, the γ measurements are relatively low, mostly in the range of 0.2–0.36, which is still far from the saturation value of transmissivity at pixel scale (in the range of 0.5–0.6). This can be explained by the scale effect, as the dominant forest type in the pixel where the sample site located is not evergreen coniferous forest, but deciduous broadleaf forest and mixed forest. In general, the validation accuracy is high for the Huma region, represented by deciduous coniferous forest, and the Xunke region, represented by deciduous broadleaf forest, with an average error of only 0.01, but the validation effect is poor for the Yichun region, represented by evergreen coniferous forest (caused by the mismatch of dominant forest types).

4.1.2. Retrieval of SD and SWE

The SD and SWE results obtained by inversion with the Pixel-wise algorithm and GlobSnow algorithm, respectively, are shown in Figure 3, and the evaluation metrics are also superimposed in Figure 3.

From the graphical results in Figure 3, all evaluation metrics show that the Pixel-wise algorithm has significantly improved accuracy compared with the GlobSnow algorithm. From the index RMSE, the accuracy of SD and SWE inversion is improved by 6.8 cm and 7.75 mm, respectively, with an effective accuracy improvement of 41% and 26%, respectively; from the index Bias, the accuracy of SD and SWE inversion is improved by 84% and 95%, respectively; in addition, the correlation coefficient r is also improved. From the histogram of angular difference (difference = measured value − estimated value) in Figure 3, it can be found that the inversion errors of the Pixel-wise algorithm are mostly concentrated near the 0 point, while the errors of the GlobSnow algorithm are significantly right-skewed, i.e., the measured value is much larger than the estimated value, and the underestimation is severe. It is also found that the SD and SWE values of the inversion performance of the Pixel-wise algorithm have almost no boundary values, so it is more usable than the GlobSnow algorithm for the application in the forest area of Northeast China. In summary, it can be found that the Pixel-wise γ Model is more advantageous than the γ algorithm applied in the GlobSnow product for the inversion of snow parameters in the forest area of Northeast China. The adaptability of the proposed model is basically similar to GlobSnow as it does not depend on regional parameters and can be applied to regions besides Northeast China.

4.2. Spatio-Temporal Distribution Characteristics of Pixel-Wise Transmissivity

The spatial distribution of γ values for the H and V polarizations at 19 GHz and 37 GHz are shown in Figure 4, where the gray grid (with a special setting of 1.00) indicates that the forest pixel does not satisfy the Pixel-wise γ Model constraint, and the retrieved value is invalid. It can be found that the overall results show γ values at 19 GHz is higher than those of 37 GHz and the H polarization is higher than the V polarization. This is consistent with the theory that the penetration depth increases with decreasing frequency [16].

Taking the winter of 2014–2015 as an example, it can be seen that the values of the Pixel-wise γ range from 0.58–0.97 at the 19 GHz H polarization, while the frequency distribution histogram (Figure 5) shows that the values are concentrated in the grouping interval (i.e., the modal interval) of 0.66–0.68, and the value ranges of the frequency and polarization in other years are summarized in

Table 3. Summary of extreme value and mode interval distribution of transmissivity estimation results.

		2014–2015	2015–2016	2016–2017	2017–2018	2018–2019
19 GHz (H)	min	0.58	0.60	0.56	0.56	0.56
	max	0.97	0.97	0.97	0.97	0.96
	Bin_left *	0.66	0.66	0.66	0.66	0.66
	Bin_right *	0.68	0.68	0.68	0.68	0.68
19 GHz (V)	min	0.44	0.48	0.43	0.46	0.46
	max	0.97	0.97	0.94	0.96	0.95
	Bin_left *	0.62	0.62	0.62	0.62	0.62
	Bin_right *	0.64	0.64	0.64	0.64	0.64
37 GHz (H)	min	0.47	0.46	0.41	0.45	0.43
	max	0.95	0.95	0.94	0.97	0.92
	Bin_left *	0.54	0.54	0.54	0.54	0.54
	Bin_right *	0.56	0.56	0.56	0.56	0.56
37 GHz (V)	min	0.39	0.38	0.34	0.41	0.39
	max	0.88	0.86	0.81	0.84	0.76
	Bin_left *	0.50	0.49	0.50	0.52	0.51
	Bin_right *	0.51	0.50	0.51	0.53	0.52

* Bin_left and Bin_right denote the left and right boundaries of the group where the mode values are located, respectively.

. In addition, it can be seen in Figure 5 that the shape of the overall γ distribution is almost the same for each year, at the same frequency and polarization. The red dashed line in the figure gives the modal values of γ for each year, and it can be seen that the modal values of γ for V and H polarizations at 19 GHz and the H polarization at 37 GHz remain stable and constant for each year, and there are fluctuations of 0.01–0.02 in the modal values of γ for the V polarization at 37 GHz, with relatively small magnitude. As can be seen in Table 3, the center of the γ modal interval undergoes a shift of roughly 0.02–0.03 in the interannual variation, and the range of extremum variations are relatively stable, with fluctuations up to 0.12 under the overall spatial range.

One of the basic premises in developing this Pixel-wise γ Model is that the γ in each winter is considered to fluctuate very little during the snow season of the year, and can therefore be considered to be almost constant because the forest grows slowly in its natural state, and the water content of the stand branches is considered to remain relatively constant under winter freezing conditions. However, it is not clear whether the change in the trend of increasing or decreasing permeability occurs from year to year. Considering the natural environment, and by excluding activities such as deforestation, the standing stock increases year by year, especially for the satellite pixel of ~25 km × 25 km, the total forest volume should theoretically increase cumulatively, but whether or how much this cumulative volume affects the change of transmissivity values needs to be further investigated.

The method of combining Theil-Sen trend analysis and the Mann–Kendall significance test was adapted to analyze the change trend and significance of γ at the pixel scale. This method has become an important means to judge the trend of time series data and is widely used in meteorology and hydrology [64,65]. The Mann–Kendall significance test is described in detail by [66]. Figure 6 shows the Theil-Sen trend analysis and Mann–Kendall (M–K) significance test plots for five consecutive years from 2014 to 2019 under dual-frequency and dual polarization, with indicator 0 representing no change, and indicators −1 and +1 representing no significant decrease and no significant increase, respectively, and indicators −2 and +2 representing more significant decrease and more significant increase. From the results in Figure 6, it can be seen that the vast majority of Theil-Sen trend indicators for the performance of transmissivity values for each frequency and polarization throughout the Northeast region are between −1 and +1, i.e., they do not pass the significance test, and the trend of interannual variation is not obvious. To further verify this interannual variability without a significant trend, it can also be discriminated by calculating the coefficient of variance [67].

$$Cov=\frac{Std}{Avg}\times 100\%$$

(17)

where Cov is equal to the standard deviation divided by the averaged value. In accordance with the common statistical criterion of discrimination, i.e., abnormal fluctuations are considered to occur when Cov > 15% (or 0.15) [68], the percentage of pixels with Cov ≤ 0.15 is calculated in Figure 7.

From Figure 7, it can intuitively be found that the percentage of pixels with γ values at each frequency/polarization satisfying Cov ≤ 0.15 ranges from 99.12% to 99.85%, indicating that the interannual variation of γ in most of the pixels does not show abnormal fluctuations. Meanwhile, considering the systematic error of the Pixel-wise γ Model, Std and Cov are used as the uncertainty of the dataset, and quality control markers are applied pixel by pixel to further improve the dataset.

The spatio-temporal distribution analysis of pixel-wise transmissivity shows that the interannual changes of forest transmissivity during the period 2014–2019 are generally small. There are two reasons for this result. One is that the forest itself grows slowly, with little change between years. The other is that in recent years, the government of China has adopted the strictest protection system for the forest area and adopted the protection policy of closing the mountain for afforestation and reasonable thinning, which has minimized the impact of manual cutting and human interference.

5. Discussion

5.1. Uncertainty Factors Analysis of Pixel-Wise Forest Transmissivity Model

The determination of the dry snow-covered permafrost state is one of the important pre-processing steps in this Pixel-wise γ Model. Despite numerous attempts to distinguish the dry snow cover state of frozen soil, following (such as [47,69]) the output of the model showed that the number of pixels satisfying the discriminating condition was relatively small, and therefore, we could not solve the γ at the pixel level. Consequently, we gave the threshold condition (Equation (6)) following a large number of trial calculations. Those pixels that do not meet the condition may fail to estimate the transmissivity, but the proportion of such pixels in the study area is very small, only about 0.4–1% in this study.

The reason for choosing 91 GHz as a screening band is due to the strong sensitivity of snow emissivity at high frequencies (e.g., 91 GHz), i.e., more sensitive to the detection of changes in snow surface properties [56]. In the process of making further improvements to the existing model, the constraint condition of the stable dry snow-covered freezing period was the top priority of the study. Moreover, the accuracy of the product of temperature elements was also an important factor affecting the model output results. In recent years, studies by [23,70] have shown that the forest’s own protection mechanism [71] leads to the presence of unfrozen water in the canopy, despite the frozen environment, which affects the dielectric constant of the substance, and thus the γ increases with decreasing temperature. This study disregards the effect of unfrozen water on γ, which may introduce some uncertainty. In addition, in forest covered regions, the influence of land cover types should not be neglected in microwave snow cover studies [38,40,43,72]. The influence of forest and non-forest cover fractions in the pixel, closely related to the transmissivity model proposed in this research, will be discussed in the next section. Despite these uncertainties, the advantage of the Pixel-wise γ Model is that it does not require the support of ground measurements, and can output a set of γ data during each winter, based on the assumption that the physical properties of snow, frozen soil, and atmosphere are the same in forest and non-forest areas, and that the variation of γ in winter is negligible (which is a reasonable assumption), when only microwave TB data and air temperature data are required.

5.2. Influence of Forestcover Fraction

Another uncertainty concerning the model is that one of the assumptions of the newly developed Pixel-wise γ Model is that the effective emissivity of snow and underlying frozen soil is considered to be the same in forested and non-forested areas, which may introduce some errors into the model system itself. This will lead to different errors in the different forest cover pixels. Therefore, the results of the inversion of snow parameters for different forest covers were carried out, and the results are shown in

Table 4. Summary of pixel SD verification results in each forest coverage rate interval.

Interval	Sample Size	RMSE (cm)
		Pixel-Wise Model	GlobSnow Model
(0.1, 0.2]	10	9.1	12.0
(0.2, 0.3]	7	9.3	8.9
(0.3, 0.4]	6	10.4	9.8
(0.4, 0.5]	3	23.0	25.8
(0.5, 0.6]	1	15.9	16.2
(0.6, 0.7]	0	-	-
(0.7, 0.8]	6	10.3	13.8
(0.8, 0.9]	3	9.4	8.2
(0.9, 1]	51	8.4	18.7

. It can be seen that the inversion results for SD are promising for validation pixels with forest cover greater than 70%, especially those at 90%. Generally speaking, the estimation uncertainty associated with the retrieval model parameterization tends to increase with forest density [9]. However, for pixels with forest cover between 40% and 60%, retrieval results are slightly worse. However, it is worth noting that the sample size of the validation pixels is not uniformly distributed: about 58.6% of the cumulative 87 pixels have more than 90% forest cover, while the proportion of pixels in other forest cover ranges is relatively small, making it difficult to further quantify this model error caused by the difference in effective permafrost emissivity between forested and non-forested areas accurately. In summary, the relatively uneven distribution of forest cover for the validation pixels at large spatial scales is also one of the limitations of model validation. Nevertheless, the validation results of the pixel-wise transmissivity model in dense forests (>90% cover) show that the model performs well and in all forest cover intervals the proposed model exhibits advantages over the GlobSnow forest transmissivity model.

In summary, the Pixel-wise γ Model does not require the support of ground measurement data, but the selection of input parameters that determine the model output becomes particularly important. One of the key conditions set in the model is the identification of a critical environmental feature, namely the state period of dry snow-covered frozen ground, and subsequent attempts will be made to include other soil freeze and thaw models and precipitation identification models to further improve the accuracy of the model. In addition, the model does not consider the effect of dominant forest types or dominant tree species within satellite pixels, nor the effect of snowpack type [73], and future work will also focus on the inclusion of these two influences.

6. Conclusions

Using time-series TB data from satellite PMW observations and temperature data from the ERA5-Land product as model inputs, a new Pixel-wise γ Model was developed with the theoretical guidance of a zeroth-order vegetation radiative transfer model. The Pixel-wise γ Model produces winter microwave transmissivity data for five consecutive years for 19 GHz and 37 GHz horizontal and vertical polarization, covering the entire for-est region of Northeast China. The model validation results showed that the accuracy of SD and SWE inversions improved by 6.8 cm and 7.75 mm in terms of the index RMSE, effectively improving the accuracy by 41% and 26%, respectively; in terms of the index Bias, the accuracy of SD and SWE inversions enhanced by 84% and 95%, respectively. In addition, the correlation coefficients also improved. It is worth mentioning that the forest transmissivity model does not depend on forest cover parameters and other ground measurement parameters, which significantly improves the model’s application scope and easy applicability.

Developing transmissivity models at large regional scales is relatively difficult, especially for model validation, which requires more data support and careful experimental design. However, from the inversion results of snow parameters, the algorithms and products of the pixel-level transmissivity models developed in this research are still expected to further expand the spatial and temporal applications, and more detailed model refinement is needed in the future. This paper provides a product of dual-frequency and dual-polarized forest transmissivity datasets for five consecutive winters, of which the application area extends beyond Northeast China. It is expected to be used for global historical snow data reproduction and climate change analysis.