Simulation of the microwave emission of multi-layered snowpacks using the Dense Media Radiative transfer theory: the DMRT-ML model

. DMRT-ML is a physically based numerical model designed to compute the thermal microwave emission of a given snowpack. Its main application is the simulation of brightness temperatures at frequencies in the range 1–200 GHz similar to those acquired routinely by space-based microwave radiometers. The model is based on the Dense Media Radiative Transfer (DMRT) theory for the computation of the snow scattering and extinction coefﬁcients and on the Discrete Ordinate Method (DISORT) to nu-merically solve the radiative transfer equation. The snowpack is modeled as a stack of multiple horizontal snow layers and an optional underlying interface representing the soil or the bottom ice. The model handles both dry and wet snow conditions. Such a general design allows the model to account for a wide range of snow conditions. Hitherto, the model has been to simulate the thermal emission of the deep ﬁrn on ice sheets, shallow overlying in Arctic and Alpine regions, and overlying ice on the large ice-sheet margins and glaciers. DMRT-ML thus been vali-dated in three very different conditions: Antarctica, Barnes Ice Cap (Canada) and Canadian tundra. It has been recently used in conjunction with inverse methods to retrieve snow grain size from remote sensing data. The model is written in Fortran90 and available to the snow remote sensing community as an open-source software. A convenient user interface is provided in Python.


Introduction
Passive microwave radiometers on-board satellites acquire useful observations for the characterization of snow-covered areas. These observations are available in these areas several times a day, are relatively insensitive to the atmosphere in many frequency bands, and are independent of the solar illumination. They are sensitive to several properties relevant for monitoring the snow cover and have been exploited in numerous algorithms to retrieve continental snow cover extent (Grody and Basist, 1996), snow depth and snow water equivalent on both land (Josberger and Mognard, 2002;Kelly and Chang, 2003;Derksen et al., 2003) and sea ice (Cavalieri et al., 2012;Brucker and Markus, 2013), snow accumulation on ice sheets (Abdalati and Steffen, 1998;Vaughan et al., 1999;Winebrenner et al., 2001;Arthern et al., 2006), melt events (Abdalati and Steffen, 1997;Picard and Fily, 2006), snow temperature (Shuman et al., 1995;Schneider, 2002;Schneider et al., 2004), and snow grain size (Brucker et al., 2010;Picard et al., 2012). Some of these studies are based on empirical relationships supported by physical interpretations (Koenig et al., 2007) and others directly use physical models and data assimilation techniques (Durand and Margulis, 2007;Picard et al., 2009;Takala et al., 2011;Toure et al., 2011;Huang et al., 2012). In both cases, understanding G. Picard et al.: Simulation of the microwave emission of multi-layered snowpacks and modeling the physical processes of the microwave emission by snow and the underlying surface are crucial.
Modeling the snow microwave emission from snow physical properties takes in general two successive steps. The first step is the computation of the electromagnetic properties (e.g., effective dielectric constant, scattering and absorption coefficients) that characterize the propagation and the single events of interaction between wave and matter. These properties are calculated from the geometric micro-structural properties of the snow, assumed to be homogeneous within a given layer. The second step is the computation of the emission and the propagation through the snowpack by accounting for the multiple interactions of microwaves within the snow as well as by refraction, reflection and transmission at the interfaces. It is often treated with the radiative transfer equation on a plane-parallel medium for which generic solutions exist (Tsang et al., 1985;Fung, 1994). In contrast, the first step is specific to the medium. In the case of snow, the main challenge for the electromagnetic calculation is the high density of scatterers. The volume occupied by the scatterers over the total volume (referred to as the fractional volume f ) is significant -typically larger than 20 % -, which implies that the scatterers strongly interact with each other and the independent scattering theory used for vegetation or clouds (Tsang et al., 1985;Ulaby et al., 1986;Chuah and Tan, 1989) is inadequate. Several empirical formulas relating the scattering and absorption coefficients to grain size and density have been proposed to solve this issue and are found in the Helsinki University of Technology model (HUT/TKK) (Pulliainen et al., 1999;Lemmetyinen et al., 2010) and in the Microwave Emission Model of Layered Snowpacks (MEMLS) . The derivation of relationships from Maxwell's equations is an attractive alternative because it is independent of particular snow conditions and the assumptions are explicit: -the strong fluctuation theory (SFT, Stogryn, 1986) allows calculations at low frequencies typically less than 20 GHz; -to cover the frequency range of 1-100 GHz, which is relevant for existing radiometers, Mätzler (1998) proposed the improved Born approximation (IBA). As in the SFT, the size of grains is given by the correlation length (Mätzler, 2002). This metric is well defined but direct measurements are only possible from 2-D or 3-D micro-structure data, which require advanced experimental techniques (Wiesmann et al., 1998). Estimations can be indirectly obtained from quantities measurable in the field like the snow specific surface area (Mätzler, 2002;Arnaud et al., 2011) or the micro-penetration profile (Löwe and van Herwijnen, 2012); -the dense media radiative transfer theory (DMRT; Tsang et al., 1985;West et al., 1993;Shih et al., 1997;Tsang et al., 2000aTsang et al., , 2007 considers snow as a collec-tion of spherical ice particles and provides analytical formulas of the effective propagation wave number, the scattering and absorption coefficients as a function of the sphere radius and the density. It is attractive because sphere radius and snow density are well defined and, when the snow grains are close to spheres like in fine grained snow (Colbeck, 1993), they are both measurable in the field. Moreover, the DMRT-ML theory has solid theoretical grounds and has been regularly improved during the last two decades.
These theories have been implemented in several models: the SFT in Surdyk and Fily (1995), the IBA in MEMLS  and the DMRT-ML in Macelloni et al. (2001); Tedesco and Kim (2006); Liang et al. (2008); Grody (2008); Brogioni et al. (2009). In addition to the underlying theory, the models differ by many aspects, including the methods for solving the radiative transfer, the presence of smooth or rough interfaces, the possible number of layers, etc. In particular, HUT uses the two-stream method and MEMLS the six-stream method while most DMRT-based models consider a larger number of streams. The comparison for a large variety of snow types (Tedesco and Kim, 2006) showed that no particular model systematically reproduces all of the experimental data. It is not yet known if the discrepancies were attributable to the electromagnetic theory, specific details of the implementations of the models, or uncertainties pertaining to input or evaluation data. Moreover, the different representations of the snow grain size in these different approaches are a major limit to the comparison. Both MEMLS and HUT/TKK models are widely used (Butt and Kelly, 2008;Durand et al., 2008;Rees et al., 2010;Brucker et al., 2011b;Gunn et al., 2011). In contrast, several groups use home-made DMRT-based models but no widely-used reference implementation exists, which limits the spread of this theory and limits the comparisons between studies.
The objective of this paper is to present the DMRT-ML (DMRT Multi-Layer) model that was released under an open source license and accompanied by a detailed documentation (http://lgge.osug.fr/ ∼ picard/dmrtml/). This model initially developed at Laboratoire de Glaciologie and Géophysique de l'Environment in Grenoble, France, for perennial snowpack (Brucker et al., 2010(Brucker et al., , 2011a and improved for seasonal snow on soil (Roy et al., 2013) and on superimposed ice (Dupont et al., 2012), is now suitable for modeling the microwave emission in a wide range of snowy continental environments. It is also designed to be extensible and uses standard For-tran90, which allows efficient computations with different types of computers and operating systems and facilitates the embedding in other models or assimilation schemes. A convenient user interface is provided by the optional bindings with the Python language. This paper is a comprehensive scientific reference for DMRT-ML users. It is structured as follows. Section 2 describes DMRT-ML in detail. Section 3 presents the sensitivity to the most important input variables and parameters and provides practical recommendations on the validity range of the input variables and parameters to use. Section 4 summarizes the results of the detailed comparisons (Brucker et al., 2011a;Roy et al., 2013;Dupont et al., 2012) with radiometric data obtained for various cold environments.

Radiative transfer equation and model architecture
The energy emanating from snowpacks is the result of the thermal emission of the snow, the substratum and the atmosphere, and of the complex propagation of this energy toward the upper snow layers. The emission and propagation can be described by the radiative transfer equation. Assuming that the medium is a stack of L plane-parallel layers containing an isotropic and homogeneous material, the equation in every layer is (Jin, 1994) By convention, matrices are written in bold, and vectors in bold and italic. I is the unit column vector. T B (z,θ,φ) is the radiance at depth z propagating along the direction with zenith angle θ and azimuth angle φ. According to the Rayleigh Jeans approximation, that is valid in the microwave range, T B (z,θ,φ) also represents the brightness temperature. Since the medium is isotropic, T B is a vector with only two non-null components (Jin, 1994, p. 19): the vertically and horizontally polarized brightness temperature. κ e and κ a are the extinction and absorption coefficients in the layer and P is the phase function. T is the physical temperature of the layer. The model solves this equation for the particular medium described by the input variables and parameters (Fig. 1). Several steps are needed: the determination of the dielectric properties of the constitutive materials (Sect. 2.2), the determination of the effective propagation constant and of the extinction and scattering coefficients for each homogeneous layer of the medium (Sect. 2.3) by using the DMRT-ML theory, the determination of the boundary conditions (Sect. 2.4) and lastly the numerical computation of the solution of the radiative transfer equation (Sect. 2.5). The result is the brightness temperature emerging from the surface in several directions. Finally, other relevant variables are accessible via additional post-computations (Sect. 2.6). For each homogeneous layer k, the input variables and parameters are grain size a k , snow density ρ k , physical temperature T k , stickiness τ k and liquid water content LWC k . The extinction and absorption coefficients κ e and κ a are calculated with the DMRT-ML theory. The boundary conditions require the input of the downwelling brightness temperature T atmos B and the variables and parameters to compute the reflection coefficient r(θ) of the underlying interface (see Table 1). The DISORT method is used to solve the radiative transfer equation in the snowpack and provides as output the brightness temperature T B (θ, z = 0) emerging from the snowpack.

Dielectric constants
The computation of the scattering and extinction coefficients in the DMRT-ML theory requires the dielectric constant of the scatterer and background materials. When snow is dry (temperature strictly below the melting point), the grains are assumed to be composed of pure ice, whose dielectric constant ε ice was measured by Mätzler and Wegmüller (1987) and Mätzler et al. (2006): where ν is the frequency in GHz, T the ice temperature in K and j 2 = −1. This model is valid for temperatures above 240 K and in the microwave range (1-200 GHz) but its accuracy certainly decreases at low frequencies (Jiang and Wu, 2004).

G. Picard et al.: Simulation of the microwave emission of multi-layered snowpacks
When the snow is wet (at 273.15 K), ice grains are coated by liquid water with an uneven distribution controlled by capillarity forces. Since such a complex particle cannot be taken explicitly into account in the DMRT-ML theory and given the low content of liquid water usually present in moist and wet snow (less than 8 %; Colbeck, 1993), DMRT-ML assumes that the grains are composed of a fictitious homogeneous material. The dielectric constant ε wet ice of this material is calculated using the mixture relationship (Borghese et al., 2007;Chopra and Reddy, 1986): where LWC is the liquid water content expressed here as the ratio between the volume of liquid water and the volume of ice present in snow (in m 3 m −3 ) and the water dielectric constant ε water is given by (4f)

Snow extinction and absorption coefficients, and phase function: the DMRT-ML theory
The extinction and absorption coefficients, as well as the form of the phase function are obtained by the so-called DMRT-ML theory. Different versions of this theory have been published over time and the underlying assumptions can significantly differ from each other. Four characteristics distinguish these versions: (i) the underlying approximation used for the DMRT-ML derivation, (ii) the particle size with respect to the wavelength, (iii) the stickiness between particles and (iv) the distribution of particle size. DMRT-ML proposes some of these versions and allows various options as detailed here. DMRT-ML uses the quasicrystalline approximation with coherent potential (QCA-CP) (Tsang et al., 2000b) and is limited to particle size smaller than the wavelength. This may be a limitation at frequencies higher than 37 GHz and with large grains commonly found in aged snow. The calculation for large particles requires a Mie-like development (Tsang et al., 2000a) and is computationally much more intensive than the QCA-CP calculation. In addition, it leads to a form of the phase function that is incompatible with the optimization of the radiative transfer equation resolution used in DMRT-ML (Sect. 2.5). To avoid the strong divergence that characterized QCA-CP with large particles, Grody (2008) proposed an empirical and computationally efficient treatment of this issue. He noticed that snow is composed of particles with a broad range of sizes, which results in a smoothing of the undulation characteristic of the Mie resonances. Hence, a good approximation of the medium scattering efficiency Q s for large particles is the asymptotic limit, i.e., Q s ≈ 2, accounting from the fact that the absorption is weak (Twomey and Bohren, 1980). If enabled by the user, DMRT-ML applies this idea by limiting the Q s maximum value to 2. Nevertheless, this ad hoc correction has not the objective to replace the rigorous solution for fine-grained studies (Tsang et al., 2000a).
Recent versions of DMRT-ML introduce the concept of stickiness. Instead of considering randomly positioned nonpenetrable spheres, the sticky spheres are attracted to one another and tend to form clusters with large voids between. This concept is meant to better represent the micro-structure of natural snow using solely spherical grains. In this case, the stickiness is also a means to account for coarse grained snow by considering that such snow is made of small clustered particles. However, DMRT-ML only implements the "short range" version of the sticky formulation, which assumes that the clusters are small with respect to the wavelength (Tsang et al., 2000b, pp. 504-505). In this simplified case, the phase function remains identical to that of the nonsticky small particle case for which the optimization of the resolution of the radiative transfer equation works (Sect. 2.5).
Another way to improve the representation of the microstructure of natural snow is to consider a collection of spheres with different sizes. In DMRT-ML, the particle sizes follow the Rayleigh distribution (Jin, 1994). In addition, only the nonsticky particle case is available because the formulation with both stickiness and size distribution leads to a quadratic system of equations that is difficult to solve and is computationally intensive (Tsang and Kong, 2001, p. 430).
In summary, DMRT-ML includes two implementations: -QCA-CP mono-disperse, with optional "short range" stickiness, and optional Grody's-based empirical treatment of large particles; -QCA-CP poly-disperse with a Rayleigh distribution, no stickiness and no large particles.
The mono-disperse version is formulated according to Shih et al. (1997). The effective dielectric constant without scattering E eff0 is obtained by solving the following quadratic Eq. (3) in Shih et al. (1997) with a = 0 or Eq. (5.3.125) in Tsang and Kong (2001): where f is the fractional volume of scatterers, b and s are the dielectric constants of the background and scatterers, respectively. The effective dielectric constant with scattering (Eq. 3 in Shih et al., 1997) combined with Eq. (5) yields: where a is the radius of the spheres and k 0 = π/λ is the wave number with λ the wavelength. t is zero for nonsticky spheres and otherwise is the largest solution of Eq. 6 in Shih et al. (1997) where τ is the stickiness parameter (Shih et al., 1997;Tsang and Kong, 2001, p. 430). At last, the extinction and absorption coefficients are respectively given by with the imaginary part indicator and κ s the scattering coefficient.
In the case of small particles, the phase function in DMRT-ML has the same form as with independent particles Rayleigh phase function (Jin, 1994, Eqs. 2-26 to 2-28).
The DMRT-ML theory inherits its name from the fact that it extends the classical radiative transfer theory, which requires the particles not to interact with each other. This condition is met only when the fractional volume f is less than a few percent (Ishimaru and Kuga, 1982). Despite its name, recent works suggest the DMRT-ML theory does not work for any large fractional volumes (see Sect. 3.3 for details). As a consequence, for icy layers or layers subject to several melt-refreeze cycles, it is recommended to use the DMRT-ML theory for a collection of air bubbles embedded in ice background. This is achieved in DMRT-ML by exchanging the dielectric constant of ice and air in Eqs. (5) and (6).

Properties of the interfaces
Once the extinction and absorption coefficients and the phase function within every layer are determined, all the variables and parameters in Eq. (1) are known and a general solution can be obtained independently for every layer. To obtain the particular solution and hence the brightness temperature emerging from the surface, it is necessary to propagate the radiation between the individual snow layers. This propagation must ensure the energy conservation and requires the reflection properties at every interface as well as any external source of energy.
For the interfaces within the snowpack and at the airsnow interface, DMRT-ML considers smooth interfaces. The roughness of this interface due to ripples, sastrugi and other dunes (Watanabe, 1978) may contribute especially at horizontal polarization and grazing incidence angles (Lacroix et al., 2009), but the existing formulations of the bi-static reflection coefficient are based on the method of moments and are therefore computationally intensive and limited to low frequencies (<18 GHz; Liang et al., 2009;Chang and Tsang, 2011). For flat interfaces, the reflection coefficient only depends on the zenith angle and the refractive indexes of both sides of the interface. In DMRT-ML, it is obtained by Fresnel's relationships (e.g. Jin, 1994, p. 59) using the effective dielectric constants calculated with Eq. (6).
At the air-snow interface, the energy coming down from the atmosphere T atmos B is an input variable given by the user. It is assumed isotropic in the current version of DRMT-ML but this assumption can be easily changed. The calculation of T atmos B from meteorological data requires an external model such as proposed by Rosenkranz (1998) and Saunders et al. (1999).
At the bottom interface, to account for the diversity of media that can be present below snow covers (e.g., soil, ice, lake ice, sea ice) and the diversity of electromagnetic modeling approaches for the soil, several substratum models are available in DMRT-ML and the addition of new models is easy. The role of the substratum model is to provide the reflection coefficient of the interface. DMRT-ML only takes into account the reflection of the coherent wave and neglects the diffuse reflections at the bottom interface. It means that the roughness of the interface is only partially taken into account. In addition to reflecting downwelling energy, the substratum is an emitter: the energy emitted and entering into the snowpack is calculated using the temperature of the substratum given as an input variable and the emissivity deduced from the reflection coefficient and Kirchhoff's law.
DMRT-ML proposes 11 models that cover soil, ice, lake and semi-infinite snowpack. The user selects the type of substratum model using an identification number and provides the required variables and parameters depending on the type of model (Table 1). Since these substratum models have been published elsewhere, their evaluation is not addressed in this paper. In the case where snow lies over ice, DMRT-ML offers two options. If the ice is free of bubbles, isothermal and semiinfinite, the best option to use is the "ice" substratum model (ID 202, or ID 1 with ice dielectric constant). In any other conditions, it is recommended to represent the underlying ice using layers with a high density (up to 917 kg m −3 for bubble-free ice) and no substratum. The number of layers to use depends on the temperature and density profile in the ice. The total depth of the ice layers must be large enough to avoid "leakage" at the bottom. Special attention is required at low frequencies (like L band, 1.4 GHz) since ice absorption can be very weak and the layers well below 100 m can significantly contribute. Table 1. Available models for the reflection coefficient of the bottom interface and required input variables and parameters. SM is the soil moisture (volume fraction), ε is the complex dielectric constant, σ is the surface root mean square height (in meters), f clay and f sand are the fractions of clay and sand, ρ orga is the density of dry organic matter (kg m −3 ), T is the temperature (in K) and Q and H are dimensionless parameters. P99 stands for Pulliainen et al. (1999), D85 for Dobson et al. (1985), M87 for Mätzler and Wegmüller (1987) and M06 for Mätzler et al. (2006

Solution of the radiative transfer equation using the DISORT method
Once the extinction and scattering coefficients of every layer and the reflection coefficients at all the interfaces are known, the radiative transfer equation is solved using the DISORT method (Chandrasekhar, 1960). This method takes into account multiple scattering within and between the layers, which is an asset with respect to the iterative method (Tsang et al., 1985;Jin, 1994;Ishimaru, 1997) for which the number of calculated order of scattering is limited. It also computes the energy propagation in an unlimited number of directions (or "streams") as opposed to the two-stream (Pulliainen et al., 1999) or six-stream  methods whose formulations are based on a fixed and small number of directions. The drawback of the DISORT method is usually the computational cost. However, in the case of passive remote sensing, isotropic media, and when the phase function has a simple analytical form, the azimuthal dependence in Eq. (1) can be analytically integrated. This simplifies the equation and reduces the numerical complexity and computation cost with respect to other cases like active remote sensing (Stamnes et al., 1988;Picard et al., 2004) or anisotropic media. For snow passive microwave modeling, the assumption of an isotropic medium is reasonable and the DMRT theory with small scatterers predicts that the phase function is the Rayleigh phase function (Jin, 1994, Eqs. 2-26 to 2-28), which only involves the cosines of the azimuth angle.
In the single layer formulation of DISORT, the integration over the zenith angle θ uses a Gaussian quadrature (Jin, 1994, p. 96), i.e., it is replaced by a discrete sum of integrand evaluation at optimal angles. The number of angles n is defined by the user. This approach cannot be seamlessly transposed for multiple layers. The variations of refractive index -mainly driven by the snow density profile -cause a change of the direction of the streams between the layers (Fig. 2). A possible approach (Liang et al., 2008) uses the same Gaussian quadrature in every layer (as in the singlelayer case), which ensures an optimal integration at the expense of complex boundary conditions since cubic spline interpolations are needed to "reconnect" the streams. We use a simpler approach in DMRT-ML similar to that proposed in Jin (1994, p. 151). The angles are determined by Gaussian quadrature only in the most refractive layer. In the other layers, the angles are determined from Snell's refraction law, which ensures the continuity of the streams between the layers. The boundary conditions are simpler at the expense of a sub-optimal integration (except in the most refractive layer). This issue is easily compensated by increasing the number of streams n.
Another consequence of using the refraction law to determine the stream angles is that the number of streams (n k , k = 1 ··· L, where L is the number of layers) varies from one layer to another. This is caused by the total reflection at large zenith angles when a stream propagates from a higher to a lower refractive index layer. Such streams in the high refractive index layer have no companion in the low index one. This also applies to the n 0 streams emerging from the snow into the air. Since snow density is much higher than air density, n 0 is usually much lower than n. The only consequence for the user is that n is not the number of emerging streams as it would be in the six-stream or two-stream methods. In practice, it is recommended to adjust n to get the wanted number of emerging streams n 0 or alternatively to increase n until the residual variations of brightness temperature are less than the wanted accuracy. Figure 2 illustrates a 4-layer snowpack (with snow density values, from top to bottom, of 50, 400, 200, 320 kg m −3 ). Only upwelling streams are represented. Fig. 2. Upwelling streams for a 4-layer snowpack with varying density. The stream directions are calculated with Snell's refraction law using the real part of the effective propagation constant E eff (Eq. 6) that mostly depends on snow density. The number of streams is 8 in the most dense layer and decreases as the density decreases. Only 4 streams emerge in the air above the snowpack. For clarity, downwelling streams are not represented.
In this particular example with extreme variations of density, the number of streams in the air is only 4 whereas it is 8 in the most refractive layer. The total reflection phenomenon is also called internal reflection in Wiesmann and Mätzler (1999) and causes energy trapping when a layer is surrounded by less dense layers.
Once the integration is replaced by the discrete sum, the differential Eq. (1) in z is written for every angle θ j (j = ±1 ··· ± n k , j > 0 and j < 0 for the upward and downward directions, respectively) in each layer.
where μ j = cos(θ j ), P is found in Jin (1994, Eq. 2-28) and I is the unit column vector. In the following, we only present the most relevant equations, intermediate calculations are given in Jin (1994, Chaps. 4 and 5). By introducing the asymmetric and symmetric brightness temperatures (T Bs (μ j ) = T B (μ j ) + T B (−μ j ) and T Ba (μ j ) = T B (μ j ) − T B (−μ j ) respectively), Eq. (9) becomes where I is the 2 × 2 identity matrix and δ jj is 1 when j = j , otherwise 0. Note that the W matrix in Jin (1994, Eq. 4-42b) is trivial in our case W = −κ e owing to the symmetry of the Rayleigh phase function. Differentiating Eq. (10a) and combining it with Eq. (10b) leads to a second-order ordinary system of differential equations. The solutions are of the form where d is the layer depth, α m = + √ m , (m = 1...2n i ) and m are the eigenvalues and T 0 Ba,m (μ j ) the eigenvectors of the 2n i × 2n i matrix whose elements are κ e μ 2 j A pp jj (j = 1...N, p = v, h and j = 1...N, p = v, h). x m and y m are 2 × 2n i unknown constants to be determined with the boundary conditions. In DMRT-ML, the eigenvalue problem is solved using LAPACK routines (Anderson et al., 1999).
The boundary conditions express the conservation of energy at every interface. For each layer k of depth d k , every direction j and polarization, there are two boundary conditions (Jin, 1994, Eqs. 5-10a and 5-10c). The boundary condition at the top interface is +T 0 Bs,m,k−1 (μ j,k ) y m,k−1 +T 0 Bs,m,k+1 (μ j,k ) x m,k+1 where r  (Table 1) and T L+1 is the substratum temperature.
The set of boundary conditions forms a linear system of N equations (N = 4 k n k ). Since the unknowns x and y in layer k are only linked to unknowns in layer k − 1 and k + 1, the system can be arranged in block-diagonal structure and can be solved using the efficient banded matrix solver in LA-PACK. The brightness temperature emerging from the snowpack is then calculated with The value from this last equation is returned to the user as the simulated top-of-snowpack brightness temperature.

Post-computation: emissivity and reflectivity
The emissivity, i.e., the coefficient measuring the departure from a black body, is often used to present modeling results in passive microwave studies instead of brightness temperature (whose value is related to the snow physical temperature). Only when the medium is strictly isothermal, that is, when the physical temperature of snow and the underlying medium is uniform and equal to T , it is possible to compute the emissivity using a single simulation and the following equation: For a nonconstant temperature profile, which is the rule for any natural snowpack, the definition of emissivity is not trivial. To mimic the in-equilibrium case and the Kirchhoff law, the emissivity can be defined as one minus the reflectivity of the medium (Peake, 1959). The calculation requires two simulations with slightly different atmospheric brightness temperatures (T atmos B and T atmos B + T atmos B ). Assuming T B and T B are the results of these simulations, the reflectivity and emissivity are given by In practice, using T atmos B = 0 and T atmos B = 1 K is recommended.

Results
This section presents the sensitivity of DMRT-ML to the most important snow properties required as inputs. It also discusses the limitations of the model and provides the range of validity of the input variables whenever possible. Figure 3a shows brightness temperatures at 18 GHz as a function of the zenith angle calculated with DMRT-ML and calculated by Kong et al. (1979) together with observations reported by Tsang and Kong (2001, Fig. 7.7.2). The medium is considered semi-infinite, the grain size was chosen to be 1.75 mm and the ice dielectric constant ε ice = 3.2 + i0.016. The result of the DMRT-ML simulation with prescribed dielectric constant (solid line) closely matches the original calculation (dotted line) by Tsang and Kong (2001) up to incidence angles of 65 • . This provides a technical validation of our implementation in the single layer case. However, the imaginary part of the ice dielectric constant used by Kong et al. (1979) was an order of magnitude higher than the one obtained with the more recent Eq. (2) from Mätzler and Wegmüller (1987). With the latter parameterization, the brightness temperatures simulated using the same grain size of 1.75 mm are much lower (dash line in Fig. 3b), leading to a strong disagreement with the observations. Modeling results and observations can be re-conciliated by using a smaller radius of 0.83 mm (solid line in Fig. 3b). This new simulation yields results very close to those of the original simulation and observations. Even if Eq. (2) is likely closer to reality than earlier formulas, the parameters of ice dielectric constants are still not  (Tsang and Kong, 2001) (dotted curve) and with experimental observations (filled dots). The medium is semi-infinite, with a density of 350 kg m −3 , a temperature of 272 K, a grain radius of 1.75 mm and an ice dielectric constant of 3.2 + j 0.016. Simulations with a more realistic dielectric constant (Mätzler and Wegmüller, 1987) are shown on the right panel with the original radius of 1.75 mm (dashed curve) and refitted radius of 0.83 mm (solid line). The gray bars represent variations of the dielectric constant imaginary part of ±20 % around in the later case.

Sensitivity to ice dielectric constant
well constrained by measurements especially at low frequencies (Jiang and Wu, 2004) and the uncertainty is unknown. The consequences on the predicted brightness temperature can be significant as illustrated in Fig. 3b (gray area) obtained by assuming an error of ±20 % of the dielectric loss (imaginary part of the dielectric constant). At a 53 • incidence angle for instance, the error is 11 and 12 K at vertical and horizontal polarizations respectively. Such an error is significant and must be taken into account in the interpretation of simulations by thermal microwave emission models.

Sensitivity to the grain size
The significant sensitivity of microwave thermal radiation to snow grain size is widely recognized (Zwally, 1977). It stems from the fact that snow grains are usually smaller than the wavelength in the microwave range and their scattering coefficient (Eq. 8c) increases as the cubic power of the sphere radius a. Figure 4 presents the variation of the scattering efficiency (Q s = 4aκ s /3f ) as a function of sphere radius. The calculations with DMRT-ML (for various densities) are plotted as solid lines.
In contrast, the absorption coefficient increases linearly with the size. It results that the emissivity and the brightness temperature of a semi-infinite medium strongly decrease with size (Fig. 5). The difference between vertical and horizontal polarizations also slightly decreases as scattering increases.
An important consequence of this cubic dependence is that a collection of spheres with different sizes is not equivalent to a collection of identical spheres with the averaged size. The contribution of the largest spheres of the collection to scattering is comparably greater than the contribution of the small-est. However, Jin (1994) shows that, under the assumption of small grains, any collection can be represented by equalradius spheres with an equivalent grain size a c . This size is always larger than the average of the distribution but depends on the distribution shape as well as snow density (Jin, 1994, Eq. 3-42 and Figs. 3-11, 3-12). In DMRT-ML, we implemented the calculation for a Rayleigh distribution of size and reached the same conclusion as Jin (1994). Unfortunately, the results are highly dependent on the choice of the distribution and especially on the slope of the upper tail of the distribution. For instance, a log-normal distribution -relevant for snow (Flanner and Zender, 2006) -features many very large grains for a reasonable mean grain size, which leads to the divergence of the DMRT-ML calculation and breaks the assumptions of small scatterers. In practice, the choice of the distribution is difficult and is related to the more general issue of the representation of snow by a collection of spherical grains. Figure 4 also illustrates the empirical correction for large particles proposed in DMRT-ML (Sect. 2.3). The corrected scattering efficiency (circles) is limited to a maximum value of 2, which corresponds to the theoretical asymptotic value for very large particles (Grody, 2008). This correction eliminates the unrealistic divergence of the scattering efficiency for large grain sizes (solid line). Nevertheless, the quality of this correction is difficult to evaluate for dense media. For a sparse medium however, the scatterers are independent and the scattering efficiency is obtained by Mie's calculation (Warren, 1982). Figure 4 (dashed blue curve) shows that the Mie scattering efficiency at 89 GHz diverges from DMRT-ML with density tending to 0 kg m −3 (i.e., independent scatterers) for grain radii larger than 0.75 mm (i.e., a ≈ λ/5 and Q s of 2-3). This is in agreement with our correction.  However, Mie efficiency reaches a maximum value of nearly 5 and remains largely above 2 in the range of realistic grain sizes. This is not reproduced by our correction. In fact, the convergence towards 2 is only observed at much larger grain sizes. This result suggests that the empirical correction would be more accurate if the efficiency limit were increased up to a value around 4. However, this result was derived in particular conditions (sparse medium, 89 GHz) and its generality is unknown. We therefore recommend to use the correction with caution and only when a very limited number of layers in the snowpack have large grains. It is worth noting that the condition Q s < 2 is valid for most types of snow at frequencies lower or equal to 89 GHz as shown in Fig. 6. Harlow and Essery (2012, Fig . 11) show using Mie-QCA with sticky spheres (τ = 0.2) that this assumption is valid up to about 60 GHz for 1 mm-radius particles and to 200 GHz for particles of 0.3 mm radius and smaller. This confirms that the Fig. 6. Range of grain sizes for which the DMRT-ML QCA-CP nonsticky is reasonably valid as a function of the frequency for snow at 260 K and 200 (red) or 300 kg m −3 (black). The upper limit is the grain size for which the scattering efficiency reaches a value of 2.
general criteria a λ/5 is adequate to evaluate the validity of QCA-CP for small scatterers.
However, even if the condition Q s < 2 is true, we note that the brightness temperatures become unrealistically low before Q s reaches a value of 2 (which occurs for a > 1 mm in Fig. 5). For example, the lowest brightness temperature ever observed at 89 GHz in Antarctica by AMSR-E is 117 K.

Influence of snow density
Snow density is involved in many components of the model. It drives (i) the proximity of the grains, thus the scattering coefficient of the medium in relation with the grain size (Eqs. 5, 6); (ii) the mass of ice, thus the absorption coefficient; (iii) the real part of the refractive index, which determines the stream angles, and the transmission and reflection coefficients of every interface (West et al., 1996).
The first two effects are illustrated in Fig. 7 at 37 GHz and for a grain size of 0.3 mm. The absorption and scattering coefficients are calculated assuming a medium of "ice spheres in air" for density less than half of the pure ice density (i.e., fractional volume of 50 %) and "air spheres in ice" otherwise. The discontinuity observed in Fig. 7 comes from the fact that both representations are not strictly equivalent even at a fractional volume of 50 %.
The absorption coefficient increases almost linearly with the density because the ice is the only absorber. In contrast, the case of the scattering coefficient is more complex. For very low density, the medium is sparse and the scattering coefficient calculated with the DMRT-ML theory increases closely to the one calculated with the independent scatterers assumption (Fig. 7, dashed line). However, for densities larger than 100 kg m −3 , the latter becomes invalid, because the scatterers are too close to be considered independent -they are in fact in the shadow of each other, which weakens their scattering efficiency (Liang et al., 2006). This effect is nicely captured by the DMRT-ML theory (Tsang et al., 2000b, Fig. 10.4.5), which predicts that the rate of variation of the scattering coefficient decreases with density. The coefficient reaches a maximum at a density around 150 kg m −3 and decreases for higher densities. Although this general behavior is expected, the accuracy of the DMRT-ML theory is not well known at such high densities. Recent work using exact numerical calculation has shown that the theory DMRT-ML QCA starts to deviate from fractional volume around 30 % (Liang et al., 2006), i.e., a density of 275 kg m −3 (validity range represented by circles in Fig. 7). However, the generality of this result for smaller grains, moderate stickiness or under the QCA-CP assumption is unexplored. If we assume that the value of 30 % is correct and applies also for "air spheres in ice", the theory would be valid in the range 640-917 kg m −3 (represented by squares in Fig. 7). Unfortunately, snow density falls in the intermediate range (275-640 kg m −3 ) in many conditions. To deal with this issue in practice, a pragmatic option is to consider the deviation above 30 % fractional volume is moderate with respect to other errors (such as grain size measurements) and to apply DMRT-ML QCA-CP using the most adequate medium representation as a function of the density as in Fig. 7. A second option is to interpolate the scattering and absorption coefficients using polynomials fitted with anchor points taken in both domains where the theory is valid. This option called "bridging" (Dierking et al., 2012) is appealing because it yields continuous relationships as a function of the density. However, the sensitivity to the choice of the polynomial order and the anchor points has to be evaluated. Therefore, the "bridging" option is not implemented in DMRT-ML yet. Figure 8 shows the scattering coefficient at 37 GHz as a function of density for different values of the stickiness parameter τ and grain radius a for a given temperature. It shows that the scattering coefficient increases as the stickiness parameter τ decreases (from blue to black to green curves). The lower values of stickiness correspond indeed to stronger attractions between spheres and a more pronounced clustering effect. Clusters are "seen" by the microwaves as larger objects than the particles they are made of. It means that they scatter more than if the particles were randomly positioned (i.e. nonsticky case). However, the stickiness parameter is not just a scaling factor of the grain size. In fact, a cluster of small particles is not equivalent to a large particle as illustrated in Fig. 8: the sensitivity to the density is different between a cluster (green curve) and a large particle (red curve). The stickiness tends to shift the maximum of the scattering coefficient toward a larger density.

Influence of the stickiness
In practice, choosing a realistic value of stickiness to represent natural snow is difficult. There is currently no means Fig. 7. Scattering (blue and plain symbols) and absorption (red and hollow symbols) coefficients at 37 GHz as a function of the density. The temperature is 260 K and the grain radius is 0.3 mm. The model "ice spheres in air" is used for densities less than half the pure ice density (458.5 kg m −3 ). For higher densities, the model "air spheres in ice" is used. Circles and squares show the domain of validity of the DMRT-ML theory for each model. The scattering coefficient under the assumption of independent "ice spheres in air" is shown as a blue dash curve. to estimate this value either from field measurements, 3-D images of snow micro-structure or snow evolution model outputs. As for choosing a grain size distribution, the core of the problem is the representation of snow by spheres.  suggest to use τ = 0.1 because it yields a frequency-dependence in better agreement with measurements. With such a low value, the clustering effect is in fact very strong and the size of the cluster approaches the wavelength (long-range effect), which explains the change of the frequency-dependence. Mätzler (1998) suggests to use a higher value, τ = 0.2, based on a comparison between the density dependence in DMRT-ML with various degrees of stickiness and in the improved Born approximation (Mätzler, 1998).
The implementation of the stickiness in DMRT-ML is limited to the "short range" version, i.e., both the grains and the cluster must be small with respect to the wavelength. In practice, τ should be larger than its theoretical minimum value of (2 − √ 2)/6 = 0.098 (Tsang et al., 2000b, p. 427). After Tsang et al. (2000b, Fig. 8.4.3) and our own calculation (not shown), the "short range" calculation is valid for τ values larger than 0.25 and grain sizes, which respects the small scatterer assumption. For lower τ , the validity depends on the grain size. Figure 9 shows the relationship between the brightness temperature at 19 GHz and horizontal polarization as a function of the total column liquid water content. The snowpack is considered homogeneous, except that the liquid water is concentrated in the top 10 cm. This calculation confirms the strong influence of the liquid water on brightness temperature, which is exploited to detect snowmelt events from passive microwave observations (e.g. Picard and Fily, 2006). It also shows that brightness temperature reaches a nearly constant value from about 0.5 kg m −2 of liquid water. This feature explains why the retrieval of the amount of liquid water from passive microwave observation is probably impossible. The threshold value, 0.5 kg m −2 , is close to value obtained with the MEMLS model (Tedesco et al., 2007) and that obtained indirectly by comparing observations and outputs of the RACMO regional snow and atmosphere model (Kuipers Munneke et al., 2012).

Layered snowpack
It is well known that the natural snowpack is composed of relatively homogeneous layers formed by episodes of accumulation and subject to metamorphisms. Hence, the density and grain size (and stickiness) can be different in every layer. Since the electromagnetic properties (i.e., the scattering and absorption coefficients and the effective constant of propagation) have a nonlinear dependence to the snow properties, the single-layer "average" snowpack is not electromagnetically equivalent to the multi-layer snowpack. In addition, the temperature is rarely uniform within the snowpack, especially close to the surface where the strong temperature gradients are caused by the daily variations of solar energy (Brandt and Warren, 1993).
Several effects can result from the layered nature of the snowpack. In general, accounting for the multi-layered nature is particularly important for -the difference between the horizontal and vertical polarizations. The difference is particularly sensitive to reflection at the air-snow interface, which is driven by the Brightness temperature (K) radius 0.5 mm radius 1.0 mm Fig. 9. Brightness temperature at 19 GHz and horizontal polarization as a function of the total amount of liquid water in the snowpack. The temperature is 273 K, density is 300 kg m −3 and there exist two grain size values (0.5 and 1 mm). The liquid water is concentrated in the top first 10 cm of the snowpack. density in the upper layer and by the reflection at the internal interface, which depend on the contrast of density between successive layers. There are several experimental and theoretical evidences of this effect (Mätzler et al., 1984;Liang et al., 2008;Champollion et al., 2013). Smoothing the profile of density in simulation results directly in a decrease of the difference between the polarizations.
-the spectrum of emissivity. It is particularly sensitive to the scattering and absorption profiles because the penetration depth highly depends on the frequency. For instance, Brucker et al. (2010) showed that the spectra observed in Antarctica cannot be explained with a homogeneous snowpack, and that the increase of grain size with depth is a necessary condition to explain the observed spectra. In some extreme cases called "anomalous spectra" by Rosenfeld and Grody (2000), the observed emissivity increases with frequency although the emissivity of a homogeneous snowpack decreases with frequency due to a stronger increase of scattering relative to that of the absorption. To illustrate the effect of the resolution of the snow parameter profiles, the calculation presented in Brucker et al. (2011a) using the original 2.5 cm high-resolution profiles of grain size and density measured at Dome C in Antarctica are presented in Fig. 10 along with calculations using lowerresolution profiles. The simulations are performed at 19 and 37 GHz with a homogeneous temperature of 219 K (the annual mean at Dome C). Only the resolution in the upper 3 m for which measurements were available is varied. Figure 10 shows that coarser-resolution results in higher brightness temperatures (up to about 8 %), except for the lowest resolution (3 m-thick layer, corresponding to the homogeneous snowpack). This complex dependence might be explained by the fact that the modeled brightness temperature is more sensitive to the resolution of the density profile than that of the grain size (simulations not shown) and that the density dependence of the electromagnetic properties is not monotonic as shown in Fig. 7. These results are specific to Dome C but they emphasize the importance of the resolution of the input parameters.
-time series of brightness temperature. Even if the grain size and density are assumed homogeneous the temperature profile is rarely uniform nor constant over time in snow and its variations can even be significant in the layer identified in the field based on the homogeneity of the grains and compactness. The variations of temperature near the surface are in general more rapid than the change of grain size and density due to metamorphism, and thus have the most important contribution to the short-term brightness temperature variations. Simulations at Dome C in Brucker et al. (2011a) show that even with grain size and density profiles taken constant over a few years (metamorphism is very slow at Dome C due to low temperature), the variations of brightness temperature are well reproduced (≈ 2 K) using measured time series of temperature profiles. Reducing the resolution of the temperature profile would result in a smoother time series.

Validation of DMRT-ML with external data
The comparison between measured brightness temperatures and results of DMRT-ML simulations using measured inputs was addressed in several studies: for a typical dry semiinfinite snowpack in Antarctica (Brucker et al., 2011a), for Arctic and sub-Arctic seasonal snowpacks (Roy et al., 2013) and for snow overlying ice as found in the ablation zone of the ice sheets (Dupont et al., 2012). In the three cases, it was necessary to estimate some parameters by optimization with respect to the measured brightness temperature. It is indeed difficult to measure all the input variables and parameters and the brightness temperatures with the same representativeness. In addition, some quantities -like the snow grain size and the soil properties -are notoriously difficult to determine in the field. The representation of snow grain in the DMRT-ML theory by spherical particles is also a conceptual difficulty. Hence, the results of the comparisons and the errors estimated in these studies depends on the methodology and are meaningless out of the context of each study. Nevertheless, these studies converge on two facts. First, the grain size derived from specific surface area measurements, i.e., optical radius, needs to be increased by a factor between 2.8 and 3.5 to be suitable as input of DMRT-ML (see discussion in   Brucker et al. (2011a) and are composed of two parts: the upper 3 m were measured with a vertical resolution of 2.5 cm (simulation marked by a symbol) and the lower part (3 to 100 m depth) is a deterministic function of the depth. The profiles at lower resolutions are generated by merging successive layers in the upper part only to emphasize the influence of measured profile resolution, the lower part remaining unchanged. The x axis reports the thickness of the layers, from 2.5 cm (original profile) to 3 m (i.e. one single layer for the upper part). Roy et al., 2013). There is no evidence that similar adjustments would be required for the density or temperature measurements. Second, the model predicts reliable dependence between horizontal and vertical polarizations near the Brewster angle (50-55 • ). Figures 11 and 12 illustrate the latter point for the polarizations at 37 and 19 GHz respectively. They show the brightness temperature at horizontal polarization versus vertical polarization for all the snow pits or pixels analyzed in Brucker et al. (2011a), Roy et al. (2013), and Dupont et al. (2012). Observations (in red) were acquired with a ground-based radiometer except at Dome C (stars), where measurements were recorded by the advanced microwave scanning radiometer for EOS (AMSR-E). DMRT-ML predictions are in blue. A large variety of environments are represented: Antarctica, ice-sheet ablation area, Arctic tundra, Arctic windy tundra, Arctic fen, and grassland. Each gray line links the observation and simulation result from the same snow pit or pixel. The length of each line illustrates the discrepancy between the observation and the simulation result of the order of 2-13 K but, as stated before, this is not an absolute error since it depends on the calibration procedure that differs between the studies. The relationship between the polarizations can be safely interpreted as it is almost independent of the calibration: Figs. 11 and 12 demonstrate that the measured brightness temperatures at both polarizations are highly correlated over a large range of about 120 K at  37 GHz and 70 K at 19 GHz. Part of this correlation stems from the linear relationship between the brightness temperature and the snow physical temperature. However the range of physical temperature -typically from −55 to 0 • C -is less than the brightness temperature ranges. It means that the emissivities at both polarizations are correlated. At 19 GHz, the correlation is lower than at 37 GHz. A probable explana-tion is the larger contribution of the substratum at the lower frequency due to the longer penetration depth. The important point is that the model nicely reproduces this general correlation. Even where the model and the observations disagree (i.e. long gray lines), the correlation between polarizations remains (i.e. the gray lines follow the general trend of the points).

Conclusions
The DMRT-ML is a physically based model used to compute brightness temperature at any frequency in the microwave range and at horizontal and vertical polarizations from input variables describing multi-layered snowpack and its environment. These variables and parameters include the profiles of snow temperature, density, grain size, stickiness, and liquid water content, the characteristics of the substratum (e.g. soil moisture, texture and temperature in the case of a soil substratum), and the downwelling atmospheric brightness temperature. The paper presents the sensitivity of the microwave emission to the most important input variables and parameters and makes recommendations on the validity ranges of these variables and parameters, either constrained by the underlying DMRT-ML theory or by the specific DMRT-ML implementation. The validation of DMRT-ML with external in situ measurements is detailed in several studies (Brucker et al., 2011a;Roy et al., 2013;Dupont et al., 2012) for various environments. The error found between predicted and observed brightness temperatures ranged between 2 and 13 K, which gives the magnitude of accessible errors but which depend on the methodology used for the comparison. In particular, the choice of the relationship to relate the measured grain size to the grain size metric relevant to the DMRT-ML theory is critical and no ideal solution exists yet. However, these studies add up to many others that have contributed to validate the DMRT-ML theory (Macelloni et al., 2001;Tsang and Kong, 2001;Tedesco et al., 2004;Grody, 2008). Currently, the most problematic point is probably the limited accuracy of the DMRT-ML theory for intermediate density values (about 300-500 kg m −3 ) that are commonly observed in natural snow and firn. An empirical correction of this problem has been recently proposed (Dierking et al., 2012). Even though it could be accurate enough, it has not the theoretical grounds that characterize the DMRT-ML theory and makes one of its merits. Further improvements are needed.
The main characteristics of the DMRT-ML implementation are the availability as an open source software, the efficiency of the computation due to the use of Fortran90 and LAPACK library, the wide range of cryospheric environments that can be modeled without any change of the source code (dry and wet snowpack, seasonal snow over soil, seasonal snow over bubbly ice, frozen lakes, perennial snow, etc.). In addition, modeling sea ice will be possible in the near future with the implementation of the dielectric constants of salted snow and water. These particularities are strong assets for the coupling of DMRT-ML with land surface models (Vionnet et al., 2012) and for the integration as an observation operator into data assimilation schemes or computationally intensive inverse methods. DMRT-ML is available for download at http://lgge.osug.fr/ ∼ picard/dmrtml/.