SMRT: an active–passive microwave radiative transfer model for snow with multiple microstructure and scattering formulations (v1.0)

. The Snow Microwave Radiative Transfer (SMRT) thermal emission and backscatter model was developed to determine uncertainties in forward modeling through intercomparison of different model ingredients. The model differs from established models by the high degree of ﬂexibility in switching between different electromagnetic theories, representations of snow microstructure, and other modules involved in various calculation steps. SMRT v1.0 includes the dense media radiative transfer theory (DMRT), the improved Born approximation (IBA), and independent Rayleigh scatterers to compute the intrinsic electromagnetic properties of a snow layer. In the case of IBA, ﬁve different formulations of the autocorrelation function to describe the snow microstructure characteristics are available, including the sticky hard sphere model, for which close equivalence between the IBA and DMRT theories has been shown here. Validation is demonstrated against established theories and models. SMRT was used to identify that several former studies conducting simulations with in situ measured snow properties are now comparable and moreover appear to be quanti-tatively nearly equivalent. This study also proves that a third parameter is needed in addition to density and speciﬁc surface area to characterize the microstructure. The paper provides a comprehensive description of the mathematical basis of SMRT and its numerical implementation in Python. Mod-ularity supports model extensions foreseen in future versions comprising other media (e.g., sea ice, frozen lakes), different scattering theories, rough surface models, or new microstructure models.


Introduction
The number and diversity of spaceborne observations from passive and active microwave sensors over snow-covered regions has considerably increased over the last 3 decades. Due to the demand for global monitoring of the cryosphere and its change, numerous algorithms have been developed to retrieve geophysical information on snow cover extent (Grody and Basist, 1996;Nghiem and Tsai, 2001), snow depth, and snow water equivalent on both land (Josberger and Mognard, 2002;Kelly and Chang, 2003;Derksen et al., 2003) and sea ice (Comiso et al., 2003;Cavalieri et al., 2012), snow accumulation on ice sheets (Abdalati and Steffen, 1998;Vaughan et al., 1999;Drinkwater et al., 2001;Winebrenner et al., 2001;Flach et al., 2005;Arthern et al., 2006;Dierking et al., 2012) wet snow (Zwally, 1977;Shi and Dozier, 1995;Abdalati and Steffen, 1997;Nagler and Rott, 2000;Steffen, 2004;Picard et al., 2007), snow temperature (Shuman et al., 1995;Schneider and Steig, 2002;Schneider et al., 2004), snow grain size (Brucker et al., 2010;Picard et al., 2012), and snow density Champollion et al., 2013). Even though many applications still rely on empirical approaches to relate snowpack properties (e.g., snow water equivalent, SWE) and measured signals, it is generally accepted that a physical understanding of the interaction between snow and electromagnetic waves is necessary to improve the accuracy and overcome inherent difficulties of the retrieval as an underdetermined problem. The retrieval of snow properties is therefore often preceded by forward modeling and data as-similation (Durand and Margulis, 2007;Picard et al., 2009;Takala et al., 2011;Toure et al., 2011;Huang et al., 2012) to predict the satellite signal from prescribed snowpack properties that can be either obtained from measurements (e.g., Rosenfeld and Grody, 2000;Brucker et al., 2011a;Rees et al., 2010;Derksen et al., 2012Derksen et al., , 2014Kontu et al., 2014) or snow models (e.g., Flach et al., 2005;Brucker et al., 2011b;Andreadis and Lettenmaier, 2012;Kang and Barros, 2012;Wójcik et al., 2008;Kontu et al., 2017). The actual modeling challenge lies in the snowpack and the underlying surface (soil, ice, or water) where the coupling of various ingredients needs to be understood with sufficient accuracy to build efficient forward models. Examples comprise scattering by snow microstructure, liquid water, salinity, ice lenses , coherent effects (Mätzler and Wegmüller, 1987;Leduc-Leballeur et al., 2015;Tan et al., 2015a), the underlying surface, and especially its roughness. All of these effects have to be taken into account by physically based snow microwave models.
Several physically based models have been developed previously mainly for passive microwave remote sensing, including HUT (Lemmetyinen et al., 2010), MEMLS , DMRT-QMS (Tsang et al., 2006;Liang et al., 2008), DMRT-ML , and other ones based on dense media radiative transfer (DMRT) (Macelloni et al., 2001;Grody, 2008;Brogioni et al., 2009). In addition, several models were tailored to low frequencies (i.e., up to a few gigahertz), such as 2S (Schwank et al., 2014), CMES (Drusch et al., 2009), WALOMIS (Leduc-Leballeur et al., 2015), and others (Tan et al., 2015a), triggered by the inception of spaceborne L-band radiometry (Barre et al., 2008). Early models for active microwave observations include only single scattering mechanisms (Bingham and Drinkwater, 2000;Flach et al., 2005;Longepe et al., 2009;Lacroix et al., 2008), which is generally sufficient at low frequencies at which scattering is weak compared to absorption. Only recently have DMRT-QMS and MEMLS been adapted to an active mode that accounts for multiple scattering Proksch et al., 2015), which is particularly relevant for high-frequency radar such as SARAL AltiKa (Verron et al., 2015). The combined active-passive capability in the same model is particularly relevant for dual-mode missions such as SMAP (Entekhabi et al., 2010). The large number of different models is a natural consequence of both the diversity of possible approaches at each stage of the calculation (e.g., effective snow permittivity, scattering, solution of the radiative transfer equation) and the wide range of applications (e.g., research versus operational use). This results in a practical difficulty of choosing the most suitable model for a given application. In addition, the scope and comparability of predictions of the same property from different models must be taken with caution, given the differences in model ingredients.
As a remedy, more and more studies include predictions from different models (e.g., Wójcik et al., 2008;Rees et al., 2010;Roy et al., 2013;Kwon et al., 2015;Sandells et al., 2017) to draw more general conclusions. Other studies directly focused on the intercomparison of different models (Tedesco and Kim, 2006;Tse et al., 2007;Tian et al., 2010;Xiong and Shi, 2013;Pan et al., 2016;Löwe and Picard, 2015;Sandells et al., 2017;Royer et al., 2017) to quantify the differences. Though insightful and necessary, these efforts did not lead to a reduction of the number of models as none of the studies considered the entirety of models and none showed a clear superiority of a single model. The latter fact was partly explained in Löwe and Picard (2015), who demonstrated the near equivalence of two approaches, namely improved Born approximation (IBA) (Mätzler, 1998) and DMRT (Tsang et al., 1985;Shih et al., 1997), which were previously considered to be different. This was achieved by relating the microstructural foundations of either approach, demonstrating the necessity to compare different microstructural formulations.
The representation of snow microstructure is critical since it immediately constrains the choice of formulation to compute the scattering coefficient. Several empirical formulations of the scattering coefficient have been developed as a function of traditional grain size (Hallikainen et al., 1987) or the exponential correlation length (Wiesmann et al., 1998). These formulations are available in the HUT and MEMLS models. But as for any empirical approach, the applicability is not guaranteed beyond the limits of calibration. This makes formulations based on fundamental principles (Maxwell equations) attractive. For instance, the DMRT theory (Tsang et al., 1985(Tsang et al., , 2000aWest et al., 1993;Shih et al., 1997) is used by several models (e.g., DMRT-ML, DMRT-QMS, Longepe et al. (2009), etc.). DMRT represents snow as a collection of ice spheres whose relative positions are constrained by the sticky hard sphere (SHS) model. Thereby a stickiness parameter controls the propensity of the spheres to stick to each other and form clusters with higher scattering power than uniformly dispersed grains. The stickiness thus has an impact on the validity of approximations when computing the scattering coefficient. Some DMRT-based models (e.g., DMRT-ML and Macelloni et al., 2001) are restricted to short-range approximation, which yields a close-form analytical solution for the scattering and absorption coefficients and the phase function. However, this approximation requires that both grain (sphere) size and the cluster size are small compared to the wavelength. While this is reasonable for snow at frequencies below 19 GHz, it is more problematic at higher frequencies (Grody, 2008). The long-range approximation relaxes the constraint on cluster size. To our knowledge, this approximation is not implemented in any available model. To additionally relax constraints on grain size, the DMRT-QCA Mie formulation is needed (Tsang et al., 2000a), allowing simulations at frequencies higher than 37-89 GHz. DMRT-QMS is the only model to implement this advanced assumption. Despite the attractive features of the DMRT theory, the representation of snow microstructure by the SHS model has a major drawback. The stickiness parameter cannot be easily retrieved from field measurements yet because microstructures of non-sticky spheres are not directly applicable to natural snow (Brucker et al., 2011b;Picard et al., 2014;Roy et al., 2013). Furthermore, estimating stickiness from highresolution microstructure images -as obtained from X-ray micro-computed tomography (µCT) -appears to be numerically unstable , leading to the conclusion that SHS is likely not a good representation for natural snow.
The IBA developed by Mätzler (1998) is an alternative approach to compute the scattering coefficient. It uses the same basic electromagnetic principles (Born approximation) as DMRT but it is not limited to a particular microstructure model. Instead of employing a particle model and characterizing their relative positions through the pair-correlation function as in DMRT, IBA uses the relative position of the ice material directly, which is mathematically captured by the autocorrelation function (ACF) of the ice indicator function (Torquato and Haslach, 2002;Löwe and Picard, 2015). In Mätzler (1998) the ACF of non-sticky overlapping spheres was investigated to obtain an analytical form for the scattering coefficient. However, in MEMLS ), the main model using IBA, the choice of ACF is limited to an exponential function that is characterized by a single parameter, the correlation length. The correlation length can be obtained from thin 2-D sections of snow samples (Wang et al., 1998;Wiesmann et al., 1998) or µCT. Even though the measurements are time-consuming, the estimation is numerically stable. On the one hand, using only a single parameter to describe the whole microstructure seems advantageous over SHS which requires two parameters, size, and stickiness. On the other hand, Mätzler (2002) had to propose different relationships between correlation length and surface area-to-volume ratio to represent different snow types, demonstrating the ambiguity of the exponential correlation length and indicating the necessity of describing snow microstructure by at least two parameters. This is also reflected by more recent attempts that use level-cut Gaussian random fields as a microstructure model for a bi-continuous medium as an alternative to the SHS model (Ding et al., 2010;Chang et al., 2014Chang et al., , 2016. This approach is very flexible, but as for SHS, the link of model parameters to natural snow microstructure and in situ measurement techniques remains to be understood (Chang et al., 2016). This requires a comparison of different microstructure models in the context of a chosen scattering theory. Due to the near equivalence of IBA and DMRT  it seems reasonable then to utilize IBA together with a library of ACFs as candidates to represent natural snow.
All examples mentioned above indicate a clear demand for a modular and extensible approach that unifies exist-ing knowledge and facilitates efficient intercomparisons of model ingredients with particular focus on the representation of microstructure. To this end we developed the Snow Microwave Radiative Transfer (SMRT) model as a versatile tool to compute backscattering and brightness temperature (active-passive mode) from multilayered media, composed of bi-continuous, random microstructures (typically snow or bubbly ice), overlying a reflective surface (typically soil, water, or ice). The originality of this new model is the flexibility for the user to select among various electromagnetic or microstructure formulations at different stages of the forward modeling problem. SMRT includes IBA, DMRT, and independent Rayleigh scattering theories to compute the scattering and absorption coefficients and the phase function. When using IBA, it is possible to choose between several representations of isotropic microstructures that are prescribed by analytical forms of the ACF. This is complemented by several soil model implementations and permittivity formulations. Additionally, language bindings are implemented to facilitate a direct comparison with widely used models (DMRT-QMS, MEMLS, and HUT) using their original code. In short, SMRT is designed to enable easy and rigorous intercomparison and exploration of electromagnetic theories, common models, and microstructure representations. SMRT version 1.0 is written in Python (https://www.python.org/, last access: 2 July 2018) and released as open source under the LGPLv3 license (https://www.gnu.org/licenses/lgpl-3.0. en.html, last access: 2 July 2018).
The paper is organized as follows. The next section gives an overview about the model architecture, the most important formulations, the code structure, and basic usage. In the third section we present an intercomparison of SMRT with other models and explore the equivalence between different microstructures. The fourth section is dedicated to the discussion of limitations and perspectives. The last section concludes the paper.

SMRT description
SMRT was designed to be easy to use and computationally efficient and to allow exploration of the various approximations or formulations available for computing snow scattering and emission in the microwave domain. Even though the goal was to maximize flexibility and versatility, some specific choices and compromises were nevertheless necessary: (i) SMRT is a radiative transfer model. This implies that interlayer interferences and coherent effects are neglected. It is not suitable for interferometric computation. (ii) SMRT considers media composed of plane-parallel, horizontally infinite, homogeneous layers and is therefore not suitable to compute 3-D effects. (iii) The current version is limited to isotropic media at the microstructure scale as well as at the scale of the snowpack. This means that microstructural anisotropy of snow is neglected (Leinss et al., 2016) and that structures formed by wind (sastrugi, dunes) are not taken into account yet. Even though SMRT is primarily designed for microwaves and snow, restrictions on spectral range and materials are not made explicit to allow for future extensions to the optical range and other random media (sea ice, layered soil, atmosphere). As a consequence of these decisions on design, the model is therefore composed of a fixed architecture, described in Sect. 2.1, and many switchable formulations described in Sect. 2.2 and 2.3 and in Appendix A.

Model architecture
The model is centered around the radiative transfer equation to compute the propagation of radiative energy in the medium produced by thermal emission in the medium (passive mode) and received from the sky (radar beam in active mode and sky thermal emission in passive mode). In addition to the radiative transfer equation, the other main components include the electromagnetic model that describes electromagnetic behavior of snow (i.e., the effective refractive index or permittivity, absorption and scattering coefficients, and phase function) and the boundary conditions between layers (called interfaces hereinafter) and at the bottom interface (called substrate hereinafter). All these components are well isolated in the code and various formulations from the literature are available. Here, only the common elements are presented; the switchable formulations are described in the following sections and appendix.
The model solves the time-independent radiative transfer equation assuming a horizontally homogeneous medium with isotropic snow at the microscopic level this is µ ∂I (µ, φ, z) ∂z = − κ e (µ, φ, z) I (µ, φ, z) (1) where I = (I V , I H , U, V ) is the reduced specific intensity defined as I = I ′ /n 2 , where I ′ is the specific intensity and n the refractive index at the same location (Mobley, 1994). P(µ, φ; µ ′ , φ ′ , z) is the 4×4 phase matrix. κ a and κ e are the absorption and extinction coefficients and the vector 1 = (1, 1, 1, 1). The extinction coefficient is given by κ e = κ s + κ a , where κ s is the scattering coefficient. Directions are defined by the cosine of the zenith angle µ and by the azimuthal angle φ. The associated solid angle is . The z axis is taken upward (as usual in Earth science), meaning that the incident beam and downwelling radiation have µ < 0, while upwelling radiation has µ > 0. This equation is valid in both active and passive modes in the microwave range. The brightness temperature T B,p , with p = H or V , is proportional to the reduced specific intensity I p = αT B,p (Rayleigh-Jeans approximation) with α = 2ν 2 k/c 2 0 , where k and c 0 are the Boltzmann constant and speed of light in the vacuum. ν is the wave frequency. In practice, for the passive mode and by using the linearity of Eq. (1), I p can be replaced by the brightness temperature and α set to 1. This is the case in our code. Further assuming that (i) the medium is azimuthally symmetric and (ii) the medium is composed of homogeneous layers ( Fig. 1), the equation becomes Here l = 1. . .L denotes the layer index ranging from the top (l = 1) to the base (l = L).
The continuity conditions at layer interfaces and the boundary condition at the bottom interface are expressed by (µ, S l,l−1 (µ))I (l−1) S l,l−1 (µ), φ, z l−1 Geosci. Model Dev., 11, 2763-2788, 2018 www.geosci-model-dev.net/11/2763/2018/ I (l) (µ > 0, φ, z l ) = R spec,bottom,(l) (µ)I (l) (−µ, φ, z l ) (4) z l is the z position of the bottom of layer l and conversely z l−1 is the height of the top of the layer l. R and T are reflectivity and transmittivity matrices. The superscript "spec" denotes the specular (a.k.a. coherent) components and "diff" is the diffuse (a.k.a incoherent) components. For a perfectly flat interface, the diffuse component is zero and the specular component is given by the Fresnel coefficients (e.g., Jin, 1994, p. 59). The "top" superscript denotes the coefficients from a layer to the one above, and "bottom" denotes coefficients to the layer below. The function S l 1 ,l 2 (µ 1 ) computes the change of beam incidence angle from layer l 1 to layer l 2 due to refraction. This function writes accordingly with the Snell-Descartes law: where n l denotes the refractive index in layer l. In case of total reflections, S l 1 ,l 2 (µ 1 ) is a purely imaginary complex number. In this case, for the sake of simplicity, we consider the transmittivity matrix, is null. Given the main governing equations (Eqs. 2, 3 and 4) it is instructive to summarize the architecture and main components of SMRT (Fig. 2). The quantities κ s , κ a , and P in the main equation (Eq. 2) are computed independently for each layer by the electromagnetic model component (Fig. 2) using one of the implemented theories (in version 1.0 IBA, DMRT, independent Rayleigh scattering) and input parameters characterizing the snow microstructure component. The interlayer reflectivity and transmittivity coefficients in Eqs. (3) and (4) are computed with the interface component (e.g., with Fresnel coefficients for flat interfaces) and with the substrate component. The effective refractive index needed for these calculations is given by the electromagnetic model component, which in turn uses material permittivity formulations of the raw materials (ice, water, air, etc.). Once fully specified, the equations are numerically solved with the radiative transfer equation solver component, which provides a numerical method adapted to the plane-parallel, multilayer configuration, and the result, that is, the intensity emerging in all or specific directions from the snowpack, is returned to the Microstructure representation Interface model (Fresnel...) Substrate model (Soil, ice, etc.) Figure 2. SMRT architecture and main components. The core components (blue) are fixed and contain no scientific code in contrast to the switchable and extensible components (orange), which define the snowpack and model configurations.
user. All formulations and methods for each component are described in the Appendix, except the IBA (one of the electromagnetic models detailed in the next section), which is essential to understand the representation of snow microstructure in SMRT.

Improved Born approximation
The implementation of the IBA in SMRT closely follows the original work of Mätzler (1998) with slight differences. The phase function in the 1-2 frame (Mätzler, 1998;Ding et al., 2010) is calculated for a two-phase medium (subscript 1 denotes the host constituent and subscript 2 denotes the scattering constituent, e.g., air and ice are used for light snow) as where the angles (ϑ, ϕ) denote the scattering direction if the incident direction is taken as the polar axis. The free-space wave number is denoted by k 0 = 2π ν/c with the wave frequency ν. The volume fraction of constituent 2 is denoted by f 2 and related to the medium density ρ by f 2 = ρ/ρ 2 . The relative permittivities of phases 1 and 2 are denoted by ǫ 1 and ǫ 2 . The temperature and frequency dependence of the permittivity is taken into account but not made explicit in the notation. Polarization information is carried in the po-larization angle χ , which is the angle between the incident electric field and scattering direction. This angle is given by sin 2 χ = 1 − sin 2 ϑcos 2 ϕ (Ishimaru, 1997, p. 21). The mean squared field ratio of field Y 2 (denoted by K 2 in Mätzler, 1998) accounts for the difference in electric field inside the scatterers and the background. This can be represented analytically for small spherical or ellipsoidal scatterers with random orientations as follows (Sihvola, 1999, Eq. 4.20): where A j are the depolarization factors along the Cartesian directions. In SMRT version 1.0, only isotropic microstructures are considered, which implies A j = 1/3. The apparent permittivity is ǫ a = 1 3 (2ǫ eff + ǫ 1 ) (Mätzler, 1998). The microstructure term M(|k d |) is a function of the difference of wave vectors in the effective medium in the incident and scattering directions, so the modulus is given by where is the scattering angle, i.e., the angle between the incident and scattering direction, and ǫ eff denotes the effective permittivity, which is by default computed with the Poldervan Santen mixing formula (Sihvola, 1999). This microstructure term can be determined from the Fourier transform C of the autocorrelation function of the medium indicator function as ) Due to the assumption of isotropy, the Fourier transform of the correlation function C(k d ) = C(|k d |) depends only on the magnitude |k d | of the scattering vector. Several analytical functions for C are implemented in SMRT, thus offering different representations of the microstructure to choose from. This is detailed in Sect. 2.3. Equations (6) to (9) fully determine the phase function in the 1-2 frame. The 4 × 4 phase matrix in the principal frame is obtained following the method of Tsang et al. (2007) and Ding et al. (2010). Co-polarization phase function matrix elements can be determined for each ϑ through calculation of p 11 = p 1-2 frame (ϑ, ϕ = π/2), and p 22 = p 1-2 frame (ϑ, ϕ = π ) and cross-polarization terms in the 1-2 frame vanish, viz. p 12 = p 21 = 0. Since the structure of the IBA phase matrix is identical to the phase matrix from Rayleigh and strong fluctuation theory (SFT) , the last two diagonal elements can be estimated as p 33 = p 44 = √ p 11 p 22 . Finally, the 4 × 4 phase matrix P in the principal frame of the radiative transfer equation (with z axis normal to the Earth surface) is obtained by rotation Mätzler et al., 2006): where α (α ′ ) is the angle of rotation from the 1-2 frame to the incident (scattering) frame. It is related to the incident and scattering zenith and azimuth angles in the principal frame by (Mätzler et al., 2006;Tsang et al., 2007, Eq. 3.23). The scattering angle is given by cos = cos θ cos θ ′ + sin θ sin θ ′ cos(φ − φ ′ ) (Mätzler et al., 2006, Eq. 3.14) so that it follows The angle α ′ is obtained by exchanging primed and nonprimed angles. Because the IBA phase matrix in the 1-2 frame is diagonal and the fourth component of the rotation matrix is orthogonal to the three others, the fourth component of the phase matrix in the main frame is also orthogonal to the three others. Except if the full Müller matrix is required by the user, the radiative transfer equation can be solved considering only the three first components, thus reducing the computational cost. This is the way it is implemented in SMRT.
The scattering coefficient κ s is, by definition, calculated from the integration of the phase matrix over all incident directions: Taking into account the isotropy of the medium , the integral can be computed in the 1-2 frame and yields a diagonal matrix with all elements equal to For the absorption coefficient, κ a is also diagonal and a multiple of the unit matrix. SMRT provides two different implementations of IBA. The first one is called "original IBA" and uses the formulation introduced in Mätzler (1998), which is used in MEMLS : where ℑ denotes the imaginary part. The second one is called just "IBA" and uses based on the Polder-van Santen mixing formula for ǫ eff . The latter is the recommended default in SMRT because the Polder-van Santen mixing formula has been shown to satisfy formal requirements (e.g., symmetry in background and inclusion permittivities) and to perform well for snow over the full range of fractional volumes (Sihvola, 1999). This allows, in particular, the representation of pure ice lenses and ice crusts in the snowpack using IBA. The effective permittivity is not only needed to compute the absorption coefficient but also implicitly to compute the boundary reflection equations (Eqs. 3 and 4) to account for the refraction (Snell's law) and the Fresnel coefficients at the interfaces. The default formulation in SMRT IBA is the Polder-van Santen mixing formula as in Mätzler (1998) and Mätzler and Wiesmann (1999). Compared to the classical Maxwell-Garnett formula, it is symmetrical between the scatterers and the background and has been shown to be slightly better for snow (Mätzler, 1996;Sihvola, 1999).

Microstructure representations
Different electromagnetic theories use different microstructure representations. In the simplest setting of Rayleigh or independent Mie scattering for a collection of spheres, the microstructure is solely characterized by the sphere radius. The positions of the scatterers are random and uncorrelated, meaning that interpenetration is possible. In DMRT the microstructure is provided in terms of the Fourier transform of the pair-correlation function (Tsang et al., 2000a) and analytical developments have been mainly given for the SHS model, which is determined by two parameters, the sphere radius and the stickiness τ . In IBA, the microstructure is provided by the ACF as shown in Sect. 2.2. Analytical expressions of ACF for independent spheres and thin shells are given in Mätzler (1998) and MEMLS proposes a generic exponential function  parametrized by the correlation length.
SMRT provides a unified and versatile vision of the microstructure representation. Any microstructure model is defined by specifying the set of required and optional parameters and by providing, at least for use with IBA, an analytical expression of ACF, either for the real-space form or its Fourier transform (or for both). Though IBA requires only the Fourier transform, see Eq. (9), some microstructure models suggested in the literature such as the level-cut Gaussian random field model (Ding et al., 2010) are rather based on real-space expressions. SMRT handles these cases using automatic Fourier transformation. Due to isotropy, required 3-D Fourier transforms can be expressed as 1-D Bessel transforms, which are numerically handled as fast (discrete) sine transforms according to Lado (1971): in terms of k d = |k d |.
Overall, the microstructure representation in SMRT closely follows a library concept as commonly employed for small angle scattering software such as in Breßler et al. (2015). In version 1.0, five different microstructure models are implemented as a starting point. Some microstructure models are defined by the Fourier transform of the ACF, and some by the real-space ACF. The most convenient characterization of a microstructure is in terms of the Fourier transform of the ACF. Presently the following models are implemented: Teubner-Strey: sticky hard spheres: in terms of the sphere volume v(a) = 4/3π a 3 , the spherical form factor P (X) defined by The SHS structure factor SF shs defined by with X = k d a and t is given by the smallest solution of the quadratic equation: under the additional condition t < (1 + 2f 2 )/(f 2 (1 − f 2 )), which guarantees SF shs (0) to be positive (Baxter, 1968;Tsang and Kong, 2001).
Note that each microstructure model comes with its own microstructure parameters. The exponential model (Eq. 19) is indeed equivalent to a real-space form C ex (r) = f 2 (1 − f 2 ) exp(−r/l ex ), which is characterized by the exponential correlation length l ex . Other models come with other parameters, which are the repeat distance d TS and the correlation length ξ TS for the Teubner-Strey (TS) model, the sphere radius a for the independent spheres (SPH) model, and sphere radius a and stickiness τ for the SHS model.
The necessity of also including models that are defined via the real-space ACF mainly originates from the use of levelcut Gaussian random field models in the context commonly termed bi-continuous DMRT (Ding et al., 2010;Chang et al., 2014). To this end we implemented a microstructure model that is defined by Gaussian random field: (26) Here r denotes the lag distance from one point to another in the medium. In the case of level-cut Gaussian random fields, the ACF of the bi-continuous medium is determined (Teubner, 1991) by the covariance C ψ (r) of an underlying zeromean, unit-variance Gaussian random field ψ from which a two-phase microstructure is obtained by "segmentation" of the continuous field values with threshold β (cut-level parameter), which is in one-to-one correspondence with the volume fraction f 2 . Our particular choice of the field correlation function C ψ in Eq. (27) was motivated by the apparent similarity to the Teubner-Strey model (Eq. 19). This particular form has been investigated by Roberts and Torquato (1999), for example, and involves the microstructure parameters d GRF and ξ GRF , similar to the TS model. However, other choices for the ACF, as used in Ding et al. (2010), for example, based on a gamma spectral density, are possible and can be implemented in the future. For running SMRT with DMRT theory, the SHS microstructure must be selected. In contrast, when using IBA, any of the above microstructure models can be selected.

Model implementation
The model implementation is highly modular to allow switching among several formulations at each stage of the computation and adding new formulations defined by users. Another feature is the extensive use of default behaviors to facilitate an easy use by beginners but still allow experts to set advanced formulations for specific investigations or sensitivity studies, for example. The code is carefully encapsulated; each "science" component (indicated by the orange color in Fig. 2 and defined in Sect. 2.1) is designed as an independent module. Table 1 summarizes the available formulations for each component in version 1.0. Additional modules contain input-output components (green boxes in Fig. 2) and core infrastructure components (blue boxes in Fig. 2). Green and blue components do not contain any science and the core component should not be modified by the users or scientific developers.
To illustrate the mode of operation of the model it is instructive to relate the instructions of a tiny but fully functioning code snippet to the model operations carried out in the background: In the code snippet first a snowpack is built (function make_snowpack) by providing the defining properties of each layer, interface, and the substrate. Layer characteristics always include density and a microstructure model to use (e.g., microstructure using exponential autocorrelation or SHS). The specification of temperature is optional, mostly relevant for the passive mode. Additional parameters depend on the selected microstructure model. For instance, the exponential function requires the exponential correlation length while SHS requires the sphere radius and stickiness. In the code example, a 100 m thick snow layer is used (i.e., to mimic a semi-infinite medium) with a density of 320 kg m −3 , correlation length of 50 µm, and temperature of 270 K. For the interfaces among snow layers, the choice is presently limited to a "flat interface", which does not require any parameter. In the future rough interfaces could be implemented. The substrate can be selected from various Ice (Mätzler and Wegmüller, 1987) Fresh water (Mätzler and Wegmüller, 1987) Wet snow (Jin, 1994) Electromagnetic DMRT-QMS (passive only) HUT models of soil, a homogeneous medium with a flat surface (e.g., bulk of isothermal ice), or a reflector with reflectivity coefficients prescribed by the user. Each model has specific parameters and all of them require temperature for the passive mode.
In the second step, the definition of the model is completed by selecting the electromagnetic theory (that computes the scattering and absorption coefficients, phase matrix, and effective permittivity) and the radiative transfer solver. As mentioned before, some electromagnetic theories are only compatible with particular microstructure models, e.g., DMRT only works with SHS and Rayleigh works with any microstructure that defines a radius but inherently considers independent spheres. For solving the radiative transfer equation, only the discrete ordinate and eigenvalue (DORT) method is currently implemented, based on Picard et al. (2004Picard et al. ( , 2013, though two-or six-flux solvers  could be implemented in the future as well. In the next step the sensor characteristics are specified (active or passive, frequencies, polarizations, etc.). For convenience, a list of predefined sensors is available (like here, AMSR-E) but sensors with arbitrary characteristics can be defined. The last step is to launch the simulation by combining the prescribed snowpack, the sensor, and the defined "model" to obtain a result (e.g., brightness temperature, backscattering coefficient, or Müller matrix). The result of this code shows a brightness temperature of 268.2 and 251.7 K at V and H polarization, respectively.
The model is implemented in Python (2.7+ and 3.4+), which makes it easy to implement switchable formulations with default and extensible behaviors. This also avoids the cumbersome step of code compilation, though at the cost of a computational overhead compared to compiled languages. To limit this drawback, the model uses common numerical libraries, such as NumPy and SciPy extensively, allowing fast and numerically accurate calculations. The code is fully documented. It also entirely uses SI units without prefix to avoid any ambiguity.
In addition, we provide different tools for convenience: (1) to facilitate convenient computation of time series or sensitivity study by a few, clear-cut lines of code the model can be run on lists of different snowpacks. (2) To foster comparisons between SMRT and other common existing models (MEMLS, DMRT-QMS, and HUT), we provide language bindings to seamlessly run these models within SMRT, which use the prescribed snowpack in SMRT and collect results as if they were produced by SMRT. This requires that the source code of these models is separately installed (they are not distributed with SMRT for licensing reasons). Note that this feature is currently limited to the passive mode.

Model validation and exploration of the microstructure
As SMRT seeks to unify formulations from other models, a natural starting point for the validation of SMRT is a model comparison (namely with DMRT-ML, DMRT-QMS and MEMLS) to assess the validity of the implementation. This is conducted in Sect. 3.1. Note that the performance of various models against observations has been extensively evaluated in the past and is not repeated here. Instead, we exploit the fact that SMRT offers the opportunity to compare different microstructure formulations within the same electromagnetic framework to investigate equivalent aspects of different microstructure models. This is addressed in Sect. 3.2. For the sake of readability, all comparisons are carried out, unless otherwise stated, for a single snowpacksensor configuration of 37 GHz in a semi-infinite medium.
The code used to build the following figures is provided as open source (see code availability section) to let the reader explore other frequencies or more complex snowpacks.

The sparse medium approximation
For a sparse medium, i.e., when density tends to zero, many formulations must show the same behavior as the independent spheres with Rayleigh or Mie theory. In SMRT, it is possible to run several combinations of microstructure and electromagnetic models as shown in Fig. 3. The results for 100 µm radius spheres show that at the origin (for f 2 → 0) the linear trend is the same for several microstructures (independent spheres, non-sticky hard spheres and sticky hard spheres) and different theories (Rayleigh, DMRT QCA-CP, IBA). These results provide a first technical validation of the SMRT implementation of several theories. However, the sparse medium approximation is valid only for very low densities in the range 10-20 kg m −3 which is unrealistic for the goal of snow modeling. It is well known that scattering in snow must be treated with dense media theories such as DMRT or IBA. The results from Fig. 3 already indicate that the influence of microstructure on deviations from the sparse medium assumption for the scattering coefficient at low densities is more severe than the electromagnetic theory. The next sections therefore consider dense media and a detailed comparison between different microstructure models.

Comparison of SMRT to DMRT-based models
We compare SMRT to results produced from original code of several DMRT variants. Figure 4 shows the angle dependence of the brightness temperature and backscattering coefficient for SMRT DMRT compared to other models for a semi-infinite medium with a sphere radius of 0.1 mm, density of 300 kg m −3 , stickiness of τ = 0.5, and temperature of 256 K. The results reveal that the closest implementation to SMRT DMRT is the model DMRT-ML . Both use exactly the same formulation for the scattering and absorption coefficient, namely DMRT QCA-CP with small, monodisperse spheres in the short-range approximation (requiring moderate stickiness, i.e., stickiness parameter should not be small). They also use a similar method to solve the radiative transfer equation, which explains the small rootmean-square difference in brightness temperature of about 0.03 K obtained at both polarizations for the angle range 0-60 • . In contrast, the comparison of SMRT to DMRT-QMS shows larger differences since the latter computes scattering by DMRT Mie QCA and implements a different connection of streams between layers in the interface conditions for solving the radiative transfer equation Liang et al., 2008). Nevertheless, the differences at both polarizations do not exceed 0.3 K RMS, which is acceptable considering the implementations are different and fully independent. We attribute this difference solely to the radiative transfer solvers because we confirmed that running an SMRT simulation with prescribed scattering and absorption coefficients and effective permittivity pre-computed from DMRT-QMS, the brightness temperature difference of 0.3 K RMS remains unchanged. In active mode (Fig. 4b), the difference is small as well, 0.65 dB RMS at HH and VV polarizations and 1.4 dB RMS at HV polarization.
The previous results were obtained for small scatterers and moderate stickiness, which is compatible with the shortrange approximation. It is therefore of interest to investigate the limits of this implementation. To this end, Fig. 5 shows two plots for brightness temperature and backscattering coefficient as a function of sphere radius and stickiness, respectively. In the first column plot, stickiness is fixed at 0.5 and in the second, the radius is set to 200 µm. DMRT-QMS is considered here as the reference because it implements DMRT QCA Mie, which has no theoretical limita-  tions on the size of the particles and on the stickiness parameter. The results show that for radii larger than 185 µm (285 µm) the error starts to exceed 1 K (5 K). Translated to surface specific area (SSA) values, this corresponds to lower bounds of 17 and 11 m 2 kg −1 , respectively, which is relatively restrictive for most snow types (Domine et al., 2007;Roy et al., 2013;Picard et al., 2014), but may still be sufficient for some applications particularly at lower frequencies.
Similarly, stickiness values lower than 1 (0.3) yield an error larger than 1 K (5 K). Even though stickiness values for natural snow are strictly unknown due to the lack of direct measurements, indirect estimates suggest that values below unity are common Roy et al., 2013). For the active mode, the DMRT QCA short range does not significantly depart from the reference simulations (less than 1 dB) but the code fails to run for a radius above 280 µm and for stickiness lower than 5. This is due to an unrealistically large scattering coefficient compared to the absorption coefficient, leading to non-real eigenvalues in the diagonalization for the DORT method.
To overcome the restrictive range of validity of the DMRT QCA short range, and considering that SMRT version 1.0 does not provide DMRT QCA in the long-range approximation, an alternative strategy is to combine IBA with the SHS microstructure model. Figure 5 shows the results, which are much closer to DMRT QMS than DMRT QCA. The difference always remains lower than 5 K for the brightness temperature and 0.5 dB for the backscattering coefficient in the explored range of input parameters. The brightness temperature becomes larger than 1 K only for radii larger than 285 µm and stickinesses lower than 0.3. The difference in backscattering coefficient does not show significant dependence on the parameters varied. This numerical result confirms the quasi-equivalence of the DMRT and IBA theories when using the same microstructure as shown theoretically by Löwe and Picard (2015). It even extends this work as only the short-range approximation was considered by Löwe and Picard (2015). Figure 6 shows the brightness temperature predicted by MEMLS along with SMRT using the original IBA and the default IBA, which computes the effective permittivity using the Polder-von Santen mixing formula. For a faithful comparison between SMRT and MEMLS, it is required to select the IBA formulation in MEMLS among the 12 available scattering formulations. In addition, MEMLS with IBA allows a choice among different grain shapes, which controls the mean squared field ratio Y 2 . As SMRT only considers spherical scatterers, MEMLS grain type must be set accordingly (grain type 2 in MEMLS code). The microstructure in SMRT is set to the exponential autocorrelation function as in MEMLS  and depends on the correlation length, which is set to 100 µm in this computation. The results show a difference of 1.2 and 1.6 K at V and H polarization, respectively, between MEMLS and SMRT with the original IBA. The cause is not the scattering and absorption coefficients, which are very close in both models (κ s = 0.2054 m −1 and κ a = 0.3092 m −1 for MEMLS and κ s = 0.2056 m −1 and κ a = 0.3087 m −1 for SMRT). Likewise, the effective permittivities are numerically close, 1.5244 in MEMLS and 1.5236. The difference is thus likely due to the different methods used to solve the radiative transfer. MEMLS uses a six-flux solver while SMRT uses the DORT method with 32 streams in the simulations presented in this paper. Similar discrepancies were observed when comparing MEMLS to DMRT-ML and DMRT-QMS ( Fig. 2 Royer et al., 2017). An implementation of the six-flux solver in SMRT would provide a route to further explore this issue. It is worth noting that setting a low number of streams  in DORT (e.g., two or six) is not recommended and is not equivalent to the two-flux and six-flux methods, which use specific stream angles and integrals of the bistatic scattering coefficient. Figure 6 also highlights the difference between original IBA and the default IBA in SMRT. The default IBA results in higher brightness temperature by 1.2 and 1.3 K on average at V and H polarization, respectively. The reason is a slightly higher absorption of κ a = 0.3426 m −1 versus κ a = 0.3092 m −1 , while all the other properties remain the same.

On the equivalence of microstructure models
Equipped with the confidence from the previous sections that SMRT is working as desired, we shall address an actual, open scientific question. Setting the correct microstructure parameters in microwave model simulations from in situ observations or snowpack simulations is notoriously difficult and nearly every study uses a different approach. To this end we demonstrate how the equivalence between different approaches can be investigated with SMRT.
The problem originates from the fact that high-level microstructural characterization in terms of the ACF is commonly not available since complete profiles of µCT or 2-D thin sections for the entire snowpack are rare. Instead, density and SSA are commonly measured or predicted by snowpack models and the initialization of microwave microstructure models relies on them. The density is unambiguous, the parameter is manifest for each microstructure model, and no problems should be expected. In contrast, using SSA is a bit more involved. Theoretically, the SSA is rigorously related to the slope of the ACF at the origin (Debye et al., 1957) and therefore parameterizes a basic size of the constituent scatterers. For microstructures comprising spheres, the SSA (m 2 kg −1 ) can be directly converted to sphere radius using a = 3/(ρ ice SSA). For an exponential ACF there is also the well-known relationship (Debye et al., 1957) l ex = 4(1 − f 2 )/(ρ ice SSA), henceforth termed the Debye relation. However most microstructure representations involve three parameters (all except the exponential autocorrelation function) and the additional parameters must be set as well. Although grain type is often observed in the field, quantitative relationships with the microstructure metrics (stickiness or autocorrelation function) have not yet been established and we do not consider this information here.
These issues have been solved in different ways in literature. For the SHS microstructure, Liang et al. (2008) suggest setting the stickiness parameter to "0.1 because it yields 2.8 for the frequency dependence of the extinction coefficient which corresponds to the experimental values" (Hallikainen et al., 1987). These experimental values are the basis of the extinction formulation in the HUT model . However setting stickiness to 0.1 is insufficient to strictly determine the power dependence as it also depends on the grain size and density, i.e., very small scatterers always show a dependence of a power of 4 (Rayleigh scattering). Another approach was elaborated in a series of empirical studies (Brucker et al., 2010;Picard et al., 2014;Roy et al., 2013Roy et al., , 2016Dupont et al., 2013). It consists of using non-sticky spheres (i.e., infinite stickiness parameter) and scaling the radius a computed from SSA by an empirical factor φ SHS (called "grain size scaling factor"). This factor is obtained by fitting model results to microwave observations. To prevent over-fitting, a single φ SHS was applied to all SSA measurements and the fit was performed using microwave observations at several frequencies, polarizations, and/or angles.
To explore if this latter approach is equivalent to choose an optimal stickiness value, we use SMRT to find the equivalent microstructure representations for non-sticky spheres with grain size scaling and sticky spheres. In the following, equivalent microstructures are interpreted as microstructures with the same density but different size parameters that produce the same electromagnetic behavior. This is exemplified by using SMRT IBA and matching brightness temperatures at V polarization and 55 • close to the Brewster angle to inte- grate properties of scattering and absorption coefficients and phase function (see also Veysoglu and Kong, 1996). Figure 7 shows the grain size scaling factor of non-sticky hard spheres as a function of the stickiness value to obtain this equivalence. For instance φ SHS = 2.1 (used by Picard et al., 2014) is equivalent to a stickiness value of around 0.13. Higher values of φ SHS up to 3.5 were used in the other studies (Brucker et al., 2010;Roy et al., 2013Roy et al., , 2016Dupont et al., 2013), corresponding to lower stickiness values approaching 0.1 as suggested by Liang et al. (2008). This confirms that despite using different approaches, these studies converge towards stickiness values in the range 0.1-0.2, in agreement with Löwe and Picard (2015), who retrieved the stickiness from µCT of snow samples. However, the relationship between stickiness and grain size factor depends on density, especially for φ SHS > 2.5 (Fig. 7), and thus the approach of scaling grain size cannot be strictly equivalent to selecting an optimal stickiness value.
Though the approach of using a stickiness close to 0.1 seems more physical compared to an empirical scaling factor, it also has weaknesses. Natural snow is composed of grains with variable size, which more resembles a collection of spheres with a distribution of radii (i.e., polydisperse spheres). Such dispersion is important and generally leads to increased scattering compared to the medium with monodisperse spheres with the mean radius of the polydisperse spheres (Tsang and Kong, 1992). However, the analytical treatment of the ACF for polydisperse SHSs is tedious (Gazzillo et al., 2006) and choosing the distribution form and its parameters is an open question. In the case of non-sticky small scatterers, Jin (1994) showed that a poly-disperse microstructure can be equivalent to a monodisperse sphere assembly with an effective radius. This effective radius was found to be about 1.4 times the radius derived from SSA when a Rayleigh distribution of sizes was taken (Jin, 1994). This factor would be slightly different for another distribution but this gives an order of magnitude of the size distribution effect. Based on this, Roy et al. (2013) proposed a pragmatic approach mixing the scaling approach and a fixed stickiness value. For this, they suggested using φ SHS = 1.4 found by Jin (1994) and optimizing the stickiness to obtain good fit with observations. This proposition has not been evaluated in other studies.
The exponential autocorrelation is a different and attractive solution because it involves only two parameters that should be fully determined by density and SSA. However, in practice a "hidden" third parameter must be introduced to empirically scale the correlation length in the Debye relation (Mätzler, 2002). Based on comparisons between simulations and observations, Mätzler (2002) suggested a scaling factor of 0.75 in the Debye relation and justified this adaptation with the necessity of fitting the exponential function to the real nature of snow, i.e., to the actual ACF of snow. However, more recently Montpetit et al. (2013) performed an optimization of the simulations with MEMLS on a large set of observations on the Arctic snowpack and found a different coefficient of 1.3. While the origin of this large discrepancy can be understood from the effect of shape (or equivalently size dispersity) of the 3-D microstructure (Krol and Löwe, 2016) it remains a practical problem, similar to the freedom of choosing an appropriate stickiness value. To this end we explore the connection between the Debye scaling factor and stickiness, or in other words, the equivalence between the exponential ACF with scaled correlation length and SHS. Figure 8 shows the scaling factor φ exp of the correlation length in the Debye relationship required to obtain the same electromagnetic behavior as with SHSs. Each curve is obtained, for a given density, by optimizing φ exp to obtain equivalence between exponential and SHS microstructure. The results show that stickiness higher than 0.2 corresponds to φ exp lower than 0.5, with little dependence on density. This range seems inadequate however for snow considering the values of stickiness and φ exp used in the literature. Conversely, the value of φ exp = 0.75 corresponds to a stickiness of 0.13 at 300 kg m −3 and lower at higher densities. This means that scaling the correlation length proposed by Mätzler (2002) is equivalent to stickiness values suggested by various studies (Liang et al., 2008;Roy et al., 2013). In contrast, φ exp = 1.3 found by Montpetit et al. (2013) is barely accessible for the scaled correlation length derived from the Debye relation, indicating the limitations of the exponential ACF for snow. Moreover, the large dependency on density indicates that a strict equivalence between SHS and an unscaled exponential model is not possible.
The numerical experiments facilitated by SMRT from this section show how different studies, which were hitherto not (1 − ρ/ρ ice )a to yield the same brightness temperature at V polarization as simulations with sticky hard spheres with a radius a = 100 µm. Other parameters are similar to Fig. 4. amenable to a comparison due to apparently different approaches, are now comparable and can be shown to be nearly equivalent for particular parameter choices. Moreover the results unambiguously show that density and SSA are not sufficient to appropriately characterize snow microstructure for microwave modeling purposes and that the sensitivity to a third parameter is highly significant. Until alternative measurement techniques or progress in modeling the microstructure evolution are available, the initialization of microstructure models relies on µCT characterization or some empiricism to infer the missing parameter.

Limitations and perspectives
SMRT version 1.0 bears some limitations that are inherent to the architecture as discussed in Sect. 2.1; others are related to the current set of available modules and their approximation as shown in Table 1. Some limitations could be simply overcome by implementing new modules or formulations. This section focuses on the latter category.
The scope of SMRT is currently limited to a snowpack over a surface (called substrate), which is a common approach for some applications such as soils, but may be inappropriate for other snow-covered environments in which volume scattering, layering within the substrate, or temperature heterogeneity may be important. For instance snowcovered sea ice or frozen lakes need to account for bubbly and salty ice with a nonuniform temperature profile. While the generic plane-parallel layered structure in SMRT and the DORT solver are readily suited for this kind of modeling, the electromagnetic behavior of these materials needs to be additionally implemented, which is technically easy due to the modular architecture. Bubbly ice  has been modeled with DMRT for fresh ice. This should also work for salted ice unless the scattering of brine becomes significant.
Considering soil as a volume scattering medium or accounting for inhomogeneous temperatures or wetness can be treated within DMRT and layered radiative transfer (Lu et al., 2009). Though promising, this approach has still been hardly explored. Likewise, the atmosphere could benefit from a multilayer representation as employed in specific, atmospheric radiative transfer models (Eriksson et al., 2011). Implementing atmospheric layers in SMRT would be of interest to deal with cases of strong surface-atmosphere coupling as observed around 60 GHz near the oxygen absorption band. A simple non-scattering bulk atmosphere can be prescribed in the current SMRT version; however this requires the down-and upwelling brightness temperatures and transmittance to be calculated externally.
Accurate simulations of snow on the ground in active mode would require more advanced surface scattering models than implemented in the current version. SMRT inherits from the soil modules implemented in DMRT-ML and previously in HUT and MEMLS, which were tailored to the passive mode. These modules mainly compute a specular reflection while a faithful backscatter computation is required for the active mode. DMRT-QMS includes an advanced rough surface treatment from independent numerical simulations (Zhou et al., 2004). In SMRT soil backscatter is prescribed in the current version, but an implementation of a numerical approximate method for rough soil surfaces such as the advanced integral equation method (AIEM; Chen et al., 2003) is foreseen in the future. Likewise, taking the roughness of the snow surface and internal snow interfaces into account is another interesting perspective (Liang et al., 2009).
A strong assumption in SMRT version 1.0 is the isotropy of the microstructure. Some types of snow have been shown to be highly anisotropic, especially due to differences between the vertical and horizontal directions (Löwe et al., 2013). This results in polarization effects in the volume (Leinss et al., 2016). Implementing anisotropic microstructures is possible in the existing architecture but requires significant developments at several locations, namely (i) the effective permittivity tensor (ii) scattering and absorption coefficients and phase function and (iii) solution of the radiative transfer equation taking into account the ordinary and extraordinary streams. Another, related assumption in the current version is the isotropy at the snowpack scale. Accounting for anisotropically reflecting interfaces would only require an improvement of the radiative transfer solver and the implementation of anisotropic surface reflections. However, to include all emergent effects (such as multiple scattering between surface and volume) a full 3-D model is required, which is not compatible with the SMRT architecture.
Some limitations of SMRT are inherent to the radiative transfer equation, which does not keep track of the absolute phase. This obviously prevents interferometric calculations and may be restrictive when the layer thickness is smaller than the wavelength of the microwaves, that is, at low frequencies (e.g., at L band; Tan et al., 2015a;Leduc-Leballeur et al., 2015) or in the case of thin ice lenses in the snowpack (Mätzler and Wegmüller, 1987). In some cases, ad hoc corrections of the radiative transfer solution can be implemented. For instance MEMLS  computes the effect of interferences between the interfaces of sub-wavelength layers for short phase differences at the condition that scattering is negligible and these thin layers are surrounded by thick layers. This correction is suitable for isolated ice lenses  but not sufficient for low frequencies. Another important case concerns the active mode in the backscatter direction -which is the most common configuration for radars. In such a configuration, some of the many possible trajectories of radiation propagation are paired, cyclical double bounces involving reflections between one of the interfaces (soil or air-snow surface) and the volume. Theses pairs constructively interfere with each other, according to wave theory. As a result, the backscatter for these bounces is increased by 3 dB compared to what is predicted by the incoherent radiative transfer theory. This phenomenon called backscattering enhancement has recently been taken into account by developing a specific solver of the radiative transfer equation able to distinguish the noncyclical and cyclical trajectories, and to apply a correction of 3 dB to the latter group (Tan et al., 2015b).
Another limitation concerns simulations of altimetric signal or frequency-modulated continuous-wave radar. The radiative transfer equation solver available in SMRT version 1.0 considers the stationary radiative transfer equation (Eq. 1), which is insufficient to simulate altimetry waveform or time-resolved radar echo. However, the SMRT architecture could accommodate such an enhancement with little change; only an adequate solver needs to be added (e.g., Lacroix et al., 2008).
Finally we acknowledge that the Python implementation of SMRT bears some peculiarities. By extensively using Python dynamic capabilities, the model computation is probably less efficient than specialized code, even though numerically critical code is delegated to optimized libraries through SciPy. Because of Python, the model may be inadequate for high-performance computation. In this case SMRT may still be useful for prototyping and determining the optimal subset of formulations that could then be implemented in compiled language since a numerical reference greatly helps to achieve such an optimization step. Moreover, it is worth noting that the Python ecosystem for high-performance computing is fast improving and that SMRT code may be parallelized in the future.
A new radiative transfer model to simulate emission or radar echo from a snowpack has been presented in this paper. It is built around the radiative transfer equation and specifically tailored to model snow but in the future also other planeparallel media in the cryosphere. SMRT differs from other models in its scope in many aspects. SMRT is not a new model with a more advanced theory, it is rather a repository of established formulations or widely used model configurations that can be easily interchanged. The novelty is thus to allow testing of different existing configurations and exploration of new ones, in particular regarding the microstructure. Using SMRT, we have highlighted the equivalence between different widely used microstructure representations (SHS and exponential autocorrelation function) and different approaches proposed in the literature to run simulations based on in situ measurements. These results show that to fully describe snow in microwave models requires at least three main metrics, the density, grain size, and another parameter characterizing larger-scale structural correlations of the ice matrix. The fact that these latter properties are presently inaccessible by other measurements or snowpack modeling contributes to the uncertainties in microwave simulations, and actually constitutes one of biggest challenges to solve.
The numerical validation of SMRT has shown the numerical equivalence with DMRT-ML for the DMRT QCA-CP electromagnetic formulation and has shown close results with DMRT-QMS under DMRT QCA under the small scatterer assumption in passive and active modes even though small differences remain unexplained. Larger differences are observed with respect to MEMLS, which we attribute to the six-flux method used by MEMLS to solve the radiative transfer equation. Regarding HUT, SMRT contains no sufficiently similar configuration to perform a validation. Nevertheless the language binding to the HUT code has been included for future comparisons with other configurations. Not all SMRT configurations and available microstructure representations have been tested in this study because of the large number of possible combinations; this is left to future work.
Several limitations of SMRT version 1.0 have been outlined that can be readily overcome by model extensions which are supported by modularity. The developed code is highly structured for each step of the radiative transfer calculation. The model is designed to facilitate future developments of existing and new formulations without changing existing code, which should foster community-based contributions and consolidate SMRT as a repository of the community knowledge. Future work includes implementation of new features to account for different media (e.g., sea ice), variants of electromagnetic models (e.g., DMRT QCA long range) or radiative transfer solver (e.g., six-flux solver or time-resolved radiative transfer equation) to increase the scope of applications. In this paper we focused on two widely used microstructure representations; SMRT already includes other representations and new ones, such as empirical autocorrelation functions derived from µCT, could be included, which opens a new promising way to characterize the microstructure.
Code availability. The code to generate the figures is released under the LGPLv3 open-source library and is available at https://github.com/smrt-model/smrt (Kluyver et al., 2016). The precise version used for this article is registered as https://doi.org/10.5281/zenodo.1249698.

Appendix A: Discrete ordinate radiative transfer method
The discrete ordinate and eigenvalue method (DORT) is a widely used method to solve the radiative transfer equation for multilayered media. It is particularly recommended when optical depth is thick and multiple scattering is significant (Liu et al., 2016). Overall, it is not the most computationally efficient, but it is robust. Many variants have been proposed in the literature for the scalar or vector (i.e., polarized) radiative transfer equations, for dense or sparse media, for passive mode only, and with different basis functions for the azimuthal angle (Liu et al., 2016). DORT in SMRT inherits from the passive-only variant used in DMRT-ML (Jin, 1994;Picard et al., 2013) and the active-only sparse medium variant developed in Picard et al. (2004) for the simulation of forest backscatter. It is also similar to the DORT method implemented in DMRT-QMS and only differs in the way the interface conditions are handled. In SMRT the streams in the different layers are directly connected; that is, their zenith angles are governed by the Snell refraction law, whereas DMRT-QMS uses the same angles for all the layers and spline interpolation to connect the streams.
The transformation of the radiative transfer equation into a system of linear ordinary differential equations requires the discretization of the azimuthal and zenith angular dependences. The φ dependence is treated by decomposition into cosine and sine functions: I c,m (µ, z) cos(mφ) + I s,m (µ, z) sin(mφ) with Fourier coefficients I c,m and I s,m . Because of the azimuthal symmetry of the medium, the elements of the intensity vector are either purely even or odd functions of φ so that and Similarly the phase matrix writes and because of the azimuthal symmetry of the medium, some elements of the phase matrix components are even functions of φ so that (2), multiplying by cos(mφ) and integrating over φ from 0 to 2π, we obtain and µ dI s,m (µ, z) dz = − κ (l) e (µ)I s,m (µ, z) (A8) +P c,(l),m (µ, µ ′ )I s,m (µ ′ , z) where δ m is 1 for m = 0 and 0 otherwise. We then introduce the even intensity and phase matrix: + δ m κ (l) a (µ)T (l) 1. In the following the superscript e is dropped.
The θ dependence is replaced by a set of discrete ordinates (samples at fixed angles), so that the integral over µ ′ is transformed into a discrete sum as i−1 )/2. It is worth noting that the number of samples N(l) varies among layers and is always lower than or equal to the number of streams in the most refringent layer.
Introducing the matrix-vector notation for the linear solver, we define an intensity vector I m containing the four Stokes components (or a subset of them) for all directions, that is, with cosine µ ranging from µ where This is a nonhomogeneous system of ordinary first-order differential equations with constant coefficients.
The general solution of the discrete equation (Eq. A13) in layer l can then be written where E (l),m is the matrix containing the eigenvectors E . .e −β (l),m 4N (l) z ) is the diagonal matrix describing the transmission of eigenvectors through the layer, and x (l),m are the unknown constants to be determined from the boundary conditions. For the purpose of the boundary conditions, we distinguish the upwelling (µ > 0) and downwelling (µ < 0) subspaces. In addition, for the numerical stability, the reference height is chosen at the top of the layer for the downwelling z

A2 The boundary conditions
At this point, the problem consists of determining the 8N(l) unknown constants x (l),m per layer, i.e., 4 L l=1 N(l) unknown constants in total. The boundary condition stems from the necessary energy conservation at the interfaces and depends on the bistatic reflection coefficient at the interface (for both upwelling and downwelling waves). The precise form of boundary conditions depends on the choice of cosines µ (l) i in the layers. Here, we make a specific choice following the strategy used in DMRT-ML and described in Picard et al. (2013). The cosines in the layer with the highest refractive index (most refringent) are taken to follow the (optimal) Gaussian quadrature rules, and the cosine in the other layers is deduced based on Snell's refraction law. This leads to a one-to-one connection among the streams in adjacent layers, except in the case of total reflection for which the more grazing streams in the more refringent layers are not connected to any stream in their less refringent neighbor layers . This choice yields a simple form of the boundary conditions. The downwelling waves at the top interface of a layer l are the sum of the reflected upwelling waves of layer l and the downwelling transmitted waves of layer l − 1: considering that l = 0 designates the air above the snowpack. In active mode, the downwelling beam with incidence angle θ 0 is represented by where i = 1. . .N (l) and the bracket shall be interpreted as taking the coordinates of the vector. Since θ 0 is in general not exactly on a stream angle, we perform a linear interpolation between the two nearest streams as follows: i 0 is the (integer) index so that cos θ 0 is between µ i 0 and µ i 0 −1 ; then α = (cos θ 0 − µ i 0 )/(µ i 0 − µ i 0 −1 ). In active mode, except if the fourth Stokes component of the incident beam is nonzero, this component remains null under the isotropic assumption used here. As a consequence only the first three components need to be computed in practice. In the passive microwave mode the situation is even simpler. The incident energy comes from the atmosphere (downwelling atmosphere emission), which is considered azimuthally symmetrical, so that only the m = 0 mode is nonzero and as both the forcing (incident radiative) and the source (emission) have an azimuthal symmetry, the m > 0 mode equations all have a zero solution. As a consequence for passive microwave, only the m = 0 mode needs to be solved and only the first two components of the Stokes vector are nonzero. The SMRT DORT code is built to accommodate a variable number of Stokes components which allow optimized computations. All boundary conditions provide 4 L l=1 N(l) equations linking the unknowns x (l),m for each of the L layers since the boundary conditions only connect successive layers; the system of equations takes the form of an almost block diagonal matrix. Picard et al. (2014) uses an almost block diagonal algorithm to solve it, whereas DMRT-ML and SMRT cast it as a band diagonal system for which efficient algorithms exist in LAPACK and SciPy. Solving this system yields x (l),m for each layer l and each mode m. The outgoing intensity from the top layer l = 1 can be calculated by inserting x (1)m into the general solution given by Eq. (A16) and using the topmost boundary conditions (Eq. A18): The last operation is to reconstruct the azimuth series, as follows: I (0)m (z 0 ) µ > 0 cos(mφ).
For the passive microwave case, only the first mode m = 0 is non-null. For the active microwave case, this equation gives the total outgoing intensity from the snowpack, which accounts for both diffuse and the reflected coherent radiations. The backscatter only includes diffuse radiation and in principle the coherent radiation is only in the forward direction. However, the truncation of the φ series at mode m contains spurious remainders of the coherent intensity in any φ direction. To overcome this numerical problem the total intensity (Eq. A26) is written as diffuse and coherent components (Ishimaru, 1997): and the diffuse part is obtained by subtracting the total intensity from the coherent part. The former is the result of Eq. (A26) while the latter is computed using the equations as for the total intensity except all sources of scattering are switched off (Picard et al., 2004). This includes volume scattering (κ s = 0 and P (l),m = 0) and diffuse reflection (R diff,top/bottom,(l),m = 0). Denoting with a superscript [c] the unknowns and eigenvectors obtained by this calculation, the diffuse intensity is then given by The intensity in the backscatter direction is finally obtained by setting φ = π and by linear interpolation of the intensities I [d] (φ) at µ i to the exact incident angle.
(1 − f 2 ) 4 (1 + 2f 2 − tf 2 (1 − f 2 )) 2 , where ℑ denotes the imaginary part of a complex number. The effective permittivity in the QCA approximation is given by and the extinction and scattering coefficients with and κ s = 2 9 k 4 0 a 3 f 2 In the short-range approximation, the phase matrix of DMRT QCA and DMRT QCA-CP is the same as for independent Rayleigh scatterers.