Microstructure representation of snow in coupled snowpack and microwave emission models

Abstract

This is the first study to encompass a wide range of coupled snow evolution and microwave emission models in a common modelling framework in order to generalise the link between snowpack microstructure predicted by the snow evolution models and microstructure required to reproduce observations of brightness temperature as simulated by snow emission models. Brightness temperatures at 18.7 and 36.5 $\mathrm{GHz}$ were simulated by 1323 ensemble members, formed from 63 Jules Investigation Model snowpack simulations, three microstructure evolution functions, and seven microwave emission model configurations. Two years of meteorological data from the Sodankylä Arctic Research Centre, Finland, were used to drive the model over the 2011–2012 and 2012–2013 winter periods. Comparisons between simulated snow grain diameters and field measurements with an IceCube instrument showed that the evolution functions from SNTHERM simulated snow grain diameters that were too large (mean error 0.12 to 0.16 $\mathrm{mm}$), whereas MOSES and SNICAR microstructure evolution functions simulated grain diameters that were too small (mean error$-$0.16 to $-$0.24 $\mathrm{mm}$ for MOSES and $-$0.14 to $-$0.18 $\mathrm{mm}$ for SNICAR). No model (HUT, MEMLS, or DMRT-ML) provided a consistently good fit across all frequencies and polarisations. The smallest absolute values of mean bias in brightness temperature over a season for a particular frequency and polarisation ranged from 0.7 to 6.9 $\mathrm{K}$.

Optimal scaling factors for the snow microstructure were presented to compare compatibility between snowpack model microstructure and emission model microstructure. Scale factors ranged between 0.3 for the SNTHERM–empirical MEMLS model combination (2011–2012) and 3.3 for DMRT-ML in conjunction with MOSES microstructure (2012–2013). Differences in scale factors between microstructure models were generally greater than the differences between microwave emission models, suggesting that more accurate simulations in coupled snowpack–microwave model systems will be achieved primarily through improvements in the snowpack microstructure representation, followed by improvements in the emission models. Other snowpack parameterisations in the snowpack model, mainly densification, led to a mean brightness temperature difference of 11 $\mathrm{K}$ at 36.5 $\mathrm{GHz}$ H-pol and 18 $\mathrm{K}$ at V-pol when the Jules Investigation Model ensemble was applied to the MOSES microstructure and empirical MEMLS emission model for the 2011–2012 season. The impact of snowpack parameterisation increases as the microwave scattering increases. Consistency between snowpack microstructure and microwave emission models, and the choice of snowpack densification algorithms should be considered in the design of snow mass retrieval systems and microwave data assimilation systems.

Introduction

Global observations of the snow cover extent from optical and microwave satellite observations combined with in situ data have shown a reduction in the spring snow cover . Observed decline in snow cover extent during 2008–2011 exceeded that predicted by climate models . Observations also indicate that duration of snow cover is also reducing, but they cannot determine whether mass or volume of snow has changed.

Microwave, altimetry, or coarser-scale gravity satellite sensors offer the only feasible way to measure snow mass or depth on a global scale, with microwave observations spanning the longest timescale of these. However, microwave algorithms such as those developed by and can result in large errors because of the high sensitivity of applied forward models to parameterization of the snow microstructure . In particular, the assumption of a fixed snow scatterer radius in the algorithm does not reflect the naturally changing snowpack structure. Errors in snow mass products derived from these algorithms mean that the products are difficult to use for evaluation of snow mass in climate models and unsuitable for assimilation into land surface models for streamflow forecasts . Development of the assimilation-based technique in GlobSnow allows changes in the snow microstructure to be taken into account through inversion of ground-based observations of snow depth and coinciding microwave brightness temperatures . Although more accurate than other global products, some errors remain, and the GlobSnow accuracy relies on the proximity and representativity of the ground stations . In addition, the intermediate retrieval of the snow “grain size” in GlobSnow is a parameter that also incorporates other land surface features, so is not a true representation of the snow effective diameter .

Snowpack evolution models offer a way to estimate temporal changes in snow microstructural parameters and stratigraphy e.g.. Intercomparison studies have shown large differences between snow evolution models driven by the same forcing data . Given that the mass inputs were the same for the 33 snow models considered in the SNOWMIP2 study of , it is differences in the internal snow physics and model structure (layering assumptions) that result in the wide range of simulated depth and snow mass. Temperature, temperature gradient, and density drive changes in the snow microstructure e.g., so it is likely that different snow physics assumptions in a coupled snowpack and emission model result in different thermal structures, microstructure parameterisations, and ultimately different microwave extinction behaviour.

Theoretical differences between specific electromagnetic models have been examined in , , and other intercomparisons carried out by . These studies are useful for interpreting differences in electromagnetic model outputs for a snapshot profile of the snowpack properties. Given the dependence of microwave scattering on snow microstructure, a satellite retrieval system needs some quantification of microstructure. Snowpack evolution modelling offers a means to quantify the metamorphic changes in snow microstructure. Indeed, snowpack evolution models have been coupled with microwave emission models to demonstrate the potential of this approach for snow remote sensing applications . These studies all examined the accuracy of a single snowpack model coupled with a single microwave emission model.

The purpose of this study is to inform future design of retrieval and assimilation systems where snowpack evolution models may be used to provide microstructural parameters for microwave emission models, by examining how particular snowpack and emission model choices lead to a variation in simulated brightness temperatures throughout the winter period, and evaluate how the simulated values compare to observations. The Jules Investigation Model JIM; has been coupled with three widely used microwave emission models: the Dense Media Radiative Transfer Multi-Layer model DMRT-ML;, the Microwave Emission Model of Multi-Layer Snow MEMLS;, and the Helsinki University of Technology (HUT) multi-layer model . Snowpack microstructure metamorphism is represented here by three different options with differing complexity for grain diameter evolution (or equivalently the specific surface area, SSA). These models are the grain growth models of SNTHERM SNT;, SNICAR SNI;, and MOSES MOS;. This allowed quantification of the seasonal variation in uncertainty in brightness temperature simulations from 1323 coupled snowpack–emission model systems, as evaluated against ground-based observations of brightness temperature.

The study approach, model descriptions, and field measurements are given in Sect. . Comparisons between simulations and between simulations and observations are presented in Sect. , and the implications for future approaches to the remote sensing of snow mass are discussed in Sect. .

Equations for options governing the representation of processes in the JIM model subset.

Option	Description	Model parameterization
Compaction: 0	Physical	$\frac{\mathrm{1}}{{\mathit{\rho }}_{\mathrm{s}}}\frac{\mathrm{d}{\mathit{\rho }}_{\mathrm{s}}}{\mathrm{d}t}=\frac{M_{\mathrm{s}}g}{\mathit{\eta }}+c_{\mathrm{1}}\exp [-c_{\mathrm{2}}\left(T_{\mathrm{m}},-,T_{\mathrm{s}}\right)-c_{\mathrm{3}}\max \left(\mathrm{0},,,{\mathit{\rho }}_{\mathrm{s}},-,{\mathit{\rho }}_{\mathrm{0}}\right)]$
1	Empirical	${\mathit{\rho }}_{\mathrm{s}}\left(t,+,\mathit{\delta },t\right)={\mathit{\rho }}_{\max }+\left[{\mathit{\rho }}_{\mathrm{s}}\right(t)-{\mathit{\rho }}_{\max }]\exp \left(-,\mathit{\delta },t,/,{\mathit{\tau }}_{\mathit{\rho }}\right)$
2	Constant	${\mathit{\rho }}_{\mathrm{s}}=250$ $\mathrm{kg}{\mathrm{m}}^{-\mathrm{3}}$
Fresh snow density: 0	Empirical	${\mathit{\rho }}_{\mathrm{f}}=\max \left[a_{\mathrm{f}},+,b_{\mathrm{f}},\left(T_{\mathrm{a}},-,T_{\mathrm{m}}\right),+,c_{\mathrm{f}},U_{\mathrm{a}}^{\mathrm{1}/\mathrm{2}},,,{\mathit{\rho }}_{\min }\right]$
1	Empirical	${\mathit{\rho }}_{\mathrm{f}}={\mathit{\rho }}_{\min }+\max \left[d_{\mathrm{f}},{\left(T_{\mathrm{a}},-,T_{\mathrm{m}},+,e_{\mathrm{f}}\right)}^{\mathrm{3}/\mathrm{2}},,,\mathrm{0}\right]$
2	Constant	${\mathit{\rho }}_{\mathrm{f}}=100$ $\mathrm{kg}{\mathrm{m}}^{-\mathrm{3}}$
Thermal conductivity: 0	Empirical	${\mathit{\lambda }}_{\mathrm{s}}={\mathit{\lambda }}_{\mathrm{a}}+\left(a_{\mathit{\lambda }},{\mathit{\rho }}_{\mathrm{s}},+,b_{\mathit{\lambda }},{\mathit{\rho }}_{\mathrm{s}}^{\mathrm{2}}\right)\left({\mathit{\lambda }}_{\mathrm{i}},-,{\mathit{\lambda }}_{\mathrm{a}}\right)$
1	Empirical	${\mathit{\lambda }}_{\mathrm{s}}=c_{\mathit{\lambda }}{\left(\frac{{\mathit{\rho }}_{\mathrm{s}}}{{\mathit{\rho }}_{\mathrm{w}}}\right)}^{n_{\mathit{\lambda }}}$
2	Constant	${\mathit{\lambda }}_{\mathrm{s}}=0.265$ $\mathrm{W}{\mathrm{m}}^{-\mathrm{1}}{\mathrm{K}}^{-\mathrm{1}}$
Maximum liquid water: 0	Empirical	$\frac{{\mathit{\gamma }}_{\mathrm{w},\max }}{{\mathit{\rho }}_{\mathrm{s}}}=r_{\min }+\left(r_{\max },-,r_{\min }\right)\max \left(\mathrm{1},-,\frac{{\mathit{\rho }}_{\mathrm{s}}}{{\mathit{\rho }}_{\mathrm{r}}},,,\mathrm{0}\right)$
1	Constant	${\mathit{\gamma }}_{\mathrm{w},\max }={\mathit{\rho }}_{\mathrm{w}}\left(\mathrm{1},-,\frac{{\mathit{\gamma }}_{\mathrm{i}}}{{\mathit{\rho }}_{\mathrm{i}}}\right)S_{\mathrm{wi}}$
2	None	${\mathit{\gamma }}_{\mathrm{w},\max }=\mathrm{0}$ $\mathrm{kg}{\mathrm{m}}^{-\mathrm{3}}$

JIM variables are snow density ${\mathit{\rho }}_{\mathrm{s}}$, overlying snow mass $M_{\mathrm{s}}$, snow temperature $T_{\mathrm{s}}$, air temperature $T_{\mathrm{a}}$, wind speed $U_{\mathrm{a}}$, snow effective thermal conductivity ${\mathit{\lambda }}_{\mathrm{s}}$, partial density of liquid water ${\mathit{\gamma }}_{\mathrm{w}}$, and partial density of ice ${\mathit{\gamma }}_{\mathrm{i}}$. Other symbols represent constants, given in .

Models and methods

This study builds on the work of , who incorporated many published snow model parameterisations within a single model framework, the JIM, which is described in Sect. . As this earlier study did not incorporate snow microstructure changes, JIM was coupled with three microstructure evolution functions for this study, described in Sect. , and three distinct snow emission models, detailed in Sect. . Steps necessary to form the model ensemble, including assumptions about the representation of the soil, are given in Sect. . A description of the field site, driving, and evaluation data for the simulations in this paper are presented in Sect. .

Snow model parameterisation

developed the JIM, a system of 1701 snowpack evolution models to provide a systematic method and common framework to examine how the range of snowpack parameterisations used in land surface models impacts the simulation of snow parameters. Based on this work, a more computationally efficient version, a factorial snowpack model has been developed that allows for 32 model configurations. JIM is based on an Eulerian grid scheme (fixed layer structure), which requires mass redistribution between layers with precipitation events. An alternative approach is a Lagrangian grid scheme: a deforming layer structure that retains much of the same snow material throughout the season e.g.. For this paper, a subset of the original JIM members was selected as these were expected to influence the parameters important for microwave modelling. The subset includes variation in the representation of compaction, the density of newly deposited snow, thermal conductivity, and liquid-water flow (snow hydrology). Table  summarises the different approaches taken. Note that a variable fresh snow density (options 0 and 1) cannot be used when the snowpack has fixed density (compaction option 2), so there are only 63 model configurations in the model subset rather than 81. For all other snowpack parameterisations, option “1” from were used for albedo, surface exchange, and snow fraction representations to form the JIM subset.

Microstructure evolution

JIM subset outputs were used to drive three microstructure models of differing complexity. SNT growth of snow grain diameter $d$ is based on the rate of vapour transport through the snow (and therefore temperature gradient), which leads to the microstructure evolution function of dry snow in SNT as $${\displaystyle \frac{\partial d}{\partial t}}={\displaystyle \frac{g_{\mathrm{1}}}{d}}{D}_{\mathrm{eos}}\left({\displaystyle \frac{1000}{P_{\mathrm{a}}}}\right){\left({\displaystyle \frac{T_{\mathrm{s}}}{T_{\mathrm{m}}}}\right)}^{\mathrm{6}}{C}_{k{T}_{\mathrm{s}}}\left|{\displaystyle \frac{\partial T_{\mathrm{s}}}{\partial z}}\right|,$$ where $g_{\mathrm{1}}$ and $D_{\mathrm{eos}}$ are empirical constants, $P_{\mathrm{a}}$ is the atmospheric pressure, $C_{k{T}_{\mathrm{s}}}$ is the variation of saturation vapour pressure with snow temperature $T_{\mathrm{s}}$, $T_{\mathrm{m}}=273.15$ $\mathrm{K}$, and $\frac{\partial T_{\mathrm{s}}}{\partial z}$ is the temperature gradient. Grain growth under wet conditions is more rapid, with empirical constant $g_{\mathrm{2}}$ and is dependent on the liquid fractional volume, $f_{\mathrm{l}}$ by $$\begin{array}{}& {\displaystyle \frac{\partial d}{\partial t}}& {\displaystyle }={\displaystyle \frac{g_{\mathrm{2}}}{d}}\left(f_{\mathrm{l}},+,0.05\right)f_{\mathrm{l}}<0.09,& {\displaystyle \frac{\partial d}{\partial t}}& {\displaystyle }={\displaystyle \frac{g_{\mathrm{2}}}{d}}\left(0.14\right)f_{\mathrm{l}}\ge \mathrm{0.09.}\end{array}$$

SNI microstructure evolution is a computationally efficient approximation to a model based on physics that uses a look-up table for empirical parameters $\mathit{\tau }$ and $\mathit{\kappa }$, as described in . These parameters are dependent on the snow density, temperature, and temperature gradient. The equation of microstructure evolution in SNI is based on snow SSA: $${\displaystyle }\mathrm{SSA}\left(t\right)={\mathrm{SSA}}_{\mathrm{0}}{\left({\displaystyle \frac{\mathit{\tau }}{t+\mathit{\tau }}}\right)}^{\mathrm{1}/\mathit{\kappa }}.$$ SSA per unit mass of ice (${\mathrm{m}}^{\mathrm{2}}{\mathrm{kg}}^{-\mathrm{1}}$) can then be converted to grain diameter with $D=\mathrm{6}/\left({\mathit{\rho }}_{\mathrm{i}},\mathrm{SSA}\right)$ .

A third microstructure model, MOS, parameterises snow evolution as a function of grain radius $r$ and snow age: $${\displaystyle }r\left(t,+,\mathrm{\Delta },t\right)={\left[r,(,t,{)}^{\mathrm{2}},+,{\displaystyle \frac{G_{\mathrm{r}}}{\mathit{\pi }}},\mathrm{\Delta },t\right]}^{\mathrm{1}/\mathrm{2}}-\left[r,\left(,t,\right),-,r_{\mathrm{0}}\right]{\displaystyle \frac{S_{\mathrm{f}}\mathrm{\Delta }t}{d_{\mathrm{0}}}},$$ where $G_{\mathrm{r}}$ is an empirical temperature-dependent grain area growth rate, $S_{\mathrm{f}}$ is the snowfall rate in time interval $\mathrm{\Delta }t$, and $d_{\mathrm{0}}$ is a constant representing the mass of fresh snow required to reset the snow albedo to its maximum value.

Other microstructure parameterisations are available, namely the Crocus and SNOWPACK microstructure evolution functions. It is not currently possible to couple these with the JIM model due to the Eulerian grid structure of JIM. Mass transfer between layers allows numerical averaging of concepts such as grain diameter and SSA, but not shape-dependent concepts such as dendricity and sphericity. Therefore the Crocus and SNOWPACK functions have not been included in this study.

Microwave emission models

The microwave models chosen for this application span a range of physical complexity in their representation of the snow. The HUT model is a semi-empirical model based on strong forward scattering assumptions, the MEMLS model is of intermediate complexity and contains the improved Born approximation , and the DMRT-ML model is the most physically complex and is based on quasi-crystalline approximation with coherent potential (QCA-CP). Many other microwave emission models have been developed, such as Mie scattering approach of , , and , strong fluctuation theory , distorted Born approximation , the quasi-crystalline approximation , other QCA-CP models , or the numerical method of Maxwell's equations in 3-D . These references are not exhaustive but do give an illustration of the range of models available. Here, we restrict the comparison to widely available multi-layer models that span a range of complexity and whose computational efficiency is such that entire seasons can be simulated rapidly.

Of the models chosen, all are multiple layer and broadly require the same information; i.e. they use layered information on snow temperature, density, and layer thickness as input but differ in their representation of the microstructure. They are all based on radiative transfer theory, which is governed by the following general equation: $$\begin{array}{cc}{\displaystyle }& {\displaystyle }\mathit{\mu }{\displaystyle \frac{\partial {\mathit{T}}_{\mathrm{B}}\left({\mathit{\theta }}_{\mathrm{s}},,,{\mathit{\varphi }}_{\mathrm{s}},,,z\right)}{\partial z}}={\mathit{\kappa }}_{\mathrm{a}}T\left(z\right)+{\displaystyle \frac{\mathrm{1}}{\mathrm{4}\mathit{\pi }}}\underset{\mathrm{4}\mathit{\pi }}{\int }\mathrm{\Psi }\left({\mathit{\theta }}_{\mathrm{s}},,,{\mathit{\varphi }}_{\mathrm{s}},;,{\mathit{\theta }}_{\mathrm{i}},,,{\mathit{\varphi }}_{\mathrm{i}}\right)\\ & {\displaystyle }& {\displaystyle }\cdot {\mathit{T}}_{\mathrm{B}}\left({\mathit{\theta }}_{\mathrm{i}},,,{\mathit{\varphi }}_{\mathrm{i}},,,z\right)\mathrm{d}{\mathrm{\Omega }}_{\mathrm{i}}-{\mathit{\kappa }}_{\mathrm{e}}\cdot {\mathit{T}}_{\mathrm{B}}\left({\mathit{\theta }}_{\mathrm{s}},,,{\mathit{\varphi }}_{\mathrm{s}},,,z\right),\end{array}$$ where $\mathit{\theta }$ and $\mathit{\varphi }$ are the zenith and azimuth angles, $\mathit{\mu }=\cos \mathit{\theta }$, ${\mathit{T}}_{\mathrm{B}}$ is the brightness temperature vector, which we will assume here to consist of horizontally and vertically polarised brightness temperature only, ${\mathit{\kappa }}_{\mathrm{a}}$ is the absorption coefficient, and ${\mathit{\kappa }}_{\mathrm{e}}$ is the extinction coefficient, which is a sum of the absorption coefficient and the scattering coefficient ${\mathit{\kappa }}_{\mathrm{s}}$. The models differ in which two coefficients determine the third. In HUT, the derived coefficient is ${\mathit{\kappa }}_{\mathrm{s}}$, whereas ${\mathit{\kappa }}_{\mathrm{e}}$ is derived in MEMLS and ${\mathit{\kappa }}_{\mathrm{a}}$ in DMRT-ML. Other differences between models include the representation of the phase function (single-stream model with separate up- and downwelling components in HUT, six-stream in MEMLS and multiple streams in DMRT-ML), specification of the absorption coefficient and the numerical techniques applied to solve the radiative transfer equation . Differences between models are not restated here, but options chosen within each model leading to different model versions are stated in the following subsections.

DMRT-ML

DMRT-ML is based on a sticky hard spheres representation of the microstructure so that the scattering coefficient given by the QCA-CP is given as $${\displaystyle }{\mathit{\kappa }}_{\mathrm{s}}={\displaystyle \frac{\mathrm{2}}{\mathrm{9}}}{k}_{\mathrm{0}}^{\mathrm{4}}{a}^{\mathrm{3}}f{\left|{\displaystyle \frac{{\mathit{ϵ}}_{\mathrm{s}}-{\mathit{ϵ}}_{\mathrm{b}}}{\mathrm{1}+\frac{{\mathit{ϵ}}_{\mathrm{s}}-{\mathit{ϵ}}_{\mathrm{b}}}{\mathrm{3}{E}_{\mathrm{eff}}}(\mathrm{1}-f)}}\right|}^{\mathrm{2}}{\displaystyle \frac{(\mathrm{1}-f{)}^{\mathrm{4}}}{(\mathrm{1}+\mathrm{2}f-tf(\mathrm{1}-f){)}^{\mathrm{2}}}},$$ where $k_{\mathrm{0}}=\mathrm{2}\mathit{\pi }/\mathit{\lambda }$ is the wave number, a is the radius of the spheres, $f$ is the fractional volume of scatterers, ${\mathit{ϵ}}_{\mathrm{s}}$ is the permittivity of the scatterers, ${\mathit{ϵ}}_{\mathrm{b}}$ is the permittivity of the background, and $E_{\mathrm{eff}}$ is the effective permittivity of the medium. $t$ is related to the stickiness factor $\mathit{\tau }$ governing the potential of particles to coalesce. For non-sticky particles $t=\mathrm{0}$ but for sticky particles, it is given by the largest of the two solutions to the quadratic equation: $${\displaystyle \frac{f}{12}}{t}^{\mathrm{2}}-\left(\mathit{\tau },+,{\displaystyle \frac{f}{\mathrm{1}-f}}\right)t+{\displaystyle \frac{\mathrm{1}+f/\mathrm{2}}{(\mathrm{1}-f{)}^{\mathrm{2}}}}=0.$$ Whilst have shown that it may be possible to determine stickiness from micro-CT measurements of the snow, an appropriate value of stickiness is not known for the field observations used in this paper. and showed that non-sticky representation in DMRT-ML is inappropriate. For this model ensemble, two DMRT-ML configurations have been chosen to capture the range of brightness temperatures simulated: “DMRT less sticky” ($\mathit{\tau }=0.2$) and “DMRT very sticky” ($\mathit{\tau }=0.1$). These two values represent reasonable values used by others e.g..

MEMLS

Within MEMLS there are a suite of options for the calculation of the scattering coefficient. Two of the options within MEMLS were selected for this study to cover both empirical and theoretical approaches: “MEMLS empirical” and “MEMLS IBA”. The empirical version of MEMLS used gives the scattering coefficient as $${\displaystyle }{\mathit{\kappa }}_{\mathrm{s}}=(9.2p_{\mathrm{ec}}-1.23\mathit{\rho }+0.54{)}^{2.5}{\left(\mathit{\nu },/,50\right)}^{2.5},$$ where the correlation length $p_{\mathrm{ec}}$ is in $\mathrm{mm}$, density $\mathit{\rho }$ is in $\mathrm{g}{\mathrm{cm}}^{-\mathrm{3}}$, and frequency $\mathit{\nu }$ is in $\mathrm{GHz}$. This is suitable for correlation lengths $0.05<p_{\mathrm{ec}}<0.3$ $\mathrm{mm}$ and density $0.1<\mathit{\rho }<0.4$ $\mathrm{g}{\mathrm{cm}}^{-\mathrm{3}}$.

MEMLS IBA uses the improved Born approximation theory given in and , where the scattering coefficient is given by the integral of the phase function for polarisation angle $\mathit{\chi }$ $${\displaystyle }{\mathit{\kappa }}_{\mathrm{s}}={\displaystyle \frac{\mathrm{1}}{\mathrm{4}\mathit{\pi }}}\underset{\mathrm{4}\mathit{\pi }}{\int }f(\mathrm{1}-f){\left({\mathit{ϵ}}_{\mathrm{s}},-,{\mathit{ϵ}}_{\mathrm{b}}\right)}^{\mathrm{2}}{K}^{\mathrm{2}}I{k}_{\mathrm{0}}^{\mathrm{4}}{\sin }^{\mathrm{2}}\mathit{\chi }\mathrm{d}\mathrm{\Omega }.$$

One further assumption applied to distinguish this MEMLS IBA configuration is that oblate grains are used rather than small spherical scatterers or thin spherical shells. This assumption governs the representation of the mean square field ratio, $K^{\mathrm{2}}$, as detailed in . The microstructure length information is contained in $I$: $${\displaystyle }I={\displaystyle \frac{\mathrm{2}{p}_{\mathrm{ec}}}{(\mathrm{1}+\mathrm{4}{\mathit{ϵ}}_{\mathrm{eff}}{k}_{\mathrm{0}}^{\mathrm{2}}{\sin }^{\mathrm{2}}(\mathit{\theta }/\mathrm{2})p_{\mathrm{ec}}^{\mathrm{2}}{)}^{\mathrm{2}}}}.$$ It should be noted that the choice of oblate grains also affects the effective permittivity in $I$, represented by an empirical, density-dependent effective permittivity Eqs. 45–47 for this case.

HUT

HUT has three options for the extinction coefficient. These are nominally suited to different grain diameter ($d_{\mathrm{0}}$) ranges, with some overlap between them. All three versions (termed HUT H87, HUT R04, HUT K10) have been included in this version of the model ensemble. HUT H87 is based on the work of : $${\displaystyle }{\mathit{\kappa }}_{\mathrm{e}}=0.0018{\mathit{\nu }}^{2.8}{d}_{\mathrm{0}}^{1.9}.$$ This is nominally appropriate for frequency range $\mathit{\nu }=$ 18–60 $\mathrm{GHz}$ and $d_{\mathrm{0}}<1.6$ $\mathrm{mm}$.

The extinction coefficient in HUT H04, with a validity range of $1.3<d_{\mathrm{0}}<\mathrm{4}$ $\mathrm{mm}$ was derived by : $${\displaystyle }{\mathit{\kappa }}_{\mathrm{e}}=\mathrm{2}{\mathit{\nu }}^{0.8}{d}_{\mathrm{0}}^{1.2}.$$

gave the extinction coefficient for maritime snow, used here in the HUT K10 simulations as $${\displaystyle }{\mathit{\kappa }}_{\mathrm{e}}=0.08{\mathit{\nu }}^{1.75}{d}_{\mathrm{0}}^{1.8}.$$

Scaling of the grain diameter by the relationship recommended in has not been applied here as it was developed for snow microstructure observations rather than simulated snowpack microstructure.

Model framework

Interfacing of the various model inputs and outputs was enabled through the development of the ensemble framework, via a combination of shell script and Octave/MATLAB code. The DMRT-ML model was run from the shell script, which subsequently calls an Octave/MATLAB script to run HUT and MEMLS. HUT and MEMLS run alternately in this framework as the soil parameters (common between DMRT-ML and HUT) are used to calculate soil reflectivity in HUT, which is then used as the lower boundary condition in MEMLS. Internal parallelisation of the MATLAB code of HUT-MEMLS means that a season-long simulation of one HUT-MEMLS combination with one grain scaling factor takes 9 $\min$ over eight cores. For the DMRT-ML FORTRAN code, external bash shell parallelisation reduces execution time from 16 to ca. 2 $\mathrm{h}$ for one grain scale factor and two parameterisations of stickiness. Over 29 million individual brightness temperatures were simulated for this study.

Flowchart showing flow of information from JIM snow evolution model outputs to outputs from the various microwave emission models.

[Figure omitted. See PDF]

For the purposes of this study, the effective sphere size in JIM, DMRT-ML, and HUT is assumed to be identical i.e. $d_{\mathrm{HUT}}=\mathrm{2}\times r_{\mathrm{DMRT}}=d_{\mathrm{JIM}}$. This may not be a good assumption as the empirical extinction coefficient model used in HUT was based on observations of the maximum grain extent rather than effective diameter, which was almost impossible to measure at the time of the original work. The exponential correlation length in MEMLS (in $\mathrm{mm}$) is calculated from the theoretical relationship to the effective grain diameter from JIM (in $\mathrm{\mu }\mathrm{m}$) as and : $${\displaystyle }{p}_{\mathrm{ec}}={\displaystyle \frac{\mathrm{2}}{\mathrm{3}}}\left(\mathrm{1},-,{\displaystyle \frac{{\mathit{\rho }}_{\mathrm{s}}}{{\mathit{\rho }}_{\mathrm{i}}}}\right){\displaystyle \frac{d_{\mathrm{JIM}}}{1000}}.$$

Electromagnetic model inputs, as a function of JIM snowpack model outputs.

	JIM	DMRT-ML	MEMLS	HUT
	output	input	input	input
Temperature	$T_{\mathrm{JIM}}$ [$\mathrm{K}$]	$T_{\mathrm{JIM}}$ [$\mathrm{K}$]	$T_{\mathrm{JIM}}$ [$\mathrm{K}$]	$T_{\mathrm{JIM}}-273.15$ [${}^{\circ }\mathrm{C}$]
Density	${\mathit{\rho }}_{\mathrm{JIM}}$ [$\mathrm{kg}{\mathrm{m}}^{-\mathrm{3}}$]	${\mathit{\rho }}_{\mathrm{JIM}}$ [$\mathrm{kg}{\mathrm{m}}^{-\mathrm{3}}$]	${\mathit{\rho }}_{\mathrm{JIM}}$ [$\mathrm{kg}{\mathrm{m}}^{-\mathrm{3}}$]	$\frac{{\mathit{\rho }}_{\mathrm{JIM}}}{1000}$ [$\mathrm{g}{\mathrm{cm}}^{-\mathrm{3}}$]
Layer size	$\mathrm{\Delta }{z}_{\mathrm{JIM}}$ [$\mathrm{m}$]	$\mathrm{\Delta }{z}_{\mathrm{JIM}}$ [$\mathrm{m}$]	$\frac{\mathrm{\Delta }{z}_{\mathrm{JIM}}}{100}$ [$\mathrm{cm}$]	${\mathit{\rho }}_{\mathrm{JIM}}\mathrm{\Delta }{z}_{\mathrm{JIM}}$ [${\mathrm{mm}}_{\mathrm{swe}}$]
Microstructure	$d_{\mathrm{JIM}}$ [$\mathrm{\mu }\mathrm{m}$]	$d_{\mathrm{JIM}}$ [$\mathrm{\mu }\mathrm{m}$]	$\frac{\mathrm{2}}{\mathrm{3}}\left(\mathrm{1},-,\frac{{\mathit{\rho }}_{\mathrm{s}}}{{\mathit{\rho }}_{\mathrm{i}}}\right)\frac{d_{\mathrm{JIM}}}{1000}$ [$\mathrm{mm}$]	$\frac{d_{\mathrm{JIM}}}{1000}$ [$\mathrm{mm}$]
Layer number	$\mathrm{1}=$ base	$\mathrm{1}=$ base	$\mathrm{1}=$ base	$\mathrm{1}=$ top
Soil permittivity	–	${\mathit{ϵ}}_{\mathrm{obs}}$	$r_{\mathrm{0},\mathrm{HUT}}$	${\mathit{ϵ}}_{\mathrm{obs}}$

Note that for the purposes of the ensemble, DMRT-ML was adapted to allow the input of diameter rather than radius. HUT was adapted to ensure Fresnel reflectivity for a smooth soil surface and to output the soil reflectivity $r_{\mathrm{0},\mathrm{HUT}}$ at both polarisations for use in MEMLS simulations.

Air temperature, solar radiation, and precipitation data measured at the Sodankylä site, used as inputs for the JIM simulations.

[Figure omitted. See PDF]

Figure  illustrates the flow of information in the model ensemble. Meteorological data are used to drive the 189 configurations of JIM (3 microstructural models for each of the 63 snowpack parameterisations). The outputs from JIM are then reformatted for each of the three electromagnetic models. Table  gives a summary of the main differences in inputs between models. The electromagnetic model inputs are then used to drive the two DMRT-ML model versions ($\mathit{\tau }=0.1$, $\mathit{\tau }=0.2$), the two MEMLS model versions (empirical, IBA with oblate grains), and the three HUT versions (three different extinction coefficient models). Meteorological and field data used to drive and evaluate the ensemble are described in the following section.

Data

Model runs for this study were performed for the Intensive Observation Area (IOA) of the Finnish Meteorological Institute Arctic Research Centre (FMI-ARC). The site provides a wealth of forcing and evaluation data, including automated soil, snow and meteorological observations, ground-based microwave radiometry, and a programme of manual snow profile observations. Air temperature, solar radiation, and precipitation observations from this site for the two seasons of simulations are shown in Fig. . November rain events occurred in both years, as well as in early December in 2011–2012. Layers with melt–freeze polycrystals and other melt forms were detected in snow observations during both seasons. Metadata and details on the meteorological instruments are given in . Dual polarisation microwave radiometers, including at frequencies of 18.7 and 36.5 $\mathrm{GHz}$, are situated on a 4 m tower pointing inwards on the edge of a large clearing surrounded by a mainly pine forest. Further details about the IOA site are given in . Details on the manual snow profile observation programme are given by .

Simulations were carried out for the winters of 2011–2012 and 2012–2013 as there were 49 approximately bi-weekly snow pit observations over these 2 years available for snowpack model evaluation. Snow samples from 31 of these pits were extracted and used to measure profiles of the SSA with the IceCube instrument . A bulk grain diameter was calculated for the analysis from the snow water equivalent (SWE)-weighted mean SSA, excluding layers without observations. For these two seasons, the real component of the soil permittivity measurements were available at 100 $\mathrm{MHz}$, measured at three locations in the observation area with Delta-T devices ML2x sensors, installed horizontally at a depth of approximately 2 $\mathrm{cm}$ beneath the organic surface layer. Mean measurements from the stable winter period (1 December–31 March) were chosen as representative for the entire season, which resulted in values of soil permittivity of 4.4 in 2011–2012 and 4.6 in 2012–2014. For the JIM simulations in this paper, a scaling factor of 1.11 was applied to the 2011–2012 precipitation data, and a scaling factor of 1.06 was applied to the 2012–2013 data to match the measured snow accumulation on the ground better. These factors differ slightly from the values used in the 7-year consolidated data set of .

Snow depth and water equivalent simulated by the Jules Investigation Model subset used in this study. Grey lines indicate individual JIM subset member simulations. Note that erroneous positive SWE observation points have been removed at the end of the season when snow depth is zero, as this is a sensor artifact related to soil moisture changes.

[Figure omitted. See PDF]

Simulation methodology

Choices in the snowpack evolution parameterisations made here lead to 189 unique JIM snowpack models. Modelled snowpack profiles of layer thickness, temperature, density, and grain diameter were output daily at noon for this study. These were then applied to the seven microwave emission model combinations, resulting in 1323 sets of brightness temperature simulations per day.

In order to illustrate and analyse the effects of assumptions regarding snowpack evolution and microwave scattering on simulated brightness temperatures over the course of the winter season, the remainder of the paper will do the following:

• Present the range of brightness temperatures expected for any generic combination of snowpack and emission model.

• Apply a range of scaling factors ($0.1\le \mathrm{\Phi }\le 5.0$) to simulated JIM snowpack diameters ($d_{\mathrm{optimal}}=\mathrm{\Phi }{d}_{\mathrm{JIM}}$) and calculate the degree of misfit between simulated and observed brightness temperatures using the following cost function (CF):$${\displaystyle }\mathrm{CF}=\sum ^{\mathrm{ndays}}\sum ^{\mathit{\nu }}\sum ^{\mathrm{pol}}{\left({\displaystyle \frac{T_{{\mathrm{B}}_{\mathrm{sim}}}-T_{{\mathrm{B}}_{\mathrm{obs}}}}{\mathrm{2}}}\right)}^{\mathrm{2}}.$$The cost function term is summed over the two polarisations (H- and V-pol) for the two frequencies (18.7 and 36.5 $\mathrm{GHz}$) over the number of days (ndays) when observations and simulations are both available. Due to the observation schedule at the Sodankylä site, the noon “observations” for comparison with the simulations were determined as the mean of the 10 am and 2 pm observations. If observations were missing from either or both of these times, the brightness temperature for that day was excluded from the CF calculation. Optimal $\mathrm{\Phi }$ were found from the minimisation of the CF.

• Isolate the effect of snowpack parameterisations on simulated brightness temperature by presenting simulation results grouped by parameterisations of densification, liquid-water flow, initial snow density, and thermal conductivity. This will determine which factors govern the spread in brightness temperature and are therefore important for the design of snow retrieval assimilation systems.

Results

Snow depth and SWE simulated by the JIM is shown in Fig. . There is a small difference between the automatic measurements and the manual field observations attributable to the spatial variability of the snow and difference in measurement location. Ultrasonic snow depth measurements were on average 12 mm deeper than the snow pit observations in 2011–2012 but were 29 mm shallower than snow pit measurements in 2012–2013. SWE measured automatically by the gamma ray sensor had a mean value of 3.6 mm SWE greater than the pit observations in 2011–2012 but 5.9 $\mathrm{mm}$ less in 2012–2013.

Although the precipitation inputs were scaled due to known sensor undercatch problems, in 2011–2012 the SWE was underestimated until the end of January, then overestimated until the melt period. Compared with snow pit observations, the SWE bias prior to 1 February was $-$14.2 $\mathrm{mm}$. Between 1 February and 31 March, the SWE bias was 19.1 $\mathrm{mm}$. From 1 April until the end of the season, the bias was $-$24.5 $\mathrm{mm}$ water equivalent. In 2012–2013, simulated SWE was overestimated for most of the season, with a mean bias of 13.7 $\mathrm{mm}$ compared with the snow pit observation. Simulated SWE is relatively insensitive to the snow parameterisation in the accumulation period, but three distinct model groups emerge in the melt period, which are due to the three different representations of the liquid-water flow. The snowpack model parameterisations have a greater impact on the snow depth, which is to be expected as this is directly affected by the representation of densification and initial snow density.

Figure  demonstrates the impact of the snow model parameterisations on snow grain diameter growth, as simulated with the MOS, SNI, and SNT microstructure evolution functions. Each microstructure model results in a spread of bulk grain diameter due to the 63 snowpack parameterisations, but in general the difference between microstructure models is greater than the difference due to snowpack parameterisations. The simulation range is greatest at the start and at the end of the season, when the snowpacks can be subject to the largest temperature gradients or liquid-water dependent growth. In both years of simulation, the mid-season bulk grain diameter is smallest with MOS and largest with SNT. MOS and SNI are similar in magnitude, but SNT bulk grain diameters were approximately twice as large on average, with a mean ratio over the season of 1.9–2.2, as shown in Table . SNT bulk grain diameter was up to 3.2 times larger than MOS bulk grain diameter. Visual estimation of the snow grain diameter gave values that were always larger than all of the simulations. Measured SSA-derived bulk grain diameters generally lay in between the SNT simulations and the simulations with SNI and MOS. The mean absolute error and mean relative difference for these simulations are presented in Table . SNT had the lowest bias (0.12 $\mathrm{mm}$) in 2011–2012, whereas SNI had the lowest bias ($-$0.14 $\mathrm{mm}$) in 2012–2013. Bulk grain diameter simulated by the microstructure models led to a mean difference of between $-$53 and $+$45 % relative to the observations.

Bulk grain diameter evolution for the MOS, SNT, and SNI microstructure evolution models and the spread in model results. Observations of bulk diameter were derived from macro-photography (Visual) and from SSA measurements from the IceCube instrument.

[Figure omitted. See PDF]

Comparison of grain diameters simulated by different microstructure models. The mean and max ratio between pairs of models is given in columns. Where the 2012–2013 values differ, these are given in parentheses.

	Mean	Max
SNI $/$ MOS	1.2 (1.1)	1.4 (1.3)
SNT $/$ MOS	2.2	3.1 (3.2)
SNT $/$ SNI	1.9 (2.0)	2.5

Simulation of mean and range of brightness temperature from the three emission models driven by all snowpack and microstructure model combinations is shown in Fig. . Note that excessively low brightness temperatures on 1 November 2011 were excluded from this figure as the snowpack for some JIM members was extremely thin with an unphysically high snow density. In general, HUT with three representations of extinction coefficient showed the smallest range of brightness temperature, whereas DMRT-ML (covering both very sticky and less sticky assumptions) had a much greater range, which was nearly as large as MEMLS (empirical representation and improved Born approximation with oblate grains). This is demonstrated by the ratio between the seasonal mean ranges of brightness temperature presented in Table , where the ranges compared with HUT had a ratio of greater than 1. MEMLS had a larger range than DMRT-ML, although at 37 $\mathrm{GHz}$ the difference was small. As illustrated in Fig. , at 19 $\mathrm{GHz}$, the mean of DMRT-ML simulations were highest and the mean of MEMLS simulations were generally lowest (with the exception of 19H in 2011–2012). At 37 $\mathrm{GHz}$, horizontal and vertical polarisation, HUT gives the highest mean brightness temperature in both years, although the mean DMRT-ML brightness temperature is within 3 K of HUT at horizonal polarisation (both years). MEMLS mean brightness temperatures are the lowest at 37 $\mathrm{GHz}$ at both horizontal and vertical polarisation in both years. All ranges exhibit a distinctive “wedge” shape, where the ranges generally increase throughout the season until the collapse of the range in the melt period.

Mean absolute error ($\mathrm{mm}$) between bulk grain diameter simulated with the microstructure models compared with observations derived from SSA measurements with IceCube. Smallest bias for each year is shown in bold. Percentages are given in parentheses.

	2011–2012	2012–2013
MOS	$-$0.24 ($-$53 %)	$-$0.16 ($-$34 %)
SNI	$-$0.18 ($-$40 %)	$\mathbf{-}$0.14 ($-$31 %)
SNT	0.12 (32 %)	0.16 (45 %)

Ratio of mean brightness temperature ranges simulated by two microwave emission models. The mean and max ratio between pairs of models is given in columns.

	19V	19H	37V	37H
2011–2012
DMRTML/HUT	1.6	1.2	3.8	3.1
MEMLS/HUT	3.7	1.9	4.5	3.2
MEMLS/DMRTML	2.4	1.5	1.2	1.0
2012–2013
DMRTML/HUT	2.0	1.6	3.3	3.0
MEMLS/HUT	3.9	2.4	3.9	3.2
MEMLS/DMRTML	2.0	1.5	1.2	1.1

Range and mean of brightness temperature over the two winter seasons as simulated with the DMRT-ML, MEMLS, and HUT models, driven by 63 JIM outputs and 3 microstructure evolution models. Black lines indicate the observed brightness temperatures. Vertical dashed lines enclose the period of analysis (1 November–31 March).

[Figure omitted. See PDF]

Compared with the brightness temperature observations, no model gives a consistently better performance across both frequencies and both polarisations. This is illustrated in Table , where mean bias and root mean square error (RMSE) for each season has been presented for each frequency and polarisation combination. The lowest bias was less than 7 $\mathrm{K}$ in magnitude, whereas the lowest RMSE for each frequency/polarisation was less than 13 $\mathrm{K}$. For both years, DMRT-ML gave the lowest bias at 19H and MEMLS gave the lowest bias at 37 $\mathrm{GHz}$ (V and H). At 19V DMRT-ML had the lowest bias in 2011–2012 whereas HUT had the lowest bias in 2012–2013. Figure  shows that the observed brightness temperature is generally within the range simulated by each of the three microwave emission models, with the exception of 19H in 2011–2012 (MEMLS and HUT) and 37H/V in 2012–2013 (HUT). End of season brightness temperature observations are not replicated in the simulations as the liquid-water content of the snowpack model is currently decoupled from the electromagnetic snow model, so the simulations only represent dry snow brightness temperature.

Mean bias and RMSE in brightness temperature (K) simulated by DMRT-ML (very sticky and less sticky), MEMLS (empirical and IBA oblate), and HUT (H87, R04, and K10) forced by 189 JIM–microstructure model combinations. Only days in the period from 1 November to 31 March, where all four frequency/polarisation measurements were available, were included in the analysis. Bold values indicate the lowest bias/RMSE for each frequency/polarisation.

		19V	19H	37V	37H
2011–2012
	DMRTML	0.7	$\mathbf{-}$5.4	10.4	8.7
Bias	MEMLS	$-$7.8	$-$16.3	$\mathbf{-}$6.9	$\mathbf{-}$6.2
	HUT	$-$1.3	$-$18.8	20.6	7.6
	DMRTML	5.5	11.4	12.2	12.2
RMSE	MEMLS	11.3	20.4	13.0	13.1
	HUT	6.1	22.4	21.1	11.7
2012–2013
	DMRTML	5.6	6.9	25.2	24.5
Bias	MEMLS	$-$9.3	$-$11.4	$\mathbf{-}$0.9	1.5
	HUT	2.9	$-$8.0	39.2	26.4
	DMRTML	6.4	9.8	26.2	27.0
RMSE	MEMLS	11.4	13.2	7.0	7.0
	HUT	4.2	9.8	40.2	28.9

Optimal microwave microstructure scale factors dependent on snow microstructure evolution function, based on minimisations of cost function between 1 November and 31 March in each year.

	DMRTML			MEMLS			HUT
	less	very		IBA	EMP		H87	R04	K10
2011–2012
SNT	1.1	0.6		0.5	0.3		0.9	0.5	0.7
MOS	3.3	1.6		1.7	1.0		2.6	1.4	2.2
SNI	2.5	1.3		1.2	0.8		1.9	1.1	1.7
2012–2013
SNT	1.3	0.7		0.7	0.5		1.2	1.1	1.1
MOS	3.1	1.7		1.6	1.1		3.2	2.7	2.9
SNI	2.8	1.5		1.5	1.1		2.9	2.3	2.6

A value of 1.0 indicates that the snow grain diameter simulated by a particular form of the snow model may be used directly in the microwave model to give the best agreement with measured brightness temperature.

Table  indicates scaling factors that would need to be applied to the grain diameter in order to allow a particular microstructure evolution function to minimise the CF given in Eq. (), i.e. the best agreement with observed brightness temperature for all four frequency and polarisation combinations. A scale factor of 1 suggests a perfect fit between snowpack microstructure and microwave microstructure. A scale factor of less than 1 indicates a snowpack grain diameter overestimate, whereas a scale factor of greater than 1 is an underestimate. For SNT microstructure, a scale factor of less than 1 was required in 2011–2012 for all emission models with the exception of the less sticky ($\mathit{\tau }=0.2$) application of DMRT-ML. This indicates that the SNT microstructure resulted in grain diameters larger than that required by the emission models for that year. In 2012–2013 SNT microstructure required slight scaling to increase the grain diameter for HUT and for less sticky DMRT-ML, but downscaling for very sticky hard spheres in DMRT-ML and for MEMLS. With the exception of the application to empirical MEMLS in 2011–2012, the SNI and MOS grain diameters were too small and required scaling upwards. A CF minimum was achieved for empirical MEMLS driven by MOS microstructure with no scaling whatsoever in 2011–2012. The pattern is consistent between years, with the greatest interannual difference in scale factor for HUT.

Mean bias and RMSE in brightness temperature (K) simulated by DMRT-ML (sticky and non-sticky), MEMLS (empirical and IBA oblate), and HUT (H87, R04, and K10) forced by 189 JIM–microstructure model combinations, with optimal microstructure scale factors from Table  applied. Only days in the period from 1 November to 31 March, where all four frequency/polarisation measurements were available, were included in the analysis. Bold values indicate the lowest bias/RMSE for each frequency/polarisation.

		19V	19H	37V	37H
2011–2012
	DMRTML	1.0	$\mathbf{-}$5.1	6.6	5.6
Bias	MEMLS	$\mathbf{-}$0.3	$-$11.3	13.1	10.5
	HUT	$-$1.3	$-$18.8	19.1	6.2
	DMRTML	5.6	11.2	9.3	9.5
RMSE	MEMLS	5.7	15.5	14.0	12.7
	HUT	6.1	22.4	19.5	10.7
2012–2013
	DMRTML	3.7	5.6	$\mathbf{-}$0.3	1.4
Bias	MEMLS	$\mathbf{-}$1.2	$\mathbf{-}$5.1	6.0	9.1
	HUT	$-$4.4	$-$14.3	18.0	7.5
	DMRTML	4.7	8.8	5.2	7.5
RMSE	MEMLS	3.9	7.6	9.3	11.7
	HUT	6.2	15.4	19.8	11.2

Once the microstructure differences have been isolated through application of the optimal scale factor, as shown in Table , DMRT-ML bias and RMSE improved, with the exception of the small increase in 19V bias in 2011–2012. For MEMLS, improvements in bias and RMSE at the lower frequency were at the expense of the higher frequency in both years. The opposite occurred for HUT in 2012–2013, whereas in 2011–2012 the bias and RMSE decreased at all frequencies and polarisations apart from a marginal ($<$ 0.04 $\mathrm{K}$) increase in RMSE at 19V.

Differences in brightness temperature also exist in the simulations due to the snowpack parameterisation (i.e. 63 JIM combinations). Empirical MEMLS with MOS microstructure in the 2011–2012 season was chosen as a test case to illustrate the effects of snowpack parameterisation on the brightness temperature because of the equivalence of snowpack and emission model microstructure (no scaling required). This subset of 63 simulations for 37H brightness temperature in 2011–2012 is shown in Fig. . There is a seasonal dependence in the range, with model divergence from mid-January onwards. 1 February and 1 May were chosen for cluster analysis to determine which parameterisations caused the split in simulations, as shown in Fig. .

Variability in brightness temperature simulated with empirical MEMLS, driven by the MOSES microstructure model and 63 JIM snowpack outputs (no scaling of microstructure was required).

[Figure omitted. See PDF]

Clear groupings of simulations in Fig. , upper left, indicate that the snowpack densification parameterisation has a distinguishable effect on the simulation of brightness temperature. A physical representation of densification (parameterization $=$ 0) gave the lowest brightness temperatures on the 1 February, but the highest by 1 May. In contrast, where no compaction is simulated, i.e. snow density is constant throughout the season (parameterization $=$ 2), the opposite is true. An empirical representation of densification (parameterization $=$ 1) results in brightness temperatures generally between those of the physical, and of no densification. Thermal conductivity has no effect on the simulation of brightness temperature, whereas subtle differences are attributable to the fresh snow density value and to the representation of snow hydrology. There is no discernible difference between fresh snow density parameterisation schemes 0 and 1, whereas 2 gives a different set of brightness temperatures. Snow hydrology has very little effect in the early season but can lead to differences in the melt period. Overall, the snowpack parameterisations with MOSES microstructure and empirical MEMLS lead to a mean difference in the 36.5 $\mathrm{GHz}$ brightness temperature of 11 $\mathrm{K}$ at H-pol and 18 $\mathrm{K}$ at V-pol. The maximum difference in 36.5 $\mathrm{GHz}$ brightness temperature was 33 $\mathrm{K}$ at H-pol and 54 $\mathrm{K}$ at V-pol for the 2011–2012 season. The maximum difference between H and V polarisation for all unscaled microstructure–electromagnetic model combinations is demonstrated in Table . Large differences in the maximum brightness temperature difference as a result of the 63 snowpack configurations occurred for the SNT microstructure. Except for DMRT-ML less sticky and HUT with MOS or SNI microstructure, the V-pol difference is greater than the H-pol difference.

Cluster analysis of brightness temperature simulated by MOSES microstructure and empirical MEMLS for 1 February 2012 (black dots) and 1 May 2012 (red dots) according to model parameterization choices. Brightness temperature simulations are split according to the different representations for each process representation. Values 0, 1, 2 relate to parameterisations given in and as described in Sect. . Where distinct clusters occur that differ between parameterisations, this indicates sensitivity to the parameterization.

[Figure omitted. See PDF]

Maximum difference in brightness temperature in K (H-pol/V-pol) due to 63 snowpack parameterisations for each unscaled snow microstructure evolution function (2011–2012 season).

	DMRTML			MEMLS			HUT
	less	very		IBA	EMP		H87	R04	K10
SNT	50/63	148/169		88/113	126/153		27/34	26/33	28/36
MOS	22/15	41/45		24/25	33/54		20/9	21/14	20/11
SNI	22/15	59/67		24/43	55/77		20/13	21/18	20/15

Discussion

The biggest difference to obtaining accurate simulations would be made by improving the microstructure evolution models within snowpack models because the optimal scale factors are generally larger between microstructure models than between emission models. SNTHERM grains tend to be too large for the emission models and generally require scaling down to smaller values. SNICAR grains are in the mid-range and require a small amount of scaling, generally upwards to larger grains. MOSES grains are the smallest and generally require higher scale factors than SNICAR. These patterns are consistent, regardless of the electromagnetic radiative transfer model used. Differences between microstructure evolution models are so large because they were developed in models with different purposes. MOSES is a large-scale land surface model, requiring snow grain size for albedo calculations . SNICAR is a snow albedo model . SNTHERM, in contrast, was developed to predict surface temperature and uses grain diameter in the simulation of liquid-water flow as well as albedo . SNICAR and MOSES grain sizes are closer to the SSA-derived grain diameter as a result. SNTHERM simulates a grain size that is closer in concept to the visual estimates of grain diameter than the other two models. The large spread when coupling snowpack evolution and microwave models, due to the differences in the modelling of snow microstructure, is consistent with the wide range of studies that have investigated how to link snowpack observations of microstructure to the microstructure parameter required in electromagnetic models e.g.

Nevertheless there are differences between microwave emission models for a particular microstructure evolution model and even differences within the same family of emission models. “Improvement” in the microstructure for a particular model combination may lead to less accurate simulations at some frequencies and polarisations, which highlights that there is more to understand. In part, this may be due to the methodology of this study as the CF is calculated per microstructure–electromagnetic model configuration, yet the bias and RMSE are presented for each electromagnetic model family. An individual contribution can influence the group in a non-intuitive way.

Here, the lowest bias and RMSE for unscaled microstructure simulations were $-$6.9 to $+$6.9 $\mathrm{K}$ and 4.2 to 12.2 $\mathrm{K}$, respectively, but depended on microwave model, frequency, and polarisation. In an attempt to put these results into context, there are a number of studies that have quantified brightness temperature simulation errors for these models. These fall into different categories, depending on sensor characteristics, the source of the evaluation data (ground-based, airborne, satellite) and presence of ice lenses , the treatment of the snow microstructure , snow type, observation angle, and the specific electromagnetic model , and the underlying substrate . Examples of unscaled field observations of microstructure compared with ground-based observations include the HUT simulations of , who found an RMSE of 10–34 $\mathrm{K}$, and , who found a bias of 34–68 $\mathrm{K}$ that was reduced to $<$ 0.6 $\mathrm{K}$ upon application of grain scale factors of 2.6–5.3. Scaling, or best-fit, relationships were used by (mean absolute error 3.1 $\mathrm{K}$ at V-pol and 9.3 $\mathrm{K}$ at H-pol), (RMSE 8–20 $\mathrm{K}$), (RMSE 1.5 $\mathrm{K}$), (RMSE 1–11 $\mathrm{K}$), and (RMSE 12–16 $\mathrm{K}$). However, in some cases the frequency-dependent results have been combined and in others kept separate.

For DMRT-ML, consideration of the stickiness is imperative. Two constant values were considered here: extremely cohesive or less sticky particles. have made progress in understanding stickiness from micro-CT data. There are theoretical limits, based on snow density Eqs. 35–36, but in general stickiness is independent of diameter and of density and a constant value should not be used, as was done here. Further research is needed in this regard.

For the HUT radiative transfer model, the optimum combinations of snowpack and microwave model are dependent on both models and, therefore, the end application. The SNT microstructure is most closely matched to the microstructure of the extinction model. Both were developed with a similar concept of microstructure. With MOS or SNI microstructure, would be most appropriate. is more broadly applicable as the scale factor always lies between R04 and H87, regardless of the microstructure model. Therefore K10 may be better choice if a range of microstructure models is considered in a data assimilation retrieval scheme but with only one observation operator.

In the case of MEMLS, there are some differences between the empirical model and IBA, but the microstructure model really matters. IBA is a more appropriate model for the larger SNT grains and endorses the recommendation of for IBA in the simulation of larger grains. The microstructural concept of MOS matches the microstructure of empirical MEMLS very well, with no scaling required in 2011–2012, although SNI is equally appropriate in 2012–2013.

There is little variation between years for the DMRT-ML (sticky) and MEMLS models, and a consistent pattern for HUT. Other studies have investigated the microstructural link between snowpack and microwave models. found that the scale factor between exponential correlation length in MEMLS and grain diameter in SNTHERM for the Weissfluhjoch site in Davos, Switzerland, was 0.16. Applying Eq. () for snow of density 250 $\mathrm{kg}{\mathrm{m}}^{-\mathrm{3}}$, the scale factor to relate the grain diameter of SNTHERM to the exponential length in MEMLS for the Sodankylä site would be 0.24 for IBA and 0.15 for empirical MEMLS for the 2011–2012 data set. This is entirely consistent with the study, in spite of the different locations and snowpack conditions.

also reported a relationship for Crocus simulations, as did . At this stage it is not possible to make comparisons of this work with those studies because the Crocus evolution model has not been included in this study due to the difficulty of applying these models to the Eulerian frame snowpack model scheme used here. These two studies are, however, consistent with each other. found a snow-type-dependent scale factor of 0.3–0.4 between MEMLS correlation length and Crocus grain diameter, whereas the range in was 0.4–0.25 for snow density between 100 and 400 $\mathrm{kg}{\mathrm{m}}^{-\mathrm{3}}$. The scaling factor between the SNOWPACK-derived correlation length and the correlation length of MEMLS was found to be 0.1 but, again, a comparison with this work is not possible as the SNOWPACK grain evolution model has similar requirements to the Crocus microstructure model as they have a common origin.

When isolating the spread in brightness temperature due to snowpack parameterisations, this spread is largely due to the snowpack model representation of the densification process, with a variable impact throughout the season. After the microstructure model, snow compaction must be considered carefully in the design of a coupled snowpack–microwave model. Liquid-water flow representation in the snowpack model may become important in the melt period, particularly for a snowpack with mid-winter melt periods or if the snowpack model is used to provide information on SWE during melt when microwave observations cannot. If fresh snow is assumed to have a constant density in a retrieval or assimilation system then that value will have an impact but is less important than compaction. Thermal conductivity has no discernable impact on the brightness temperature simulations so the choice of its representation is largely irrelevant for snow mass retrieval and assimilation systems.

Although empirical MEMLS driven by MOSES was chosen as an example to demonstrate the impact of parameterisations, this was purely because of the apparent consistency between the MOSES grain diameter converted to exponential correlation length and MEMLS simulations for 2011–2012 at this site. This is not a general endorsement of empirical models, as those based on physics are expected to be more universally applicable, but the specific application of these models will dictate the balance of accuracy versus simplicity. Extending the analysis beyond this example, snow parameterisations affect other unscaled model combinations to varying degrees. Microstructure scaling factors in Table  can be used as a proxy for the degree of scattering in the unscaled simulations. A higher scale factor acts to increase the simulated scattering, so for a scale factor $<\mathrm{1}$ too much scattering occurs in the unscaled simulations. Snowpack parameterisations have a greater impact for a higher degree of scattering, larger at V-pol than H-pol. This is because scattering is already greater at H-pol so the spread in H-pol simulations as a result of snowpack parameterisations is suppressed by the existing level of scattering. The converse applies for high scaling factors (e.g. MOS with less sticky DMRT-ML).

Although the differences in scale factors between microstructure models are larger than the differences in scale factors between microwave models, this does not negate the need for developments in the microwave models. This is highlighted by the treatment of field observations of SSA to derive optical diameter as even these require some form of scaling e.g.. Use of scale factors can improve brightness temperature accuracy at some frequency and polarisations but may decrease the accuracy at others. The necessity of scale factors indicate the need for a deeper understanding in the role of microstructure in the microwave models. Much of this is discussed from a theoretical perspective by . With a sticky hard sphere model of the microstructure, even if the stickiness is known, showed that a scale factor to relate the measured optical diameter to microwave diameter depends on the type of metamorphism the snow has been subjected to. Indeed, here, constant scale factors have been applied with no attempt to assess how these may change over the season. They do not account for the anisotropic nature of the snow, which adds to the complexity both in the modelling of the snowpack and in microwave scattering . Some of the fundamental questions on how to relate snowpack and microwave microstructure may be addressed with a better microstructure descriptor of the snowpack rather than a single length scale and would benefit from easy interchangeability between different microwave models and different snowpack evolution models. Ultimately a consistent microstructural treatment will be needed in both snowpack evolution and microwave models.

Conclusions

Future snow mass and depth retrievals systems may rely on snowpack models to provide snow microstructural parameters. To improve accuracy in seasonal simulations of brightness temperature, the largest gains will be achieved by improving the microstructural representation within snowpack models, followed by improvements in the emission models to use accurate microstructural information and reduce bias and RMSE at all frequencies and polarisations simultaneously. For the design of retrieval systems with current capabilities, particular model combinations may be more suitable than others, and careful consideration must be given to snow compaction processes. Snow process representation becomes increasingly important as the snowpack scatters more. The future lies in a better and consistent treatment of snow microstructure in both snowpack and emission model developments.

Code availability

Code to analyse the data is available on GitHub: https://github.com/mjsandells/TC_Sandells_2017.

Data availability

Data are available at the repository 10.6084/m9.figshare.4552822 . These include JIM outputs, the non-scaled grain diameter and optimal grain diameter brightness temperature simulations, and observations of brightness temperature and grain size. The full brightness temperature data set (all scaling factors) is too large to place in a repository but can be made available upon request.

The authors declare that they have no conflict of interest.

Acknowledgements

This work was funded in part by the NERC National Centre for Earth Observation, and supported by the European Space Agency projects “Technical assistance for the deployment of an X- to Ku-band scatterometer during the NoSREx experiment” (ESA ESTEC contract 22671/09/NL/JA/ef) and “Microstructural origin of electromagnetic signatures in microwave remote sensing of snow” (ESA ESTEC contract 4000112698/14/NL/LvH). We thank the staff of FMI Arctic Research Centre in Sodankylä for performing the in situ measurements. We thank G. Picard and an anonymous reviewer for their comments to improve this paper. Edited by: M. Tedesco Reviewed by: G. Picard and one anonymous referee