Evaluation of wet snow dielectric mixing models for L-band radiometric measurement of liquid water content in Greenland's percolation zone

The effective permittivity of wet snow and firn links the snow microphysics to its radiometric signature, making it essential for accurately estimating the liquid water amount (LWA) in the snowpack. Here, we compare ten commonly used microwave dielectric mixing models for estimating LWA in wet snow and firn using L-band radiometry. We specifically focus on the percolation zone of the Greenland Ice Sheet (GrIS), where the average volume fraction of liquid water is between 0 % and 6 %. We used L-band brightness temperature (TB) observations from the NASA Soil Moisture Active Passive (SMAP) mission in an inversion-based framework to estimate LWA, applying different dielectric mixing formulations in the forward simulation. We compared the effective permittivities of the mixing models over a range of conditions and evaluated their impact on the LWA retrieval. We also compared the LWA retrievals to the corresponding LWA from two state-of-the-art Surface Energy and Mass Balance (SEMB) models. Both SEMB models were forced with in situ measurements from automatic weather stations (AWS) of the Programme for Monitoring of the Greenland Ice Sheet (PROMICE) and Greenland Climate Network (GC-Net) located in the percolation zone of the GrIS and initialized with relevant in situ profiles of density, stratigraphy, and sub-surface temperature measurements. The results show that the mixing models produce substantially different real and imaginary parts of the dielectric constant, significantly impacting the LWA retrieved from the TB. The correspondence with the SEMB-derived LWA varied by model and site, with correlation coefficients ranging from 0.67 to 0.98 and RMSD values between 5.4 and 23.9 mm. Overall, the power law-based empirical models demonstrated better performance for 2023 melt season. The analysis supports informed selection of dielectric mixing models for improved LWA retrieval accuracy.

1 Introduction

Surface melting and consequent runoff/refreezing play an increasingly major role in the Greenland Ice Sheet (GrIS) surface mass balance (SMB) and its contribution to the global sea-level rise (Greene et al., 2024; Khan et al., 2022; Khan et al., 2015; Mouginot et al., 2019; Otosaka et al., 2023; Shepherd et al., 2020). The column-integrated amount of liquid water (LWA) at the surface and within layers of the surface snowpack is a key variable for understanding processes related to meltwater on the ice sheet, and is thus an important quantity for diagnostic study, modeling, and prediction. Currently, there is no direct means of measuring LWA in the ice sheet. In situ AWS provides surface meteorological observations for limited locations over the ice sheet, which are translated to LWA estimates using coupled surface energy balance and sub-surface hydrology and heat transfer models (Fausto et al., 2021; Samimi et al., 2021; Vandecrux et al., 2020). Regional climate models provide pan-ice sheet estimates of LWA (Fettweis et al., 2020), but uncertainty results from the significant differences in the configuration and physical process representations in these models (Thompson-Munson et al., 2023; Fettweis et al., 2020; Vandecrux et al., 2020; Verjans et al., 2019). Spaceborne microwave radiometers have also been used for large-scale mapping of polar ice sheet melt (Picard et al., 2022; Tedesco, 2007; Tedesco et al., 2007; Abdalati and Steffen, 1997; Mote and Anderson, 1995; Zwally and Fiegles, 1994). However, shallow penetration into the wet snow restricts the conventional high frequency radiometers (i.e., greater than 6 GHz) to providing only surface and near-surface binary melt status, and not the actual volumetric amount of liquid water in the snow/firn (Leduc-Leballeur et al., 2025; Colliander et al., 2022a, b, 2023; Mousavi et al., 2022).

The higher penetration of L-band radiometry offers a promising new tool for quantifying the total surface-to-subsurface LWA in the firn, in addition to providing the areal extent and duration of seasonal surface snow melt (Houtz et al., 2019, 2021; Mousavi et al., 2021; Schwank and Naderpour, 2018; Colliander et al., 2022a, b; Miller et al., 2020a, 2022a, b; Mousavi et al., 2022). Houtz et al. (2019, 2021) used L-band brightness temperature (TB) from the European Space Agency's Soil Moisture and Ocean Salinity (SMOS) mission for simultaneous estimation of snow liquid water content and density in the GrIS. They used the Microwave Emission Model of Layered Snowpacks Version 3 (MEMLS V3; Mätzler and Wiesmann, 2012) with L-band specific modifications (LS-MEMLS; Schwank et al., 2014) in an inversion-based retrieval framework. By default, MEMLS V3 uses the Mätzler (1996) and Mätzler and Wiesmann (2012) formulations for dielectric mixing of dry and wet snow, respectively. Naderpour et al. (2021) used the same algorithm to quantify LWA at the Swiss Camp location (70° N, 49° W) with close-range (CR) single-angle L-band microwave radiometer measurements. Mousavi et al. (2021) developed an L-band specific snow/firn radiative transfer model that uses the Mätzler (2006) and Ulaby et al. (2014, p. 140–145) dielectric mixing model for dry and wet snow, respectively, to estimate LWA. Hossan et al. (2025) used the same approach to quantify and validate the LWA with two surface energy balance models forced with in situ observations and reanalysis data products. Additionally, Moon et al. (2025) compared the Hossan et al. (2025) retrieval with LWA values derived from subsurface thermal measurements. The studies showed mixed correspondence of the L-band retrievals to the alternative LWA estimates.

The L-band TB responds to the real and imaginary parts of the firn dielectric constant, which increases markedly with volumetric liquid water content (snow wetness or the volume fraction of liquid water in the snow mixture, v_w hereafter, expressed as a percentage) in the firn (Picard et al., 2022; Samimi et al., 2021; and references therein). The measured dielectric constant is translated into LWA using a model between snow v_w and the dielectric constant. The formulation of the effective dielectric constant of the ice, air, and water mixture is key to accurately quantifying LWA. As it is independent of the radiometer measurement, it adds an uncertainty component to the LWA retrieval that is solely dependent on the accuracy of this dielectric mixing model. Picard et al. (2022) demonstrated large differences in commonly used wet snow dielectric mixing models for both the real and imaginary parts over a wide range of v_w.

In this manuscript, we assess the performance of ten commonly used microwave dielectric mixing models in quantifying the seasonal LWA using L-band (1.4 GHz) enhanced-resolution (rSIR) TB data products from the Soil Moisture Active Passive (SMAP) mission. For this, we confine our attention to the GrIS percolation zone where the average volume fraction liquid water inclusions in the snow/firn environment is less than about 6 % of the total volume (Colbeck, 1974; Coléou and Lesaffre, 1998), and L-band TB during melting is mainly dominated by absorption (increasing trend compared to frozen season; Hossan et al., 2025).

2 Methods

Snow and firn (the transitional snow that survives at least a summer season) is generally a three-phase porous dielectric mixture of air, ice, and liquid water, where dry snow can be considered as a special case having no liquid water. Here we briefly discuss the dielectric properties of snow and firn.

2.1 Dielectric Mixing Formulas

Dielectric mixing rules attempt to approximate the effective macroscopic dielectric constant (also known as effective permittivity) ε_eff of snow/firn mixture that relates the average electric flux density D to the incident or emitting electric field E and the average polarization P of the mixer (Sihvola, 1999; Jones and Friedman, 2000),

$$\begin{array}{}\text{(1)}& \mathbf{D}={\mathit{\epsilon }}_{\mathrm{eff}}\mathbf{E}={\mathit{\epsilon }}_{\mathrm{e}}\mathbf{E}+\mathbf{P}\end{array}$$

where ε_e is the complex dielectric constant of the background medium. The polarization term, P represents the number of dipole moments per unit volume in the mixture and is a function of inclusion geometry. ε_eff thus, depends on the individual dielectric constant of the constituent materials, their respective volume fractions, and sizes, shapes, and orientations. For a two phase mixture, a generalized mixing formula can be derived from the Maxwell Garnett (MG) mixing rule (Garnett, 1904),

$$\begin{array}{}\text{(2)}& {\mathit{\epsilon }}_{\mathrm{eff}}={\mathit{\epsilon }}_{\mathrm{e}}+\mathrm{3}v{\mathit{\epsilon }}_{\mathrm{e}}{\displaystyle \frac{{\mathit{\epsilon }}_{\mathrm{i}}-{\mathit{\epsilon }}_{\mathrm{e}}}{{\mathit{\epsilon }}_{\mathrm{i}}+\mathrm{2}{\mathit{\epsilon }}_{\mathrm{e}}-v({\mathit{\epsilon }}_{\mathrm{i}}-{\mathit{\epsilon }}_{\mathrm{e}})}}\end{array}$$

where ε_i is the complex dielectric constant of the inclusion (guest), and v is the volume fraction of the inclusion. However, the fundamental MG rule considers only dilute concentration of spherical inclusions (v≪1) (Jones and Friedman, 2000) and assumes homogeneous mixtures, ignoring the second-order effects due to the mutual interactions between the inclusions. When the inclusions are arbitrarily spread within the host material, the fields within the inclusions are a function of the mutual interactions of the inclusions (through their polarization field P). The interactions among the inclusions are dependent on their relative geometry and alignment, which are both taken into account by a parameter called the depolarization factor. For an ellipsoid, the depolarization factor along its u axis is given by (Jones and Friedman, 2000; Sihvola, 1999),

$$\begin{array}{}\text{(3)}& \begin{array}{rl}& N_u={\displaystyle \frac{abc}{\mathrm{2}}}\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}{\displaystyle \frac{\mathrm{d}s}{\left(s+u^{\mathrm{2}}\right)\sqrt{\left(s+a^{\mathrm{2}}\right)\left(s+b^{\mathrm{2}}\right)(s+c^{\mathrm{2}})}}}\\ & u=a,b,c\end{array}\end{array}$$

where $N_a+N_b+N_c=\mathrm{1}$. Aspect ratios of the axial dimensions (u) of the particle describe the shape of the particle, $a:b=\mathrm{1}$ defines the spherical inclusions whereas $a:b<\mathrm{1}$ and $a:b>\mathrm{1}$ describe oblate and prolate inclusions, respectively. The depolarization factor N_u is commonly included in the more general mixing formula described below.

The Bruggeman mixing rule (Ulaby et al., 2014; Sihvola, 1999) considers mixing phases to be symmetric to describe the effective permittivity of a mixture as an implicit function of unknown effective permittivity of the mixture, which for the case of randomly oriented ellipsoidal inclusions (Sihvola, 1999) reads as,

$$\begin{array}{}\text{(4)}& {\mathit{\epsilon }}_{\mathrm{eff}}={\mathit{\epsilon }}_{\mathrm{e}}+{\displaystyle \frac{v}{\mathrm{3}}}\left({\mathit{\epsilon }}_{\mathrm{i}}-{\mathit{\epsilon }}_{\mathrm{e}}\right)\sum _{u=a,b,c}\left[{\displaystyle \frac{{\mathit{\epsilon }}_{\mathrm{eff}}}{{\mathit{\epsilon }}_{\mathrm{eff}}+N_u({\mathit{\epsilon }}_{\mathrm{i}}-{\mathit{\epsilon }}_{\mathrm{eff}})}}\right]\end{array}$$

Polder and van Santen (1946) (PVS hereafter) derived similar implicit formulation for a two-phase mixture with randomly oriented ellipsoidal inclusions. de Loor (1968) extended the work as follows,

$$\begin{array}{}\text{(5)}& {\mathit{\epsilon }}_{\mathrm{eff}}={\mathit{\epsilon }}_{\mathrm{e}}+{\displaystyle \frac{v}{\mathrm{3}}}\left({\mathit{\epsilon }}_{\mathrm{i}}-{\mathit{\epsilon }}_{\mathrm{e}}\right)\sum _{u=a,b,c}\left[{\displaystyle \frac{\mathrm{1}}{\mathrm{1}+N_u\left(\frac{{\mathit{\epsilon }}_{\mathrm{i}}}{{\mathit{\epsilon }}^{\ast }}-\mathrm{1}\right)}}\right]\end{array}$$

where ε^∗ is the effective dielectric constant of the region surrounding inclusions. For v≤0.1, ${\mathit{\epsilon }}^{\ast }\approx {\mathit{\epsilon }}_{\mathrm{e}}$ and for higher value of f, ε^∗ approaches to ε_i (de Loor, 1968).

Another widely used mixing rule is the Coherent Potential (CP here after) formula (Tsang et al., 1985), which for the case of randomly oriented ellipsoidal inclusions is given by (Sihvola, 1999),

$$\begin{array}{}\text{(6)}& {\mathit{\epsilon }}_{\mathrm{eff}}={\mathit{\epsilon }}_{\mathrm{e}}+{\displaystyle \frac{v}{\mathrm{3}}}\left({\mathit{\epsilon }}_{\mathrm{i}}-{\mathit{\epsilon }}_{\mathrm{e}}\right)\sum _{u=a,b,c}\left[{\displaystyle \frac{\left(\mathrm{1}+N_u\right){\mathit{\epsilon }}_{\mathrm{eff}}-N_u{\mathit{\epsilon }}_{\mathrm{e}}}{{\mathit{\epsilon }}_{\mathrm{eff}}+N_u({\mathit{\epsilon }}_{\mathrm{i}}-{\mathit{\epsilon }}_{\mathrm{eff}})}}\right]\end{array}$$

For dilute inclusion (v≪1), all of these formulas (Eqs. 4–6) provide the same results as the MG mixing rule (Sihvola, 1999). However, as v increases the MG formula usually predicts ε_eff closer to ε_e (dielectric constant of host or background environment) which is lower than that estimated by both the PVS and CP formulations. This is because, as mentioned above, MG neglects the second-order effects due to the mutual interactions between neighbouring inclusions. The CP formula considers the effective medium instead of the background to find the local field and, therefore, estimates higher ε_eff compared to ε_e. The PVS formula, on the other hand, represents a balance between MG and CP, as it treats both the inclusions and the surrounding environment symmetrically, resulting in an equal influence from the permittivity of the two phases (Jones and Friedman, 2000; Sihvola, 1999). However, a computational difficulty of PVS and CP formulae is that they are both implicit in ε_eff.

Tinga et al. (1973) derived an explicit formula for ε_eff by considering a two-phase mixture composing of two randomly oriented confocal ellipsoids with an inner ellipsoid representing the inclusion and an outer ellipsoidal shell representing the host material.

$$\begin{array}{}\text{(7)}& \begin{array}{rl}{\mathit{\epsilon }}_{\mathrm{eff}}& ={\mathit{\epsilon }}_{\mathrm{e}}+{\displaystyle \frac{v}{\mathrm{3}}}\left({\mathit{\epsilon }}_{\mathrm{i}}-{\mathit{\epsilon }}_{\mathrm{e}}\right)\\ & \cdot \sum _{u=a,b,c}\left[{\displaystyle \frac{\mathrm{1}}{\mathrm{1}+(N_{u\mathrm{2}}-f{N}_{u\mathrm{1}})(\frac{{\mathit{\epsilon }}_{\mathrm{i}}}{{\mathit{\epsilon }}_{\mathrm{e}}}-\mathrm{1})}}\right]\end{array}\end{array}$$

where N_u2 and N_u1 are the depolarization factors of the inner and outer ellipsoids, respectively. It is noted that for spherical inclusions ($N_{u\mathrm{2}}=N_{u\mathrm{1}}=\frac{\mathrm{1}}{\mathrm{3}}$), Eq. (7) also reduces to the MG mixing rule, but any deviations from the spherical shape increase the ε_eff if ε_i>ε_e (and vice versa).

Another group of mixing formulas follows empirical power-law relations, where a certain power of effective permittivity of a multi-phase mixture relates to the linear combination of components raised to the same power and weighted by their respective volume fractions, v (Sihvola et al., 1985; Sihvola, 1999). These exponential models do not explicitly consider the microstructure shapes (i.e. through depolarization factor, N) but take into account the higher order mutual interactions through the power-law averaging. The general form of these models takes the form,

$$\begin{array}{}\text{(8)}& {\mathit{\epsilon }}_{\mathrm{eff}}^{\mathit{\beta }}=\sum _j{v}_j{\mathit{\epsilon }}_j^{\mathit{\beta }}\end{array}$$

where v_j and ε_j are the volume fraction and dielectric constant of the jth constituent, respectively, and $\sum _j{v}_j=\mathrm{1}$. The exponent β controls the degree of nonlinearity of the model (Sihvola et al., 1985), which is bounded by $\mathrm{0}<\mathit{\beta }\le \mathrm{1}$. The lower the value of β, the higher the influence of the background (dominant volume fraction).

The effective dielectric constant of a mixture is also calculated in the frequency domain using dispersion models. In case of inclusions with permanent electric dipole moments, like liquid water in wet snow, the Debye relaxation model is best suited (Hallikainen et al., 1986; Sihvola, 1999). Wet snow shows a distinct Debye relaxation spectrum in the microwave range (Ulaby et al., 2014, p. 140–145). The Debye-like semi-empirical models are of the form (Hallikainen et al., 1986),

$$\begin{array}{}\text{(9a)}& {\displaystyle }{{\mathit{\epsilon }}^{\prime }}_{\mathrm{eff}}=A+{\displaystyle \frac{B{v}_{\mathrm{w}}^x}{\mathrm{1}+(\frac{f}{f_{\mathrm{0}}}{)}^{\mathrm{2}}}}\text{(9b)}& {\displaystyle }{{\mathit{\epsilon }}^{\prime \prime }}_{\mathrm{eff}}={\displaystyle \frac{C\left(\frac{f}{f_{\mathrm{0}}}\right)v_{\mathrm{w}}^x}{\mathrm{1}+{\left(\frac{f}{f_{\mathrm{0}}}\right)}^{\mathrm{2}}}}\end{array}$$

where ${\mathit{\epsilon }}_{\mathrm{eff}}^{\prime }$ and ${\mathit{\epsilon }}_{\mathrm{eff}}^{\prime \prime }$ are the real and imaginary parts of the effective dielectric constant of the mixture and v_w is the volume fraction of liquid water in snow. f and f₀ are the operational and relaxation frequencies respectively and $A,B,C$, and x are constants that are determined empirically by fitting experimental data.

There are numerous models and formulas in the literature describing the dielectric behaviour of mixtures. Comprehensive reviews on the topic can be found in Sihvola (1999) and Mätzler (2006), and the references therein. Many of these formulas are special cases or modifications of the above basic mixing rules. Some others are empirical in nature. In the following section, we will briefly describe some specific wet snow mixing models that we evaluated in this study.

2.2 Dielectric Constant of Dry Snow

Dry snow is a two-phase mixture of ice and air. Since the real part of the dielectric constant of ice, ${\mathit{\epsilon }}_{\mathrm{i}}^{\prime }$, is independent of frequency and almost independent of temperature, it is assumed that the real part of the dielectric constant of dry snow, ${\mathit{\epsilon }}_{\mathrm{ds}}^{\prime }$ is also independent of both frequency and temperature (Hallikainen et al., 1986). ${\mathit{\epsilon }}_{\mathrm{ds}}^{\prime }$ is thus fully determined by the density of the dry snow (Denoth, 1989; Denoth et al., 1984; Tiuri et al., 1984). However, the imaginary part of ice dielectric constant, ${\mathit{\epsilon }}_{\mathrm{i}}^{\prime \prime }$, and thus the dry snow ${{\mathit{\epsilon }}^{\prime \prime }}_{\mathrm{ds}}$ are strongly sensitive to both frequency and temperature (Ulaby et al., 2014). With known dielectric constants of air and ice, the above dielectric mixing models, such as two phase PVS mixing rule, can be applied to find the effective dielectric constant of dry snow. Empirical formulations based on experimental data also provide good results (e.g., Mätzler, 2006). For the real part of dry snow permittivity, we follow the empirical relation presented in Mätzler (2006),

$$\begin{array}{}\text{(10)}& {{\mathit{\epsilon }}^{\prime }}_{\mathrm{ds}}=\left\{\begin{array}{ll}\mathrm{1}+\mathrm{1.4667}{v}_{\mathrm{i}}+\mathrm{1.435}{v}_{\mathrm{i}}^{\mathrm{3}}& \text{for }\mathrm{0}\le v_{\mathrm{i}}\le \mathrm{0.45}\\ {\left(\mathrm{1}+\mathrm{0.4759}{v}_{\mathrm{i}}\right)}^{\mathrm{3}}& \text{for }{v}_{\mathrm{i}}\ge \mathrm{0.45}\end{array}\right.\end{array}$$

where v_i is the volume fraction of ice in snow given by the ratio of snow and ice density.

For the imaginary part, ${{\mathit{\epsilon }}^{\prime \prime }}_{\mathrm{ds}}$, we follow the Hallikainen et al., (1986) formulation based on Tinga mixing model (Eq. 7),

$$\begin{array}{}\text{(11)}& {{\mathit{\epsilon }}^{\prime \prime }}_{\mathrm{ds}}={\displaystyle \frac{\mathrm{0.34}{v}_{\mathrm{i}}{{\mathit{\epsilon }}^{\prime \prime }}_{\mathrm{i}}}{{\left(\mathrm{1}-\mathrm{0.42}{v}_{\mathrm{i}}\right)}^{\mathrm{2}}}}\end{array}$$

where ${{\mathit{\epsilon }}^{\prime \prime }}_{\mathrm{i}}$ is the imaginary part of the dielectric constant of ice determined based on Mätzler (2006).

Various formulations are available for the effective dielectric constant of dry snow. Since our focus in this manuscript is wet snow mixing models, we apply the same dry snow model (Mätzler, 2006) for the mixture where the host is dry snow.

2.3 Dielectric Constant of Wet Snow

In this subsection, we discuss ten commonly used dielectric mixing models for estimating the complex dielectric constant of wet snow. Mätzler model (Mätzler and Wiesmann, 2012; used in MEMLS V3) uses the MG rule to compute the effective dielectric constant of wet snow as a two-phase mixture of liquid water inclusions in a dry snow background (host) with experimentally determined depolarization factors N_a=0.005, $N_b=N_c=\mathrm{0.4975}$ (representing a prolate spheroidal shape of the inclusion) from Hallikainen et al. (1986) and Matzler et al. (1984),

$$\begin{array}{}\text{(12a)}& {\displaystyle }{\mathit{\epsilon }}_{\mathrm{eff}}={\displaystyle \frac{\left(\mathrm{1}-v_{\mathrm{w}}\right){\mathit{\epsilon }}_{\mathrm{ds}}+v_{\mathrm{w}}{\mathit{\epsilon }}_{\mathrm{w}}K}{\left(\mathrm{1}-v_{\mathrm{w}}\right)+v_{\mathrm{w}}K}}\text{(12b)}& {\displaystyle }K={\displaystyle \frac{\mathrm{1}}{\mathrm{3}}}(K_{\mathrm{a}}+K_b+K_c)\text{(12c)}& {\displaystyle }{K}_u={\displaystyle \frac{{\mathit{\epsilon }}_{\mathrm{ds}}}{{\mathit{\epsilon }}_{\mathrm{ds}}+N_u({\mathit{\epsilon }}_{\mathrm{w}}-{\mathit{\epsilon }}_{\mathrm{ds}})}}u=a,b,c\end{array}$$

where ε_w is the dielectric constant of pure water which is determined from Liebe et al. (1991) relaxation formula.

Tinga model (Tinga et al., 1973) also considers wet snow as a three-phase mixture (air-snow-liquid water) where air is the background and water is a spherical shell surrounding another confocal shell (ice). The effective dielectric constant is then determined following Eq. (7). Hallikainen et al. (1986) first proposed a Debye-like model for wet snow with the following expression for the constants in Eq. (9).

$$\begin{array}{}\text{(13a)}& {\displaystyle }A=\mathrm{1}+\mathrm{1.83}{\mathit{\rho }}_{\mathrm{ds}}+\mathrm{0.02}{A}_{\mathrm{1}}{v}_{\mathrm{w}}^{\mathrm{1.015}}+B_{\mathrm{1}}\text{(13b)}& {\displaystyle }B=\mathrm{0.073}{A}_{\mathrm{1}}\text{(13c)}& {\displaystyle }C=\mathrm{0.073}{A}_{\mathrm{2}}\text{(13d)}& {\displaystyle }x=\mathrm{1.31}\text{(13e)}& {\displaystyle }{f}_{\mathrm{0}}=\mathrm{9.07}\mathrm{GHz}\end{array}$$

where for the original Debye-like model, $A_{\mathrm{1}}=A_{\mathrm{2}}=\mathrm{1}$, and B₁=0.

Hallikainen et al. (1986) derived the expressions below for the constants A₁, A₂, and B₁ in Eq. (13) as a function of frequency by fitting Eq. (13) to the field measurements of v_w with a range 0 %–12 %, a density range of 0.09–0.42 g cm⁻³, a temperature range of −15–0 °C, grain radius covering 0.5–1.5 mm at the frequency range of 3 to 37 GHz (see Hallikainen et al., 1986).

$$\begin{array}{}\text{(14a)}& {\displaystyle }{A}_{\mathrm{1}}=\mathrm{0.78}+\mathrm{0.03}f-\mathrm{0.58}\times {\mathrm{10}}^{-\mathrm{3}}{f}^{\mathrm{2}}\text{(14b)}& {\displaystyle }{A}_{\mathrm{2}}=\mathrm{0.97}-\mathrm{0.39}f\times {\mathrm{10}}^{-\mathrm{2}}+\mathrm{0.39}\times {\mathrm{10}}^{-\mathrm{3}}{f}^{\mathrm{2}}\text{(14c)}& {\displaystyle }{B}_{\mathrm{1}}=\mathrm{0.31}-\mathrm{0.05}f+\mathrm{0.87}\times {\mathrm{10}}^{-\mathrm{3}}{f}^{\mathrm{2}}\end{array}$$

where f is in GHz. Here, we test the applicability of this model for L-band and full possible density range in the percolation zone of the GrIS. For f=1.4 GHz (L-band), the values of A₁, A₂, and B₁ are 0.82, 0.96, and 0.24 respectively. Hereafter, we refer to this model as Hallikainen model for simplicity. Ulaby et al. (2014, p. 140–145) used the same formulation (Eqs. 13–14) from Hallikainen et al. (1986), except they scaled the A parameter in Eq. (13a) (i.e., the real part of the ε_eff) with the A₁ factor (<1) from Eq. (14a), as follows. The imaginary part, ${{\mathit{\epsilon }}^{\prime \prime }}_{\mathrm{eff}}$, of Ulaby et al. (2014, p. 140–145) and Hallikainen et al. (1986) remained the same. Here, this model is referred to as the “Ulaby model”.

$$\begin{array}{}\text{(15)}& A=A_{\mathrm{1}}(\mathrm{1}+\mathrm{1.83}{\mathit{\rho }}_{\mathrm{ds}}+\mathrm{0.02}{A}_{\mathrm{1}}{v}_{\mathrm{w}}^{\mathrm{1.015}})+B_{\mathrm{1}}\end{array}$$

Based on observations, Colbeck (1980) revised the PVS mixing theory to derive the dielectric constant of wet snow for three distinct cases. The salient feature of this model is that it permits air, ice, and liquid water to form the continuum environment depending on their volume fraction. When both the density and liquid inclusion are low (ρ_ds<550 kg m⁻³ and v_w<7 %), air is the continuous environment throughout the medium. This regime (Colbeck (1980) case I) called the `pendular regime' where ice grains form clusters and isolated liquid water resides in the fillets and veins of the grain contacts, describes well the liquid water inclusion in the percolation zone of GrIS. The shape of the fillets (thin and longer) and veins (shorter) is represented by their aspect ratio ($n=\frac{c}{a}=\frac{c}{b}$), which can lie between 1 (spherical) and 10 (needle shaped) Colbeck (1980). However, comparing with the measurements, Colbeck (1980) suggested an average value of n=3.5 for this case (v_w<7 %).

As the liquid water inclusions increase beyond 7 % and there is enough pore space (ρ_ds<550 kg m⁻³ and thus porosity, ϕ>0.4), grain clusters break down, and a transition from the pendular regime to another regime called the “funicular regime” occurs when liquid water becomes continuous throughout the pore space containing isolated air bubbles and rounded ice grains. Although other studies (Denoth, 1982, 1989, 1994; Denoth et al., 1984) report that this transition can occur at lower v_w (<7 %), we do not consider this case (Colbeck (1980) case II) as it is not representative of typical melt conditions in the percolation zone. It may represent saturated snow or slush in the GrIS ablation zone.

However, if the density is high (ρ_ds>550 kg m⁻³, high ice fraction and thus low porosity, ϕ<0.4), ice forms the continuum medium and air becomes spherical isolated bubbles, while liquid water still resides in the fillets and veins of the grain clusters for low liquid inclusions (v_w<7 %). This case (Colbeck, 1980; case III) is also relevant to percolation zone firn, especially at depth below the seasonal snow. Therefore, we implement Colbeck (1980) cases I (pendular regime) and III (low porosity) using 3 component PVS mixing theory in the following form.

$$\begin{array}{}\text{(16a)}& {\mathit{\epsilon }}_{\mathrm{eff}}=\begin{array}{ll}\text{Three-phase PVS mixing}& {\mathit{\rho }}_{\mathrm{ds}}\le \mathrm{550}\mathrm{kg}{\mathrm{m}}^{-\mathrm{3}}\\ \text{with air background}& \\ \text{Three-phase PVS mixing}& {\mathit{\rho }}_{\mathrm{ds}}>\mathrm{550}\mathrm{kg}{\mathrm{m}}^{-\mathrm{3}}\\ \text{ with ice background}& \end{array}\end{array}$$

For n=3.5, we used m=0.072 following Picard et al. (2022), where m is the ratio of the depolarization factors.

$$\begin{array}{}\text{(16b)}& m={\displaystyle \frac{N_c}{N_{\mathrm{a}}}}={\displaystyle \frac{N_c}{N_b}}\end{array}$$

Among the power-law based models, we considered following three well known models. The Birchak et al. (1974) model uses an exponent, $\mathit{\beta }=\frac{\mathrm{1}}{\mathrm{2}}$ in the power-law relation of Eq. (8). Sihvola et al. (1985) used a similar model based on Cummings (1952) results and obtained a best fit with β=0.4. The Looyenga model (Looyenga, 1965) also follows the power-law relations of the form in Eq. (8). Specifically, it uses $\mathit{\beta }=\frac{\mathrm{1}}{\mathrm{3}}$. We will refer these models as Birchak, Sihvola, and Looyenga models, respectively. To determine the effective dielectric constant of wet snow, Eq. (8) can be extended for a three-component mixing with respective volume fraction, but we used Mätzler (2006) model for dry snow, then used Eq. (8) for water inclusion in a dry snow environment.

Aside from above models, we considered purely empirical Tiuri et al. (1984) model that used experimental data to measure the complex dielectric constant of wet snow for frequencies 859 MHz–12.6 GHz. Their results suggest that the complex dielectric constant of snow is largely unaffected by the snow structure. In dry snow, the dielectric constant is primarily determined by the density. For wet snow, both the imaginary part (${\mathit{\epsilon }}_{\mathrm{s}}^{\prime \prime }$) and the increase in the real part due to liquid water (${\mathit{\epsilon }}_{\mathrm{s}}^{\prime })$ show a similar dependence on volumetric wetness, which were empirically modelled as,

$$\begin{array}{}\text{(17a)}& {\displaystyle }{\mathit{\epsilon }}_{\mathrm{s}}^{\prime }=\mathrm{1}+\mathrm{1.7}{\mathit{\rho }}_{\mathrm{ds}}+\mathrm{0.7}{\mathit{\rho }}_{\mathrm{ds}}^{\mathrm{2}}+\mathrm{8.7}{v}_{\mathrm{w}}+\mathrm{70}{v}_{\mathrm{w}}^{\mathrm{2}}\text{(17b)}& {\displaystyle }{\mathit{\epsilon }}_{\mathrm{s}}^{\prime \prime }={\displaystyle \frac{f}{{\mathrm{10}}^{\mathrm{9}}}}(\mathrm{0.9}{v}_{\mathrm{w}}+\mathrm{7.5}{v}_{\mathrm{w}}^{\mathrm{2}}),f=\mathrm{500}-\mathrm{1000}\mathrm{MHz}\end{array}$$

A summary of the above-mentioned wet snow dielectric mixing models is given in Table 1.

Table 1 Salient features of the ten selected dielectric mixing models.

2.4 Theoretical Penetration Depth

An important quantity of interest for liquid water quantification is the depth of penetration (also known as e-folding depth, δ_p), a depth at which the signal power drops to $\frac{\mathrm{1}}{e}$ times (∼37 %) of its initial power at a reference location due to absorption and scattering in the snow and firn. The effective depth from which microwave radiometers receive emissions is usually higher depending on the medium properties although the signal strengths progressively diminish (less than 3 % and 5 % of their initial value at depths 5δ_p and 3δ_p, respectively). The actual depth of penetration also depends on the signal to noise ratio (SNR) as well as the precision of the radiometer instruments. To estimate the total liquid water amount, the radiometer should receive emissions from the full wet layer. The signal power of an electromagnetic wave propagating through snow/firn is determined by the extinction coefficient κ_e of the medium, which is the sum of the volume absorption and scattering coefficient respectively, ${\mathit{\kappa }}_{\mathrm{e}}={\mathit{\kappa }}_{\mathrm{a}}+{\mathit{\kappa }}_{\mathrm{s}}$. Since the snow grains are much smaller than the L-band wavelength, κ_s≪κ_a, the absorption coefficient dominates the extinction, κ_a≈κ_e. Therefore, neglecting scattering losses at L-band for a low volume fraction of liquid water (${\mathit{\epsilon }}^{\prime \prime }\ll {\mathit{\epsilon }}^{\prime }$) in snow/firn, the penetration depth can be approximated following Elachi and Zyl (2021) and Ulaby et al. (2014, p. 140–145) as:

$$\begin{array}{}\text{(18)}& {\mathit{\delta }}_p={\displaystyle \frac{\mathrm{1}}{{\mathit{\kappa }}_{\mathrm{a}}}}\end{array}$$

where κ_a is the wet snow power absorption coefficient given by κ_a=2α and α is the attenuation coefficient (Np m⁻¹) defined by,

$$\begin{array}{}\text{(19)}& \mathit{\alpha }=-k_{\mathrm{0}}\cdot \text{Im}\left(\sqrt{{\mathit{\epsilon }}_{\mathrm{eff}}}\right)\end{array}$$

where k₀ is the wave number in vacuum, $k_{\mathrm{0}}=\frac{\mathrm{2}\mathit{\pi }f}{c}$, c is the speed of light, and f is the frequency in Hz. Therefore, for a given frequency, δ is determined by the effective dielectric constant, depending on the average volume fraction of liquid water content and the density of snow/firn. For L-band, the penetration depth in dry snow is significantly higher (>100 m) (Matzler et al., 1984) depending on the density. However, in wet snow, the liquid water inclusion significantly increases absorption (ε^′′) thus decreasing δ.

2.5 Liquid Water Amount

For a volume fraction of liquid water v_w (percent) with a wet layer thickness of t_wet, (m), the LWA is calculated by the product of the two,

$$\begin{array}{}\text{(20)}& \text{LWA}=v_{\mathrm{w}}{t}_{\mathrm{wet}}\mathrm{m}\mathrm{w}.\mathrm{e}\end{array}$$

We chose to express the LWA in [mm], which is equivalent to [kg m⁻²] (because the density of water is 1000 kg m⁻³). The LWA represents the vertically integrated liquid water content within the SMAP grid point at that timestamp, corresponding to the SMAP effective sensing depth.

2.6 Liquid Water Retrieval Algorithm

For the LWA retrieval, we iteratively used an inversion-based framework, first minimizing a cost function between the simulated and mean observed TB measured at vertical (p=V) polarization during the frozen season which we considered to span 1 January–31 March and 1 November–31 December for the pre- and post-summer seasons, respectively, across the GrIS percolation zone. For the observations, we used the SMAP enhanced-resolution data products generated using the radiometer form of the Scatterometer Image Reconstruction (rSIR) algorithm (Long et al., 2019; Brodzik et al., 2021) posted on the EASE-2 3.125 km grid (Brodzik et al., 2012, 2014). The rSIR technique utilizes the measurement response function (MRF) of each sample and combines the overlapping (at close but different acquisition times) MRFs to reconstruct an enhanced-resolution TB image (Early and Long, 2001; Long, 2019; Long et al., 1993; Long and Brodzik, 2016; Long and Daum, 1998). The data product provides the twice daily sampling of GrIS in the form of combined morning and evening passes. The advantage of this rSIR processing is that it improves the overall effective resolution of the measurements by about 30 % as compared to the original data products (Long et al., 2023; Zeiger et al., 2024). The radiometric precision of the SMAP original data is within 0.5 K (Chaubell et al., 2018, 2020; Piepmeier et al., 2017).

Using the average measured density from the top 3 m of snow as recorded from the PROMICE or GC-Net AWS, the algorithm first optimizes the pre-summer baseline emission using the pre-summer mean frozen season (1 January–31 March) TB. If the post-summer mean frozen-season TB (1 November–31 December) is lower than the pre-summer mean frozen-season TB – such as due to crust formation caused by refreezing – the post-summer frozen-season TB is used in the optimization to determine a separate reference for the late melt season. The transition between frozen references is identified by the day on which the maximum TB is observed, as refreezing becomes the dominant process thereafter. However, if the post-summer mean frozen-season TB is higher than the pre-summer mean, indicating a warmer background from remnant melt or latent heat, the pre-melt frozen reference is used throughout the year. With these initial conditions, the melting TB in the summer season is a nonlinear function of the wet layer thickness (t_wet), liquid water content (v_w), and melt-related and other snow firn metamorphisms. Here, we did not account for melt induced snow metamorphism in the forward simulation, except for the adjustments in the reference TBs. We then used the melt season observed TBs to derive an average wet layer thickness (t_wet) and liquid water content (v_w) in a two-step optimization process. To remove any spurious melt during frozen season, we also derived a threshold-based binary melt flag. The threshold was determined by an algorithm following Torinesi et al. (2003):

$$\begin{array}{}\text{(21)}& \text{Th}={\text{TB}}_{\mathrm{ref}}+m\cdot \mathit{\sigma }\end{array}$$

where TB_ref is frozen reference, σ is the standard deviation of the TB during the pre-summer reference period, and m is an empirically derived constant. For SMAP V-pol TBs over the Greenland percolation zone, we found m=10 to be optimal.

Figure 1 Configuration of a simplified three-layer ice sheet model to represent equivalent snow and firn stratigraphy for forward modeling of the brightness temperature.

The near-surface density profile in the percolation zone is highly variable and characterized by multi-scale fluctuations (Rennermalm et al., 2022; Johnson et al., 2014, and the references therein). Seasonal melting and refreezing further complicate this structure, leading to the formation of features such as random ice layers and ice pipes. Accurately modeling these effects across the percolation zone remains a significant challenge due to the lack of detailed ancillary data (e.g., temperature, density) and is an ongoing area of research. These sub-grid-scale structural variabilities contribute to the substantial scattering of L-band TBs, particularly during the frozen season (Hossan et al., 2025). Even with measured ice core profiles, simulating L-band TBs remains difficult due to these complexities.

We implemented a simplified three-layer ice sheet configuration (Fig. 1) to simulate TB using MEMLS V3 (Mätzler and Wiesmann, 2012). In this configuration (Fig. 1), the top layer represents the seasonal snowpack, dry snow/firn during frozen season and wet snow during melt. The bottom layer of the three-layer configuration is treated as semi-infinite ice. To account for the influence of internal firn stratigraphy, we introduced a middle layer modelled as an effective dielectric slab. This layer captures the bulk reflective effects of complex subsurface layering in background emission, following a similar approach used by Mousavi et al. (2022). Rather than representing any specific physical layer, it serves to model the cumulative dielectric contrast beneath the seasonal snowpack. Its properties remain fixed over time at each grid point and have insignificant influence on liquid water retrievals, which are governed by the top layer's dynamic properties.

In MEMLS each layer is defined by its thickness (d), physical temperature (T), density (ρ), volumetric liquid-water content (v_w), exponential correlation length (l_ex), and salinity (S). Since our objective is to evaluate the relative behaviours of dielectric mixing models, we tried to make the model as simple as possible by assuming reasonable ranges and fixed values of possible parameters. The top layer has thickness 0.1–20 m, volumetric liquid-water content of 0 %–6 %, and a fixed temperature of 250 and 273.15 K when m_v=0 (dry) and m_v>0 (melt) respectively.

While, L-band TB is sensitive to dry snow density (Hossan et al., 2025; Houtz et al., 2019), we used the average measured density in the top 3 m to better constrain the retrievals. Considering the insensitivity of L-band measurements to snow microstructure, we set l_ex=0 mm for all three layers (Schwank et al., 2014). Planar interfaces and specular reflection are assumed. For ice sheets, the salinity can also be set to 0 ppt. The range of volume fraction of liquid water was determined based on earlier Experiments (Coléou and Lesaffre, 1998; Colbeck, 1974), which suggest that the irreducible water saturation because of capillary retention ranges between 6.5 %–8.5 % of the pore volume depending on the density. Considering snow/firn density in the percolation zone, we determined the maximum volume fraction is within 6 %. Exceptions to this are the saturated water such as buried and open lakes (Dunmire et al., 2021), firn aquifers (Miller et al., 2022b, 2020b), which are atypical for this area and we did not include these cases into our consideration

The middle layer has a fixed thickness of 5 m. It shares the same physical temperature as the top layer during the frozen season and uses a fixed temperature of 265 K when the top layer contains liquid water (at melting point). The real part of its relative dielectric constant is tuned to match the simulated TB with the observed TB during frozen season while the imaginary part is fixed 0.0002 (same as ice loss factor). The bottom layer is semi-infinite ice with fixed density (917 kg m⁻³) and physical temperature of 255 K regardless of dry and melt season. We also considered the cosmic background radiation (T_c=2.7 K). However, we did not consider any correction for the atmospheric contribution, because it is small compared to the melt signal (approximately at most 2 K at L-band frequencies (Houtz, et al. 2019)). For a fixed location, the same layer characteristics were applied to all the models.

MEMLS utilizes the six-flux theory to model volume scattering and absorption. It also accounts for the effects of radiation trapping due to internal reflection, along with the coherent and incoherent reflections at the layer interfaces (Mätzler and Wiesmann, 2012). The model uses an empirical approach to compute the scattering coefficient, while the absorption coefficient, refraction, and reflection at layer boundaries are derived using physical models. For effective permittivity, by default, the latest version (V3) of MEMLS considers wet snow as a two-phase mixture of prolate ellipsoidal liquid water inclusions in a dry snow background and uses the MG dielectric mixing rule with depolarization factors from Hallikainen et al. (1986) as mentioned in the previous section. We refer to this default configuration of MEMLS V3 as simply Mätzler model. For other models, we used the same setup and input parameters, except we changed the formulations for the complex wet snow dielectric constants (Sect. 2.3).

2.7 LWA Estimates from a Surface Energy and Mass Balance Model

As references, we considered independent LWA estimates from two ice sheet SEMB models, namely the Samimi et al. (2021) SEMB model, which was locally calibrated (Samimi et al., 2020, 2021; Ebrahimi and Marshall, 2016) and the Glacier Energy and Mass Balance (GEMB) model (Gardner et al., 2023), within the NASA Ice-sheet and Sea-Level System Model (ISSM). Both models were forced with the in situ automatic weather stations (AWS) measurements from the Programme for Monitoring of the Greenland Ice Sheet (PROMICE) and Greenland Climate Network (GC-Net) located in the percolation zone of the GrIS. These models used averaged hourly observations of air temperature, air pressure, upwelling and downwelling short and longwave radiation fluxes, snow-surface height, wind speeds (Fausto et al., 2021) along with subsurface profiles of temperature, density, and stratigraphy for initializations (Vandecrux et al., 2024). The SEMB models determine the net energy available for melting if the surface temperature is at the melting point, otherwise for warming or cooling the snow in the upper layer. The subsurface temperature and density then evolve within a one-dimensional model, which is coupled with hydrological processes like meltwater infiltration, refreezing, and retention within the firn. Although the two models under consideration used the same forcing, they use separate parameterizations for these physical processes and a separate model configuration. We refer the reader to Samimi et al. (2021), and Gardner et al. (2023), for specific model details. Despite limitations, these SEMB models are currently the most viable way of validating satellite retrievals. For comparison, we consider these two models individually as well as their ensemble (average).

3 Results

3.1 Effective Complex Dielectric Constant

Liquid water increases both the real and imaginary part of the dielectric constant of snow/firn. Since the dry snow has negligible loss factors at microwave frequencies, almost all the changes in imaginary parts come from liquid water inclusion. Figure 2 shows the change of the complex dielectric constant at a fixed density of 400 kg m⁻³ as function of the volume fraction (v_w) of liquid water inclusions for up to 6 % as appropriate for the percolation zone.

There are large spreads between the models for both the real (Fig. 2a and b) and imaginary parts (Fig. 2c and d); this spread increases as v_w increases. The deviations between the models are higher for the imaginary part than in the real part. For v_w<2 %, the models' agreements for the real part (Fig. 2b) are consistent within two tenths of the relative dielectric constant, except the Hallikainen model (dashed black line), which appeared an outlier. The Ulaby model (solid blue line) uses a scaling factor of A<1, resulting in the lowest estimate among the models up to about 2 % (Fig. 2b); at higher v_w values it provides an intermediate estimate (Fig. 2a). The Ulaby, Tinga, and Debye-like models provide real part of dielectric constant lower than that of even dry snow (the dashed grey line indicates the permittivity of dry snow at −0.5 °C with the Mätzler, 2006 model) for up to v_w 1.2 %, 0.7 %, and 0.4 %, respectively. The Debye-like model and the low frequency approximate of the Hallikainen model closely agree with the empirical Tiuri model, and both lie in the upper end for the v_w>2 % range, while the Mätzler model (solid black line) predicts an intermediate result for the entire v_w range. The exponential models (Looyenga, Sihvola, and Birchak) reasonably agree with the median model (e.g., Mätzler) for low v_w (<2 %) and stay in the lower end of the curves for the higher v_w. Among the structure dependent models, the Colbeck and Tinga models provide relatively lower estimates of the real part of the dielectric constant and agree with the exponential models, especially for v_w>2 %.

Figure 2 Change of real (a, b) and imaginary (c, d) parts of complex dielectric constant of snow/firn of a fixed density of 400 kg m⁻³ as function of volume fraction (v_w) of liquid water content. The right panel shows a zoomed version of the left panel for v_w range 0 %–2 %. As reference, real (a, b) and imaginary (c, d) parts of complex dielectric constant of dry snow at the same density and snow temperature of −0.5 °C are shown by horizontal grey dashed lines.

For the imaginary part of the snow/firn dielectric constant, the Debye-like, Hallikainen, and Ulaby models, which are the same group of models with modifications for frequency dependencies, generally follow the empirically derived Tiuri model almost for the entire range of v_w under consideration. However, it is worth noting that small differences in the loss factor, especially in the lower end, can result in significant differences in terms of TB and depth of penetration. The Hallikainen and Ulaby models are the same for the imaginary part, overlapping with each other. The Tinga model provides the highest estimate of the loss factor for v_w up to about 2 %, then it falls exponentially for the higher end. Mätzler model again provides an intermediate estimate of the loss factor for the entire v_w range under consideration. The exponential models stay in the lower end, as in the case of the real part, with the lower value of β, giving the lower estimate. The Colbeck model results consistently in the lowest value of the imaginary part of the dielectric constant for the entire range of v_w. Since the density is less than 550 kg m⁻³ for these curves, it included only Case I (pendular regime) of Colbeck (1980).

Figure 3 Penetration depth at L-band (1.41 GHz) for a snow/firn density of 400 kg m⁻³ as function of volume fraction (v_w) of liquid water content (a). (b) a zoomed version of (a) for v_w range 3 %–6 %.

3.2 Penetration Depth

The differences in the imaginary part of the dielectric constants are manifested in the penetration depth, an important variable for liquid water quantification. Figure 3 illustrates the penetration depth of L-band (1.41 GHz) signals estimated by the models for a snow/firn density of 400 kg m⁻³ as a function of v_w in the 0 %–6 % range (Fig. 3a) and the 3 %–6 % range (Fig. 3b). All models predict an exponential decay of penetration, but they exhibit substantial differences with respect to one another, though the range generally reduces with increasing v_w. For v_w of 1 %, 3 %, and 5 %, the model estimates of penetration depth range between 2.8–12.8, 1–4, 0.5–2.3 m respectively. The Tinga and Ulaby models provide the lowest estimate of penetration depth for v_w <2 %, and v_w>2 % ranges respectively, while the Colbeck model gives the highest estimate for the entire v_w range as it estimates the lowest loss factor among all the models. The Debye-like, Hallikainen, and Ulaby models closely follow each other. Mätzler model provides an intermediate estimate of penetration depth. The empirical Tiuri model aligns with Mätzler model for low v_w (<1 %); however, it matches better with Debye-like, Hallikainen, and Ulaby models for higher v_w (>2 %). The exponential models, consistent with their complex dielectric constant, lie in between.

3.3 Simulated Brightness Temperature

TBs simulated with these models for a typical representative snowpack in the percolation zone are illustrated as function of v_w (bottom x axis) and LWA (top x axis) in Fig. 4 for wet layer thicknesses of 1 m (Fig. 4a, d), 2 m (Fig. 4b, e), and 3 m (Fig. 4c, f), respectively. For the V-pol, the Tinga model appears to be the most sensitive for low v_w and LWA, then it gradually slows down as the LWA increases, when the Ulaby model provides the highest TB. The Debye-like model closely follows the Ulaby model; the Hallikainen model, which uses the same imaginary part as the Ulaby model, but a higher real part of the dielectric constant, shows lower TB estimates. The difference that is also function of LWA, is more obvious in the H-pol results (Fig. 4d–f). The Tiuri and Mätzler models produce higher TB projections than the Ulaby and Debye-like models for the lower range of LWA, but it flipped for higher range of LWA, and the transition depends on the thickness of the wet layer. In line with the complex dielectric constant, the Colbeck model provides the lowest estimates for almost the entire LWA range under consideration (except for LWA >150 mm, H-pol), and the Birchak, Sihvola, and Looyenga models offer moderate values for all cases. Although, the changes are more pronounced in the H-pol TB, the trends with the LWA are similar except that the saturation in TB occurs for relatively lower LWA compared to the V-pol signals, especially for thicker wet layers (Fig. 4e–f).

The results depend on the density of the dry snow (porosity), which are shown for three different densities (200 kg m⁻³ (Fig. 5a, d), 400 kg m⁻³ (Fig. 5b, e), and 600 kg m⁻³ (Fig. 5c, f)) for a fixed thickness of wet layer (2 m) in Fig. 5. The Debye-based models (Debye-like, Hallikainen, and Ulaby) along with the Tiuri model show significantly higher sensitivity with v_w and thus provide lower estimates of LWA than the Mätzler and Tinga models for the low snow density (200 kg m⁻³) at both V- and H-pol results. However, for high snow density (600 kg m⁻³), this is reversed, Tinga and Mätzler models exhibit higher sensitivity and provide lower estimates of LWA than the Debye-based and Tiuri models, while at intermediate density (400 kg m⁻³), they agree closer for both V- and H-pol TB. Although the sensitivity of the rest of the models varies with the dry snow density, they consistently demonstrate lower sensitivity and provide higher estimates of LWA than the above-mentioned models across the complete density range.

Figure 6 depicts simulated TB, like Fig. 4, but as a function of wet layer thickness in (bottom x axis) and LWA (top x axis) at a fixed snow/firn density of 400 kg m⁻³ for three cases of fixed v_w of 1 % (Fig. 6a and d), 2 % (Fig. 6b and e), and 3 % (Fig. 6c and f). In a broader perspective, the trends of TB with the thickness of the wet layer at a fixed v_w are similar to the TB trends with v_w at the fixed thickness of the wet layer, as presented in Fig. 4. TB grows exponentially with both v_w and t_wet, where each model has a different growth factor, which also depends on v_w and t_wet themselves along with dry snow density and other background conditions.

Figure 4 Vertically (a–c) and horizontally (d–f) polarized brightness temperature at L-band for a snow/firn density of 400 kg m⁻³ as function of volume fraction (v_w) of liquid water content in percent (bottom x axis) and total liquid water amount in mm (top x axis) simulated with MEMLS using different wet snow mixing models for three wet layer thickness: 1 m (a, d), 2 m (b, e), and 3 m (c, f).

Figure 5 Vertically (a–c) and horizontally (d–f) polarized brightness temperature at L-band for a wet layer thickness of 2 m as function of volume fraction (v_w) of liquid water content in percent (bottom x axis) and total liquid water amount in mm (top x axis) simulated with MEMLS using different wet snow mixing models for three snow/firn densities: 200 kg m⁻³ (a, d), 400 kg m⁻³ (b, e), and 600 kg m⁻³ (c, f).

Figure 6 Vertically (a–c) and horizontally (d–f) polarized brightness temperature at L-band for a snow/firn density of 400 kg m⁻³ as function of wet layer thickness (bottom x axis) and total liquid water amount in mm (top x axis) simulated with MEMLS using different wet snow mixing models for three fixed volume fraction (v_w) of liquid water content: 1 % (a, d), 2 % (b, e), and 3 % (c, f).

3.4 Brightness Temperature Sensitivity to Liquid Water Change

The sensitivity of TB to LWA change decreases with increasing LWA. We compute the change of TB for every 1 mm change in LWA, which is shown for V- and H-pol TB in Fig. 7. Here we considered a v_w of 3 %, and increased the wet layer thickness from 0.1 to 5 m. The sensitivities at H-pol are higher than at V-pol for all models. The sensitivity of TB to the change of LWA decays exponentially across all models, falling below 1 K mm⁻¹ at <50 mm of LWA for V- and H-pol. For models that demonstrate higher sensitivity for lower LWA, the sensitivity declines more sharply and saturates at relatively lower LWA (after which they show negative sensitivity, i.e., TB decreases as LWA increases; however, we did not consider negative sensitivity regime in this manuscript as this happens at oversaturated LWA not typical for the percolation zone of the Greenland Ice Sheet, see Sect. 1).

The Ulaby and Tiuri models show the highest sensitivity at lower LWA, closely followed by the Hallikainen, Tinga, and Mätzler models. However, the sensitivities of these models fall below tenths of K mm⁻¹ for LWA >70 mm for both V- and H-pol. The Birchak, Sihvola, and Looyenga models demonstrate moderate sensitivities, while the Colbeck model presents the lowest sensitivity among the models for lower LWA; however, their sensitivities also decrease slowly with LWA, and they remain reasonably sensitive for the higher end of the LWA. Although, the first group of models shows almost negligible sensitivities close to or beyond 100 mm of LWA, no models showed perfect 0 or negative sensitivities within 150 mm of LWA. However, it is obvious that for a majority of the models the uncertainty of the retrievals at LWA >60–70 mm will be significantly higher.

Figure 7 L-band vertically (a) and horizontally (b) polarized brightness temperature sensitivity (change of TB in K per mm change in liquid water amount) at a snow/firn density of 400 kg m⁻³ as of total liquid water amount in mm simulated with MEMLS using different wet snow mixing models for a fixed 3 % volume fraction (v_w) of liquid water content with varying wet layer thickness. Top panel shows the results for a LWA range 0–50 mm, and the bottom panel shows the results for an extended range (0–100 mm).

3.5 LWA Retrievals

Figure 8 compares the simulated V-pol TBs with the observed V-pol SMAP TB time series (blue dash-dotted line) along with the frozen season TB references (dashed grey lines), during the 2023 melt season at six PROMICE and GC-Net AWS. The AWS were selected based on their varying geographic locations in the percolation zone, climatic record of melt, and in situ data availability for the validation. Table 2 shows the geographic locations, elevations along with the mean 2 m air temperature, the mean dry snow density of the upper 3 m, and the mean pre-summer frozen season V-pol L-band TB at these sites.

Figure 8 L-band vertically polarized observed and simulated brightness temperature time series during 2023 melt season at six selected PROMICE and GC-Net within the GrIS automatic weather stations (AWS) locations within the percolation zone. Frozen season TB references are shown as dashed grey lines. In panel (f), the large downward spikes are processing artifacts.

The 2023 melt year experienced above average melting (Poinar et al., 2023), though the annual average temperatures reflect both a cooler accumulation period and warmer melt season (Poinar et al., 2023). All sites have very stable but different frozen season TB, representative of their different subsurface backgrounds. Some occasionally decreasing spikes before and after the melt season at the SDM site are artifacts of the rSIR processing (Long and Brodzik, 2024, confirmed through personal communications). Since the simulated TBs were optimized with the mean observed TBs during frozen season, they are less affected by these artifacts during frozen season. Throughout the summer season at each of the AWS location, the simulations with different models closely align with the observations with no significant bias. The agreements are so close that the simulated TBs are almost overlaid by the observed TBs (Fig. 8). But this was obtained with different combinations of t_wet and v_w, resulting different LWA.

CP1 and DY2 are both perfect representatives of the Greenland percolation zone, with moderate upper layer density and numerous ice layers and pipes due to annual refreezing of seasonal melt (Jezek et al., 2018; Vandecrux et al., 2024). These ice layers significantly attenuate microwave emissions from deeper layers giving very low frozen season TBs (∼148 and 144 K, respectively). At the same time, these sites also provide high TB sensitivity to liquid water during peak melt season as the effective emissions approaches close to unity (>0.95, see Fig. 8a–b). At both AWS locations, the TBs remained elevated through October compared to their pre-summer frozen references.

Table 2 Summary of the six selected PROMICE and GC-Net automatic weather station (AWS) locations. The location information was adopted from PROMICE and GC-Net Automatic Weather Station metadata, available at https://promice.dk (last access: 27 November 2024). Mean annual, summer (June–August) air temperature, mean dry snow density, and mean frozen season brightness temperature were calculated using 2023 data.

Table 3 Average thickness of wet layer (in cm) during 2023 melt season (May–September) retrieved by ten dielectric mixing models with SMAP observations and two surface energy and mass balance models forced by in situ observations at the six selected PROMICE and GC-Net automatic weather station (AWS) locations.

Table 4 Maximum summer melt (in mm) during 2023 melt season (May–September) estimated by ten dielectric mixing models with SMAP observations and two surface energy and mass balance models forced by in situ observations at the six selected PROMICE and GC-Net automatic weather station (AWS) locations.

KAN_U is located at the lowest elevation (1848 m) of all the sites close to the equilibrium line in the southwestern Greenland. Air temperatures are often above freezing during the melt season (Table 2), and the region experiences substantial surface melting. As a result, ice layers are thicker, and the near surface densities are high with low variability. Frozen season TB is the highest of all sites examined, and during 2023 melt season (Table 2), TB is possibly saturated due to extensive persistent melt that keeps TB elevated beyond the end of October (Fig. 8c).

NSE, SDL, and SDM are located at high elevation in southeast Greenland. These locations generally receive more accumulation and less melt than the other stations examined here (Fausto et al., 2021). Upper layer densities are low to moderate with lesser number of ice layers. This is revealed by their moderate means frozen season TBs (Table 2). Contemporaneous large summer peaks of TB at these three sites in 2023 are indicative of melt events (Hossan et al., 2025). However, the duration of the melt events at these AWS locations is substantially shorter than at the previous three AWS locations. In addition, the post-summer mean frozen-season TBs drop below their pre-summer frozen references (note the changes in the reference TBs, Fig. 8d–f), unlike at the previous three sites. This is because the post-melt temperature at these locations drops sharply, possibly due to crust formation that enhances internal reflections.

Table 5 Pearson linear correlation coefficient between LWA estimate by each of the dielectric mixing model and their ensemble with SMAP observations during 2023 melt season (May–September) and corresponding LWA estimate obtained by averaging Samimi and GEMB surface energy and mass balance models forced by in situ observations at the six selected PROMICE and GC-Net automatic weather station (AWS) locations.

LWA retrieved from L-band TB using different dielectric mixing models and two SEMB models for 2023 melt season at the selected AWS locations are presented in Fig. 9. For all the mixing models, the retrieved average thicknesses of wet layer used in the LWA retrieval are given in Table 3.

Regarding the onset of melt season at CP1, the satellite retrievals with different mixing models are the same (24 June) and reasonably contemporary with the SEMB models considering noise levels of in situ instruments. We used the SMAP melt flags (see Sect. 2) to remove spurious melts in winter. Similar flagging for SEMB models is difficult to find. Some studies have used thresholds in LWA (e.g., >2 mm in Hossan et al., 2025; Zhang et al., 2023; Leduc-Leballeur et al., 2020; Van Den Broeke et al., 2010). However, we did not use such threshold because dielectric mixing models that show higher sensitivity for light melt events can result LWA within 2 mm. Hence, here we rather focus on relative intensities and durations of significant and consistent melt events qualitatively, ignoring suspicious melt events, like one in early June at CP1 where GEMB estimated a short-lived melt (∼2.5 mm in amount) but the SAMIMI model like all the satellite retrievals indicates no melt even when both models used the same in situ measurements.

Table 6 Mean root mean squared differences (RMSD in mm) between LWA estimate by each of the dielectric mixing model and their ensemble with SMAP observations during 2023 melt season (May–September) and corresponding LWA estimate obtained by averaging Samimi and GEMB surface energy and mass balance models forced by in situ observations at the six selected PROMICE and GC-Net automatic weather station (AWS) locations.

All three methods (SAMIMI, GEMB, and satellite retrievals) agreed qualitatively on three main persistent melt events in terms of relative intensities and duration at CP1 during the melt season: a small one in late June, followed by the major melt event that sustained whole July through late August, and a moderate one that begun in late August (on top of the sub-surface remnant melt), lasting through mid-September based on SMAP or early October based on the SEMB models (Fig. 9a). For the early and late season melt events (light to moderate amount), the SEMB models estimated more LWA than any of the satellite retrievals with GEMB surpassing SAMIMI. However, during the peak melt event, we observed mixed results with the dielectric mixing models when comparing them to the SEMB models: the Colbeck and Looyenga models estimate higher LWA (Table 4) than the SEMB models, while the Tiuri, Ulaby, Debye-like, Hallikainen, Mätzler, and Birchak models show lower estimates of LWA compared to the SEMB models (Table 4). These retrievals are lower than the ones Hossan et al. (2025) presented using the Ulaby model at this site during the same melt season. Different constrains in frozen and melt season parameters, mainly density and background temperatures, explain some of these differences.

The Sihvola model was found to be in closest agreement with the SEMB models at CP1 until the peak of the melt season. Afterward, when active melting at the surface stops (as evidenced by gradual loss of LWA) and meltwater percolates deeper and refreezes, the satellite retrievals and the SEMB models exhibit more substantial differences (Table 6), which impacted the overall correlation and RMSD ($\mathrm{0.74}\le r\le \mathrm{0.84}$ and 18 mm ≤ RMSD ≤ 28 mm; see Tables 5–6). All the satellite retrievals consistently indicate a faster refreezing rate of the subsurface liquid water than in case of both the SEMB models, where the refreezing is determined by the evolution of the density profile and thermal conduction.

Figure 9 Comparison of the total daily liquid water amount (LWA) estimated from SMAP L-band TB observations using ten dielectric mixing models, with estimates from the SAMIMI EBM (orange dash-dotted line) and GEMB (pink dotted line), both forced with PROMICE and GC-Net AWS in situ measurements for 1 June–31 October 2023.

The melt trends at DY2 (Fig. 9b) are similar to CP1 with some differences (Tables 3–6). The early season melt in late June is minor, while the late season peak is relatively higher with gradual loss of persistent subsurface melt that extended even beyond the end of September in consensus with the mixing models and with the SEMB models. Under such persistent liquid water and warmer subsurface background, while the surface recommences melting, liquid water is expected to infiltrate deeper and form a thicker wet layer. Satellite retrievals of average wet layer thickness by majority of the mixing models support this (Table 3). The overall agreement between the SEMB models and satellite retrieval with different mixing models is better at this site ($\mathrm{0.95}\le r\le \mathrm{0.98}$ and 5 mm ≤ RMSD ≤ 17 mm; see Tables 5–6).

KAN_U is known to undergo extensive LWA throughout the summer (Hossan et al., 2025). However, all the mixing models, including SEMB models, report relatively lower LWA but prolonged melting conditions (Fig. 9c). All the mixing models, except the Tinga model estimated higher average thickness of the wet layer than the SEMB models (Table 3). The snowfall at this site is climatologically lower (MacFerrin et al., 2019; Machguth et al., 2016) and thicker ice layers underneath ought to prevent deeper infiltration. Regarding LWA, the SEMB models showed better alignment with the Colbeck, Looyenga, and Sihvola model-based retrievals at the beginning of the peak melt (early July). However, as the melt season progresses and the snow/firn profile evolves, the differences intensify with the previous models, rather they better match with the rest of the models which claim lower estimates of LWA throughout. Firn models push the liquid water out of the system (called runoff) if the water balance exceeds certain limit (irreducible water saturation) determined by the available pore space. In reality, this excluded liquid water must still exist somewhere, which may explain some of these misalignments (Tables 4–6). Nevertheless, the spread of maximum summer melts between the mixing models (satellite retrievals) is large (min 21 mm by Tinga and Mätzler models to max 60 mm by Colbeck model with SD = 15 mm; see Table 4).

Since NSE is located at higher elevation (Table 1), historically it receives less frequent and less intense melt. In the 2023 melt season, we, however, observed similar melt trends with shorter duration compared to the previous three sites (Fig. 9d). Only the GEMB model detected the early season melt in late June. Both GEMB and SAMIMI models estimate the presence of liquid water (max. 18 and 7 mm respectively) in the late melt season, which were also detected by all the satellite retrievals (mixing models) that better align with the estimate of SAMIMI model. Compared to the pre-summer mean frozen TB, the post-melt mean frozen TB dropped around 5 K at this site possibly from ice layers formed by refreezing of summer melt, common in this higher elevation areas (Hossan et al., 2025). Without dynamic thresholding, these late-season, less intense melt events would have been missed as false negatives. Other than this, although lesser in intensities, the order of magnitude of the satellite retrievals remains the same. However, the results of the SAMIMI model better match with the group of mixing models that provide a lower estimate of LWA while the GEMB model provides the upper limit and better aligns with the models that indicate higher LWA. The GEMB model also refreezes noticeably more slowly than the SAMIMI model, refreezing faster and aligning better with the satellite retrievals.

Though all the satellite retrievals and SEMB models estimated an overall slightly higher LWA at SDL site in 2023, the trends closely replicated that of NSE (Fig. 9e). The late season melt is however now stronger at this site, and the threshold algorithm detected melt and LWA was quantified by all the mixing models with the similar sequence – the Tinga model giving the lowest and the Colbeck model giving the highest maximum summer LWA (Table 4). Compared to the maximum summer LWA, the average thickness of the wet layer (Table 3) was found to be overall higher at this site, similar to NSE, which is anticipated since there is enough pore space but a lower number of ice layers allowing liquid water to percolate deeper. The SEMB models, despite differences between them, better align with the Sihvola, Looyenga, and Colbeck models than the rest of the others (Table 6), but both retain subsurface liquid water for elongated periods in the fall but indicate a thinner thickness of wet layer compared to the satellite retrieval with the majority of the mixing models (Table 3).

At the SDM site, no significant melt was determined until the beginning of July, neither by the satellite retrievals nor by the SEMB models. The late season melt was also insignificant. The two SEMB models closely resemble at this site in phase and magnitude, but again they disagree with all mixing models in the satellite retrievals in the rate of subsurface refreezing in the late melt season, which impacted the comparison metrics in Tables 5 and 6.

4 Discussion

The strong dielectric contrast between dry snow background and liquid water inclusion causes the complex dielectric constant of wet snow to significantly change with v_w, which can vary over a wide range depending on the density (porosity) of the dry snow. Here, we focused on the GrIS percolation zone where typical v_w is known to be approximately 0 %–6 %. Selected dielectric mixing models were found to vary widely over this narrow range, giving large uncertainties in modelling the effective depth of penetrations, TB, and consequently, in quantifying LWA based on the dielectric constant retrieved from satellite measurements. Differences of depolarization factors that describe the shape and orientation of the liquid water inclusion with respect to the emitting EM field mainly contribute to these differences for the structure dependent models. For the power law-based models, their degree controls the higher order local interactions in lieu of depolarization factor. There are significant uncertainties in the penetration depth, or emission contribution depths, between the models, especially when the v_w (and hence the absorption) is low (pendular regime).

MWR TB is on the other hand a non-linear function of multiple parameters that gradually or abruptly vary with depth. Therefore, with limited knowledge of detailed snowpack properties and their evolution, modelling and interpretation of snow/firn microwave radiation, especially at L-band, which is sensitive from the surface to deeper layers, is difficult. Here, we used a simplified profile of temperature, density, and stratigraphy to simulate frozen and melt season TB at L-band. Nevertheless, LWA estimation with a single frequency is an underdetermined problem – even with constrained frozen background parameters, numerous combinations of volume fraction of melt and wet layer thickness can produce the same or close TB, but with different LWA. By using average measured dry snow density and parameterizing the wet layer thickness with an average retrieved thickness over melting days, we attempted to minimize some of the uncertainties.

The Debye-like, Hallikainen, and Ulaby models show higher TB sensitivity to lower v_w at low density snow/firn background and provide lower estimates of overall LWA. The opposite results were observed when the background density is higher (see Fig. 4). This is counterintuitive because when the background density is low, there should be enough pore space in the snow to either hold more v_w or to support deeper percolation or vice versa. The results of these models, however, align most closely with the empirical Tiuri model. This is encouraging in the sense that the Hallikainen and Ulaby models were originally derived and validated for measurements in the 3–37 GHz range and for a limited density range; these agreements support their applicability to L-band applications. However, the Tiuri model assumes that the effective dielectric constant of wet snow is independent of the snow microstructure, whereas these models depend on specific microstructural characteristics. In reality, snow microstructures are often irregular, aggregated, and dynamically evolving.

The Tinga model, on the other hand, shows more consistency with the change of dry snow density, yet it is highly nonlinear with v_w (Figs. 4 and 5). The Mätzler model provided an overall intermediate result in terms of effective dielectric constants (these results are in line with Picard et al., 2022, results), depth of penetration, TB, and LWA, and showed a reasonable fidelity over a wide range of density (Fig. 4). But compared to the SEMB models, which were forced by AWS measurements, Mätzler model, along with the above-mentioned models, provided the lower LWA in all six AWS (Fig. 9 and Table 4).

So, if these models represent TB realistically, saturation at relatively low LWA would limit the liquid water estimation at L-band within a certain limit (approximately no more than 60–70 mm). The Colbeck model has a convincing theoretical and experimental basis, nevertheless this model consistently stayed apart from other models and provided the low end of the effective dielectric constant (and TB), and the high end of the LWA and penetration depth. The power-law dependent models, Birchak, Sihvola, and Looyenga, provided consistent estimates of the LWA and penetration depth throughout (dielectric constant and TB w.r.t. v_w) in the order of lower to higher (higher to lower), respectively. Sihvola model was found to be the best match with SEMB models for the AWS and melt season considered (the overall RMSD at six AWS was ∼11 mm; see Table 6). Although these empirical models are heuristic in nature and lack rigorous physical foundations, they offer the advantage of being easily configurable through the degree of the model (a single parameter to fit – β) and can easily be fitted to the available ground truth (or SEMB estimates).

However, we refrain from recommending any particular model in this article, except exploring and demonstrating their individual and comparative characteristics under different v_w, density, and other firn conditions, because caution should be taken when considering SEMB models as the reference for validating LWA estimates since they have their own limitations. Difference between the SEMB models when forced with the same inputs partly explains this. Consistent delayed refreezing in both the models (SAMIMI and GEMB) when melting at the upper surface ceases is an apparent indication of inaccurate thermal conduction, affecting the overall LWA. Thermal conductivity of wet snow and wet-dry interface is probably lower. Future work should better parametrize these processes to refine the models.

Satellite derived LWA can be attributed to infiltrating water, occurring through two distinct modes of unsaturated flow: the downward propagation of a wetting front and the movement of water through preferential flow paths (also called pipes or flow fingers) (Marsh and Woo, 1984). Under sustained melting, wetting fronts typically form and propagate downward from the surface into the underlying cold firn (Colbeck, 1975), although their advancement can be hindered by nightly refreezing or by snowfall events. Additionally, structural features in the firn, such as ice layers or microstructure contrasts, can trap water. However, as water accumulates, it occasionally breaks through in highly heterogenous locations. The resulting preferential flow paths allow large volumes of water to infiltrate deep below the wetting front, bypassing the cold firn layers (Marsh and Woo, 1984; Pfeffer and Humphrey, 1996). The partitioning of meltwater between infiltration via wetting fronts versus preferential flow paths is highly variable and inherently difficult to predict due to its sensitivity to subtle firn structural and thermal conditions.

Several factors suggest that twice-daily LWA retrievals are more likely to reflect the water associated with the surface wetting front than liquid water contained within deeper preferential flow paths. First, the signal from water in subsurface pipes must propagate through the overlying wet layer, which has a stronger and more coherent L-band response due to its proximity to the surface and higher spatial continuity. Second, flow through preferential pathways is typically event-driven, with the breakthrough of accumulated water quickly penetrating deep into colder firn, where it often refreezes within hours rather than persisting for days or weeks (Humphrey et al., 2012). Finally, the L-band signal inherently averages over broad spatial footprints on the order of kilometres, favouring detection of the spatially extensive and homogeneous surface wet layer over the centimetre- to metre-scale, highly heterogeneous pipe structures. Therefore, future work should better understand and parametrize these processes to refine the algorithm.

Parametrizations in the retrieval framework – such as the assumption of simplistic stratigraphy and liquid water distributions – may affect the absolute LWA estimates. However, these factors are likely to impact all models in a similar proportion. Therefore, the relative differences between the model estimates are more likely attributable to the specific formulations and assumptions of each model. Nevertheless, future work should aim to incorporate more advanced algorithms capable of resolving vertical profiles. Additionally, a spatially and temporally dependent threshold should also be considered in future work to account for not only the ice layers due to refreezing but also the seasonal evolution of the snowpack that obviously contributed to the uncertainty in the results. The challenge would be the sensitivity and saturation of TB with increasing LWA. As shown, the signal power/field intensity decreases exponentially with depth and v_w at a rate determined by the absorption and scattering coefficients (Figs. 3 and 7). Over absorption dominated regions, LWA estimates beyond ∼60–70 mm would be highly uncertain. Future retrievals should consider a wider range of v_w to incorporate the negative sensitivity (scattering dominated regions). To handle the inherent nonlinearities and dimensions of the problem, advanced techniques (such as deep learning) may be beneficial. Upcoming lower frequency missions (e.g., CryoRad: Macelloni et al., 2018) would also offer new opportunities to sense deeper and enhance the capabilities.

5 Conclusion

We investigated performances of ten dielectric mixing models for modelling wet snow TB at L-band to estimate the LWA in snow/firn column in the percolation zone of Greenland ice sheet. Six of the models (Tinga, Debye-like, Hallikainen, Ulaby, Mätzler, and Colbeck) are derivatives of fundamental mixing models (either MG or more generally, PVS) with empirically derived depolarization factors that account for the shape and orientation of the liquid water inclusions in dry snow background with respect to the emitting electromagnetic field. Except for the Colbeck model in this group, all models show overall relatively higher sensitivity of the effective dielectric constant to LWA, and thus TB, and generally produce lower estimates of LWA compared to the SEMB models. Colbeck model displays the lowest sensitivity of LWA to the effective dielectric constant and TB and, hence, yields the highest LWA of all models.

The differences among these models mostly originate from their depolarization factors that depend on multiple factors including density and LWA; they are deemed to be very difficult to quantify. Another group of models that follows power law relationships (Birchak, Sihvola, and Looyenga), not explicitly considering the depolarization factor, exhibit intermediate sensitivity of LWA to the effective dielectric constant and TB, and offer higher LWA than the former group of models (Tinga, Debye-like, Hallikainen, Ulaby, Mätzler, and Colbeck). A lower exponent (β) results in a lower sensitivity and a higher LWA, since a lower exponent allows the background to dominate in the mixing model. The results of the Tiuri model, which is fully based on empirical fitting to the field measurements at around 1 GHz, generally lie with that of the former group of models that explicitly consider the higher order interactions between the liquid water inclusions through depolarization factors. While the Hallikainen and Ulaby models were originally derived and verified for measurements made at 3–37 GHz, these agreements with Tiuri model support their applicability to L-band applications.

Compared to the SEMB models (SAMIMI and GEMB) driven by in situ observations, the first group of mixing models (Tinga, Debye-like, Hallikainen, Ulaby, Mätzler, and Tiuri) estimated consistently lower LWA in five of the six PROMICE GC-Net sites (except KAN_U) which are more typical of percolation zone snow/firn physical conditions. Colbeck and Looyenga models measured consistently higher LWA than the SAMIMI and GEMB models in all six sites. In general, the Sihvola model aligned best with the SAMIMI and GEMB models for 2023 melt season. However, the best match does not imply correctness in the absence of an actual truth estimate. The SAMIMI and GEMB models disagree widely in certain cases; in general, the SEMB models have been found to produce diverging LWA estimates with the same in situ meteorological measurements (e.g., Vandecrux et al., 2020; Hossan et al., 2025; Moon et al., 2025).

Despite the satellite retrievals using the different mixing models showing a wide variance in the total and maximum summer LWA, no significant discrepancies were observed in the timing of the onset and refreezing melt, which is based on the observed TB change. However, although satellite retrieval agrees well with the onset of melt with SEMB models, significant disagreements were found in timing of complete refreezing of sub-surface liquid water in snow/firn. While all L-band retrievals indicate a sharper refreezing in all sites except KAN_U, the SAMIMI and GEMB models seemed to refreeze slowly and retain sub-surface liquid water for an elongated period in post melt season. This was attributed to low thermal conductivity and slow heat transmission in the firn models. The differences between SAMIMI and GEMB models, even when they were run by the same set of in situ observations, are also indicative of the differences in their process representations.

This study sheds light on the behaviour of wet snow dielectric mixing models and consequent TB in presence of low liquid water (v_w: 0 %–6 %, or LWA: 0–120 mm). The sensitivity and saturation behaviour of the models were broadly explored that gives an idea about the uncertainty associated with translating the L-band retrieved effective dielectric constant to LWA. Further work is required to better understand the melt water process in the snow and firn and their interactions with the microwave emissions. More sites specific in situ measurements of firn profiles under various conditions will be the next step to calibrate and validate these models to make better recommendation about using a group of specific mixing models.