Inference of Soil Freezing Front Depth During the Freezing Period From the L-Band Passive Microwave Brightness Temperature

The freezing front depth (<italic>z</italic><sub>ff</sub>) of annual freeze–thaw cycles is critical for monitoring the dynamics of the cryosphere under climate change because <italic>z</italic><sub>ff</sub> is a sensitive indicator of the heat balance over the atmosphere-cryosphere interface. Meanwhile, although it is very promising for acquiring global soil moisture distribution, the <italic>L</italic>-band microwave remote sensing products over seasonally frozen grounds and permafrost is much less than in wet soil. This study develops an algorithm, i.e., the brightness temperature inferred freezing front (BT-FF) model, for retrieving the interannual <italic>z</italic><sub>ff</sub> with the diurnal amplitude variation of <italic>L</italic>-band brightness temperature (<italic>ΔT<sub>B</sub></italic>) during the freezing period. The new algorithm assumes first, the daily-scale solar radiation heating/cooling effect causes the daily surface thawing depth (<italic>z</italic><sub>tf</sub>) variation, which leads further to <italic>ΔT<sub>B</sub></italic>; second, Δ<italic>T<sub>B</sub></italic> can be captured by an <italic>L</italic>-band radiometer; third, <italic>z</italic><sub>tf</sub> and <italic>z</italic><sub>ff</sub> are negatively linear correlated and their relation can be quantified using the Stefan equation. In this study, the modeled soil temperature profiles from the land surface model (STEMMUS-FT, i.e., simultaneous transfer of energy, mass, and momentum in unsaturated soil with freeze and thaw) and <italic>T<sub>B</sub> </italic>observations from a tower-based <italic>L</italic>-band radiometer (ELBARA-III) at Maqu are used to validate the BT-FF model. It shows that, first, <italic>ΔT<sub>B</sub></italic> can be precisely estimated from <italic>z</italic><sub>tf</sub> during the daytime; second, the decreasing of <italic>z</italic><sub>tf</sub> is linearly related to the increase of <italic>z</italic><sub>ff</sub> with the Stefan equation; third, the accuracy of retrieved <italic>z</italic><sub>ff</sub> is about 5–25 cm; fourth, the proposed model is applicable during the freezing period. The study is expected to extend the application of <italic>L</italic>-band <italic>T<sub>B</sub></italic> data in cryosphere/meteorology and construct global freezing depth dataset in the future.


I. INTRODUCTION
A PPROXIMATELY 56% of topsoil in the northern hemisphere's landmass freezes at least for some time each year, including permafrost, which accounts for about 24% relatively [1]. Unlike permafrost in deeper soil depths, the soil freeze-thaw (FT) occurs at depths of centimeters to meters, drastically altering water and energy exchange between the land surface and the atmosphere [2], [3], [4]. In this term, the FT state transition refers to seasonal albedo changes that firmly shift the radiation balance over the land [5], [6], [7]. Besides albedo, the latent/sensible heat fluxes affected by FT have considerable changes at diurnal scales and impact the soil's water and energy states in the upper few centimeters [8], [9], [10], [11], [12], [13], [14]. Thus, FT state monitoring is sensitive and critical to weather prediction, climate forecasting, hydrology, and agriculture [15], [16]. However, due to high expense and poor representativeness, FT state monitoring is hard to achieve by regular measurements at in situ. On the other hand, the FT states are also critical to the forward modeling of microwave signals [17], [18], [19], [20], [21], which can be further applied to studies about land-atmosphere interactions with data assimilation [22], [23].
Since the launch of Skylab in the 1970s [24], passive remote sensing of the land surface condition, including FT state in microwave frequencies, has become a promising tool. Studies on the FT state detection started with shorter microwave bands, such as the Nimbus-7 Scanning Multichannel Microwave Radiometer SMMR [25], the Special Sensor Microwave/Imager SSM/I [26], and the Advanced Microwave Scanning Radiometer for EOS AMSR-E [27].
This study analyzes ΔT B observed by the ground-based Lband microwave radiometer ELBARA-III (the ETH L-Band radiometer) at Maqu in the Tibetan Plateau [16]. Then, it develops a brightness temperature inferred freezing front (BT-FF) model that can retrieve the freezing front depth (z ff ) from ΔT B . z ff is defined as the annual FT dynamics front in the soil profile. The rest of this article is organized as follows. Section II-A introduces the tower-based L-band radiometer and the Maqu Experiment; Section II-B describes the soil temperature simulation results from the land surface model and analyzes the impact of the FT transitions occurring on the soil surface or in deeper layers; Section II-C derives the relation between ΔT B and z ff from microwave transfer model [47] and Stefan's Equation [48]. Section III validates the BT-FF model with T B observations and the modeled soil temperature data. Finally, Section IV concludes this article.

II. OBSERVATIONS AND METHODOLOGY
The data to validate the BT-FF model contains T B observation from ELBARA-III and modeled soil temperature data from a land surface model named simultaneous transfer of energy, mass, and momentum in unsaturated soil with FT (STEMMUS-FT).

A. T B Observations
The Maqu monitoring network was established in 2008 in the eastern part of the Tibetan plateau [49], [50], [51] [see Fig. 1(a)]. The Yellow River bounds the network at its eastern and northern brinks, and its landscape is covered by meadow interspersed by a few trees or bushes. Its elevation is about 3300 m, and the rainy summer period caused by the East Asia Summer Monsoon is relatively short (July/August). On an annual scale, the wintertime can last for more than half a year, from October to March. Since December 2015, the L-band (1.41 GHz) radiometer ELBARA-III installed at Maqu's central station observes the typical meadow vegetation with a 40 m × 25 m footprint at incidence angles from 40°to 70°varying at 5°s teps [see Fig. 1 In this study, we analyze the observations at 40°to stay consistent with the SMAP observations and focus on a two-year period from August 2016 to August 2018 [see Fig. 2(a)]. The period contains two complete annual FT cycles and the freezing and thawing periods, respectively.
Furtherly, we separate the observations into periods according to daily soil temperature at 2.5 cm (T 2.5 cm ) characterized by unfrozen soil (I, a minimum of T 2.5 cm > 0°C), freezing soil (II, from I to III), frozen soil (III, a maximum of T 2.5 cm <0°C), and thawing soil (IV, from III to I) periods [the vertical dashed lines in Fig. 2(a)]. A lot of studies done with these dataset has proved that 2.5 cm is the most effective layer sensed by the radiometer [23], [33]. Even though the sensing depth varies for the frozen soil, 2.5 cm can be taken as an average and representative depth in this study.
Generally, ΔT B during Period I is dominated by rainfall and evaporation via changing θ and T. T B has higher values in the early afternoon and lowest values in the morning, and T B varies by about 10 K daily, respectively. In contrast, ΔT B is minimal during Period III because the high microwave penetration depth returns T B characteristic for a relatively large soil column; here, ΔT B stays less than 5 K. For Period II [see Fig. 2  interpret. The large |ΔT B | cannot be fully explained by either daily soil temperature or penetration depth variations. As shown in Fig. 2(b) and (c), |ΔT B | can reach 20 K daily. Only Subdaily FT transitions can explain it, e.g., the radiative heating over the daytime thaws the upper frozen soil, which means unfrozen soil appears at the surface [39]. Besides, in Period IV, the ponded water can drastically reduce the soil surface emissivity and penetration depth, mixing the signal of an FT transition and a water/ice phase transition.

B. STEMMUS-FT and Soil Temperature Data
With a 1 cm resolution in the vertical direction achieved by STEMMUS-FT, we can capture the daily thawing process caused by the daytime heating effect. The minimum vertical resolution that can be done by measurement is 2.5 cm, which cannot satisfy the requirements of this study. Furthermore, the sensors using TDR methods are hard to be installed over the soil surface where the depth is shallower than 2.5 cm due to their cylindrical sensing area. STEMMUS-FT [52] is a land model coded in MATLAB with extra modules based on STEMMUS [53], enhancing its capacity in FT simulations. In STEMMUS-FT, frozen soil is the thermal equilibrium system of soil grains, liquid water, ice, water vapor, and dry air [54], [55], [56], [57], [58]. The dynamics of the individual soil components, including soil\ice\water\vaper, interacting with each other were explicitly considered and resolved. Soil water and heat transfer are tightly coupled and realized by solving the equations of water conservation and energy conservation. Additional details on the STEMMUS(-FT) can be found in [59], [60]. Fig. 3 shows STEMMUS-FT's T simulations at Maqu, identifying freezing soil depth (z ff ) and thawing soil depth (z tf ). z ff is defined as T at z ff , where T upper layer < 0°C< T lower layer and reversely for z tf , where T upper layer > 0°C > T lower layer . During Period II, the amplitude of z tf due to the daily FT cycles is large, with a maximum of about 10 cm from the surface to the deeper layer. In Fig. 3, the amplitude of z tf corresponds to the development of the z ff . For instance, the daily freezing and thawing cycles were damped in Period II, and the growth of z ff was decreased in Period III. In Period IV, z tf goes deeper rapidly until the end of z ff in about April.
Compared with Fig. 2, |ΔT B | changes visually correlate well with the daily z tf changes affected by T variations between 0 and 0.2 m in Fig. 3. When the surface soil freezes at night, Fig. 4. Illustration of the freezing process assumptions in the BT-FF model. z ff is the freezing front, reflecting the annual FT state cycle, and z tf is the daily thawing front caused, e.g., by radiative heating during daytime. Symbol m is the total lasting days of the annual freezing process, i.e., the last day of period II. the deeper soil layers-frozen or thawed-will contribute to T B , but with a weaker attenuating effect from the frozen soil above than at the daytime. In Figs. 2 and 3, T B does not change much at night, even when the soil profile shows a frozen-moist (downwards, hereafter) structure. In contrast, daytime thawing of the surface soil significantly reduces the penetration depth, leading to a moist-frozen-moist structure in the soil profile. Especially during Periods II and IV, this moist-frozen-moist structure appears shortly during the daytime and completely disappears at night. Besides the FT state vertical structures, another reason that may cause |ΔT B | is open water due to melted snow or ice. It is noted that during Period IV, the frozen soil layer below the surface impedes draining and increases the residual water at the surface. Open surface water cannot appear in period frozen since snow is much more common than rainfall.

C. Derivation of the BT-FF Model
The proposed BT-FF model covers both the annual FT cycle regarding z ff and the daily FT processes, leading to ΔT B via z tf . The BT-FF model is based on the 0th-order microwave transfer model (see Section II-C-1) and relates ΔT B between two observation times within 24 h to the respective change of z tf (and z tf (n) for z tf on the nth day after the beginning of the yearly FT cycles). We estimate the temporal evolution of z ff from z tf /z tf (n)) via Stefan's equation (see Section II-C-2). To clarify, we refer z tf /z tf (n) to daily FT cycles while z ff refers to annual FT dynamics, e.g., during Period II (See symbols in the Appendix).
1) 0th-Order Microwave Transfer Model: The 0th-order incoherent model extensively used in passive microwave remote sensing of soil moisture is formulated as where T B is the observed brightness temperature in K, E is the soil slab's emissivity contributing to T B , and T eff is the so-called effective temperature of the soil slab. E depends on θ and other factors affecting the dielectric constant of a soil column. T eff relies on the dielectric constant as well. The penetration (e-folding) depth-as a measure of the soil slab depth mainly contributing to T B -depends on the dielectric constant profile and can reach a depth of meters down in extremely dry but also frozen soil at the L-band. For instance, Fig. 3 shows that in the first few days of Period II, the soil is frozen down to about 20-40 cm, and the soil below the frozen layer does not contribute significantly to the observed T B (e.g., a short period before z tf grows). We further assume that both T and dielectric constant are vertically constant in the frozen layer; then (1) can be rewritten as with T the frozen layer's temperature and E f the emissivity of the upper thawed soil layer. According to the red dots in Fig. 3, a thin upper soil layer (less than 20 cm) has thawed at 6 pm due to the daytime solar heating, giving rise to a thawing front at depth z tf . Considering the optical depth (τ tf ) [61], [62] down to z tf , (2) can be formulated as where E t is the emissivity of the upper thawed soil layer. By e −τ tf , the opacity stands for reflectivity/scattering of the internal/upper soil layer. We further assume that the temperature of the icewater mixture is always close to 0°C, which is the normal case during freezing and thawing, then ΔT B between morning and late afternoon is τ tf can be transferred to the corresponding z tf [62], [63] is given by with ε t and ε t the real and imaginary part of the dielectric constant of the thawed soil calculated from [64]. With a = 2πε t the penetration or e-folding depth in the thawed soil, we get which relates ΔT B to z tf . Thus, we can infer z tf with given E t and E f . z tf is quantified by the simulation result from STEMMUS-FT.

2) Stefan's Equation:
Both ΔT B and z tf are caused by surface thawing and refreezing (see Section II-C-1), and they vary daily. The z ff propagation (see Section II-C-2) evolves at an annual scale. We use Stefan's equation [48], [65], [66] to model both the propagation of z tf and z ff ( here, both fronts are indicated by z i ), based on the soil skin temperature time series via where Subscript j stands for either thawed (Subscript t) or frozen soil (Subscript f), c j = t 1 |T s (t)|dt the temporal integral of the soil skin temperature (T s ) in°C during freezing (T s < 0°C) or thawing (T s > 0°C) between the Time t 1 and t 2 , with t 1 /t 2 being the beginning/ending of either the annual FT cycle or the daily FT cycle. For a specific day in Period II, T s might increase above 0°C, leading to the surface soil ice thawing over time Δt = t 2 −t 1 . Accordingly, freezing (T s < 0°C) prevails during t d − Δt with t d the seconds for a day ( = 86 400 s). Equation (7) covers annual and diurnal FT cycles that depend on the integral time interval. If the integration stands for the annual cycle, we can obtain z ff by (7); otherwise, z tf is calculated by considering the daytime radiation heating effect over the frozen soil. Thus, we get z ff and z tf on day n after the start of Period II, z ff (n) and z tf (n), respectively, following (7) as where T st (n) is average T s at daily thawing time on Day n, Δt(n) the duration of the thawing on nth day, T sf (i) the average soil skin temperature during daily freezing lasting for t d − Δt(i) seconds on Day i, T sf (n) is the average skin temperature during freezing from day 1 to day n, and the overall freezing duration from day 1 to day n. z tf (n) and z ff (n) on day n are related via with c t c f assumed constant during the FT cycles since k i and θ are assumed to be constant. The average soil skin temperature on any day n (T s (n)) and its average until and including day n is given by Equation (10) means T s (n) can be separated into freezing (i.e., T sf (n)(t d −Δt(n)) t d ) and thawing process (i.e., T st (n)Δt(n) t d ), which can be used in (7). For instance, if Δt(n) = 1 2 t d and T sf (n) = −T st (n) on the first day (n = 1), then T s1 = T sf (n) + T st (n) = 0°C.
The T propagation with time and depth can be described as a sine wave by using a Fourier-serie [67], [68] as Authorized licensed use limited to the terms of the applicable license agreement with IEEE. Restrictions apply.
where z a /z b is the vertical coordinate positive downward, t is the time; T is the temporal average of the soil temperature, and A is half of the difference between the maximum and minimum at z b ; ω = 2π t d = 2π 86 400 s is the angular velocity of the Earth's rotation; φ is the initial phases of the soil temperature at depth z b ; and 1/ ω/2k i is the damping depth of the diurnal temperature wave. Thus, z a -z b = 0 and z a = 0. The freezing period starts from daily T s = 0°C, so we get T = T s = 0°C. We assume T s time-series whose maximum/minimum is T s (1)/−T s (1) on the day n = 1, then A = T s (1). In this case, a surface soil temperature time series without trend can be simplified as To further simply (11), the freezing period ends when the daily maximum T s = 0. Equation (12) is expented with a linear trend to mimic the daily T s variation during Period II as shown in [69] and [70] where m is the lasting days of Period II (i.e., the total number of days starting from daily T s = 0 to the ending at daily maximum T s = 0), and t is the accumulated seconds after Day i = 1. Equation (13) represents a T s time-series whose maximum/minimum is T s (1)/−T s (1) on the day n = 1 and gets below 0°C on the day i = n. Then, (13) is written as (14) and (10) is written as T s (n), T sf (n), T st (n) are defined by T x = 1 Δt T s (t)dt, where x can be i/sf/st or others depending on the averaging time interval for T. Thus, in the BT-FF model, T s (n) decreases from 0°C (i By considering T s (1), (9) can be rewritten as Since T s (1) is a constant, Fig. 5(a) focuses on the [sin(ωt) − 2 t mt d ] part to infer the ratio of z tf (n) z ff (n) in (16).  Fig. 5(b) shows the slope when m = 5 days, m = 25 days, and m = 100 days. Generally, the slope can be simplified as a constant that depends on m, i.e., how many days Period II can last.
Thus, the hypothesis that z ff (n) and z tf (n) in (16) and (6) are nearly linear related is solid. We now propose a linear regression function α(m) and β(m) can be identified by linear fitting of z ff (n) and z tf (n) from STEMMUS-FT simulation results [see Fig. 5(b)]. It is seen that the error of linear approximation is larger when the vertical coordinate is smaller. The relative error is about 10% for all m values. In other words, the shorter Period II for a given site, the smaller α(m) (α(m) < 0). Equation (17) is adopted in the following to infer z ff (n) from z tf (n). By using (17), the linear curve must fit for the first day in Period II when the L-band can detect the daily FT cycles, which means z ff (1) ≈ z tf (1), i.e., For m days of the L-band T B detectable in Period II, i.e., the daily maximum T s < 0°C, we get z tf (m) = 0. Then By considering (18) and (19), we get 3) BT-FF Model: With the linear relation given by (17) and its appropriateness demonstrated in Fig. 5(b), we can substitute (6) into (17), we get the BT-FF model as Equation (21) relates the ΔT B with the annual frozen front z ff (n). As mentioned above, we assumed the soil moisture content of thawed soil in the annual freezing process does not change with depth and time, which leads to b t = Constant. Thus, Parameters a, b t , are constant in the BT-FF model.

III. RESULTS
Equation (6) explains the relationship between z tf and ΔT B in theory. The parameters, i.e., a and b t , are defined explicitly. Thus, a and b t can be estimated with the dielectric constant. In detail, a dielectric constant model requires only clay, sand fraction, soil moisture (soil water content after thawing), soil temperature, and some optional inputs. The clay/sand fractions are constant for a fixed location, and T, if not crossing the freezing point of water, has less impact on the dielectric constant than θ and the FT states. According to [71], the clay fraction over the top is 9.85%, and the sand fraction is 26.95%. θ = 0.275 cm 3 /cm 3 is the STEMMUS simulated soil moisture content at the moment when Period II started in 2017. Fig. 6(a) and (b) illustrates the variation of a/b t along with θ using the Wu's model for frozen soil to caculate E f [64] and the Mironov's model for unfrozen soil' E t , ε t , and ε t [72].
Parameter a ranges from −20 K to 100 K while soil moisture increases from 0.05 cm 3 /cm 3 to 0.6 cm 3 /cm 3 [see Fig. 6(a)]. The negative a means the unfrozen arid soil has smaller emissivity than after frozen one. Parameter a also means that the maximum ΔT B can be caused by the FT transition so that the higher θ leads to larger ΔT B . It should be noted that Parameter a also depends on the polarization besides θ, but the difference between H-pol and V-pol is less than 10 K in Fig. 6(a). On the other hand, b t decreases with soil moisture exponentially. A certain z tf means the soil optical depth is reduced while freezing.
Equation (16) contains the variables c t /c f for the heat capacity when the soil bulk is unfrozen/frozen. For a site like Maqu, c t /c f could be considered constant once the soil moisture (or soil water content in either liquid or solid) is fixed. For Maqu, c t c f = 0.76 is adapted here. Thus, Parameter α(m) is determined by the lasting days of Period II (n). For Parameter β(m), it depends on n as well as the amplitude of the soil skin temperature time-series (i.e., T s (1)). Fig. 7(a) shows the theoretical values of α(m) along with the freezing days according to (20). The increasing α n means the z ff grows more slowly if the freezing process lasts longer. In Figs. 2(a)  While Figs. 6 and 7 give a full scale of a/b t /α(m) /β(m) estimates, we can also get these parameters from T B and STEMMUS-FT's simulations. In Fig. 8, a fitting curve in the form of (6) is used to get a relationship between ΔT B and z tf . ΔT B is mainly caused by the temperature difference between 6 A.M. and 6 P.M., and the emissivity changes of F/T status at the upper soil layer, i.e., it can be related to the thawing front. Thus, a = 68.26 is derived, which fits the range illustrated in Fig. 6(a). And b t = 0.06 is also obtained, which is equivalent to the value at the θ of about 0.35 m 3 /m 3 . The regression estimates of a/b t fit the BT-FF model. In the same way of estimating a/b t , Fig. 9 shows α(m) /β(m) by making a linear regression as shown in (19). It should be noted that the linear regression is made with only the maximum z tf in a day, not all z tf values we get from STEMMUS-FT. We  ELBARA-III's T B observations in Period II from 2016 to 2018. The z ff and z tf are compared with z ff and z tf inferred from T simulations in STEMMUS-FT in Fig. 10. Regarding z ff , the RMSE is 12 cm (5-25 cm for 95% confidence). For z tf , the RMSE is 7 cm.

IV. DISCUSSION
Although the BT-FF model has a clear mathematic expression, the uncertainty mostly comes from estimating parameters, i.e., a/b t from the field measurement and the assumption in derivation.
The BT-FF model is proposed to consider two major assumptions. The first one is that the frozen layer's T and dielectric constant are vertically constant. It is found that the evaporation and draining in the frozen soil are relatively weaker than in the unfrozen soil. Although the vapor diffusion affects the simulation of soil temperature, the effects are relatively small in degrees. The total amount of water removed from the frozen soil in winter is few [46], [73]. This grantee the assumption because the BT-FF model considers only period freezing when the phase change of water cannot be very active due to low daily air temperature.
Regarding soil temperature, the thermal inertia due to radiation heat is consumed as the heat of condensation and sublimation. T at the upper parts of the soil is around 0°C, while T at the bottom can always be considered constant in the daily time scale [74]. Thus, the second assumption is that the temperature of the ice-water mixture is always close to 0°C. This assumption is reasonable because the BT-FF model only applies period freezing when a rapid FT transfer happens daily. Regarding the depth interval of z tf (0-10 cm) and the huge thermal inertia of ice-water mixtures, the T gradient must be slight, and its daily average is almost 0°C. In contrast, the abovementioned two assumptions are not satisfied in period thawing. The annual freezing and thawing periods have different dynamics impacting ΔT B . For unfrozen soil, θ and T profiles change T B /ΔT B , especially in the top 10 cm. However, for frozen soil, ΔT B is dominated by the surface FT state transitions, while θ and T on the top 10 cm are not crucial for T B . These conditions simplify the microwave transfer model, leading to the BT-FF model and limiting its application to period freezing.
When it comes to modeling data, applying the BT-FF model at a global scale is not easy. It is challenging to infer z ff even from the land-surface models. For instance, SMOS needs the Hydrology Tiled ECMWF Scheme for Surface Exchanges over Land to provide land surface inputs for retrieval or forward simulations, while SMAP uses the modern-era retrospective analysis for research and applications, Version 2. Both land-surface models have too few layers for computing z ff , and thus, can hardly be used in evaluating the BT-FF model.

V. CONCLUSION
Previous works lack a theory that can bridge the gap between T B observation and the FT state of the soil profile during Period II when z ff propagates into deeper soil layers. This theoretical gap hampers the information we can retrieve from the L-band satellite, such as SMAP or SMOS. We propose the BT-FF model, which retrieves the freezing front depth (z ff ) from brightness temperature observations in the annual freezing period. The BT-FF model utilizes the Stefan Equation and zeroth-order microwave transfer model to achieve this goal. By taking T B observations from ELBARA-III and soil temperature simulations from STEMMUS-FT at Maqu, we conclude that ΔT B between 6 A.M. and 6 P.M. is an optimal pair in retrieving z ff . In the BT-FF model, θ and T profiles data are unnecessary because we assume θ and T in deeper layers vary little during Period II.
The concept of soil optical depth (τ ) is seldom used in analyzing T B at L-band. We think τ has advantages over the geometric depth frame (z) in analyzing T B signals under certain circumstances, such as the freezing/thawing period. The BT-FF model presented in this study is an application of τ in interpreting the T B signals at the L-band. The BT-FF model is expected to improve our understanding of the freezing/thawing process.