Effects of decimetre-scale surface roughness on L-band brightness temperature of sea ice

Sea ice thickness is an essential climate variable. Current L-Band sea ice thickness retrieval methods do not account for sea ice surface roughness that is hypothesised to be not relevant to the process. This study attempts to validate this hypothesis that has not been tested yet. To test this hypothesis, we created a physical model of sea ice roughness based on geometrical optics and merged it into the L-band emissivity model of sea ice that is similar to the one used in the operational sea ice thickness retrieval algorithm. The facet description of sea ice surface used in geometrical optics is derived from 2-D surface elevation measurements. Subsequently the new model was tested with T_B measurements performed during the SMOSice 2014 field campaign. Our simulation results corroborate the hypothesis that sea ice surface roughness has a marginal impact on near-nadir T_B (used in the current operational retrieval). We demonstrate that the probability distribution function of surface slopes can be approximated with a parametric function whose single parameter can be used to characterise the degree of roughness. Facet azimuth orientation is isotropic at scales greater than 4.3 km. The simulation results indicate that surface roughness is a minor factor in modelling the sea ice brightness temperature. The change in T_B is most pronounced at incidence angles greater than 40^∘ and can reach up to 8 K for vertical polarisation at 60^∘. Therefore current and future L-band missions (SMOS, SMAP, CIMR, SMOS-HR) measuring at such angles can be affected. Comparison of the brightness temperature simulations with the SMOSice 2014 radiometer data does not yield definite results.

1 Introduction

The L-band brightness temperature (T_B) is sensitive to sea ice thickness. This feature is used for sea ice thickness retrieval from L-band T_B (over thin ice, < 1.5 m) (Tian-Kunze et al., 2014; Huntemann et al., 2014; Kaleschke et al., 2016). Several factors influence T_B measurements over ice-covered regions, among them ice concentration, ice temperature, snow cover, sea ice surface roughness, and the shape of the interfaces between the snow and ice layers (Maaß et al., 2013; Ulaby et al., 2014, p. 422).

Here, we investigate the effects of surface roughness on the L-band T_B, specifically the large-scale roughness. So far this factor is not included in the modelling of sea ice emissions in operational sea ice thickness retrieval. Electromagnetic scattering theory assumes that the roughness of a random surface is characterised by statistical parameters including standard deviation of surface height (σ_z) and correlation function (R(ξ)) measured in units of wavelength (Ulaby et al., 2014, p. 422). These roughness statistical parameters are derived from measurements of surface elevation (z), conducted with altimeters that are characterised by their accuracy (δ) and sampling distance (Δx). As usual, the measurement method has an impact on the outcome, in this case by filtering out both high and low spatial frequencies of the surface roughness. Sea ice elevation measurements obtained from airborne altimeters (Ketchum, 1971; Dierking, 1995) and supplemented with terrestrial laser scanners (Landy et al., 2015) draw a picture of sea ice roughness as a multi-scale feature covering several orders of magnitude from large floes and pressure ridges (from tens to hundreds of metres) to frost flowers and small ripples (on the centimetre to millimetre scale). The incident radiation of a given wavelength (λ) reacts differently with individual components of the superimposed roughness of many scales (Ulaby et al., 2014, p. 252).

At small end of the roughness spectrum i.e. when the change of surface elevation over sampling distance (Δz∕Δx) is much smaller than λ ($\mathrm{\Delta }z/\mathrm{\Delta }x\ll \mathit{\lambda }$), the surface roughness is negligible. It means that the angular characteristics of scattered radiation are the same as for the secular surface. As a rule of thumb, Δx should be smaller than 0.1λ (Dierking, 2000). Sea ice roughness measurements with terrestrial lidar carried out by Landy et al. (2015) show that standard deviation of surface height σ_z ranges from 0.001 to 0.0064 m, after high-pass filtering (cut-off at 4 m⁻¹, Δx=0.002 m). These results indicate that, according to the Fraunhofer's smoothness criterion (${\mathit{\sigma }}_z<\mathit{\lambda }/\mathrm{32}\cos \mathit{\theta }$, where θ is the angle of incidence), most sea ice types (except artificially grown frost flowers) can be treated as a smooth surface for the L-band at scales lower than 0.25 m.

In this study, we focus on the other side of the roughness spectrum, i.e. the large-scale surface roughness of sea ice ($\mathrm{\Delta }z/\mathrm{\Delta }x\gg \mathit{\lambda }$). In this case, changes in surface elevation are not negligible and alter the local incidence angle (θ_i). Studies of surface scattering by Lawrence et al. (2011, 2013) conclude that region of 8λ×8λ is sufficient to model small-scale roughness. Here, we assume that at larger spatial scales (larger than 8λ×8λ) the surface roughness can be characterised in terms of geometrical optics (GO), for sea ice with ϵ_ice=4.1, ${\mathit{\lambda }}_{\text{ice}}=\mathit{\lambda }/\sqrt{{\mathit{ϵ}}_{\text{ice}}}\approx \mathrm{0.1}$ m.

GO approximation describes the surface as a set of facets (Ulaby et al., 2014, pages 564–567). This approach was applied for modelling the effective emissivities of mountainous terrain (Matzler and Standley, 2000) and ocean surface (Prigent and Abba, 1990). The last study involved probability distribution of slopes in crosswind and downwind directions. A similar method was used in the context of sea ice to assess the uncertainties caused by roughness in sea ice concentration products derived from passive microwaves (Stroeve et al., 2006). Liu et al. (2014) measured ice surface slopes and other roughness statistics in the Bohai Sea. Their result was obtained with linear (1-D) scans under the assumption of isotropic roughness characteristics. The study by Beckers et al. (2015) has demonstrated that the statistics of sea ice roughness (mean z, σ_z, kurtosis, and skewness) obtained from 1-D altimeter and 2-D laser scanner converge, on the condition that the surface is not strongly heterogeneous. Nonetheless, the 1-D altimeter data cannot properly represent the spatial orientation of surface facets. The surface facet orientation is characterised by both the slope (α) and the azimuthal angle in which it is facing (γ).

In this paper, we address the knowledge gap regarding the influence of large-scale surface roughness on L-band T_B. The paper comprises four main sections. The introduction, presenting the context, is in Sect. 1.

Section 2 introduces the experimental data collected during the SMOSice 2014 campaign (Sect. 2.1). Among them are the EMIRAD2 L-band radiometer T_B measurements, which will serve as a reference for the T_B simulations. Another instrument is the airborne laser scanner (ALS) used for surface elevation measurements. The surface elevation measurements are used to construct a digital elevation model (DEM) of sea ice surface. From the DEM we derive the facet surface slopes and their orientation. In Sect. 2.2 we analyse the statistics of the facet orientation. Based on facet orientation statistics, we derive a parameterisation of the probability distribution function of surface slopes (PDF_α), which will serve as surface roughness representation in T_B simulations.

Section 2.3 presents the setup used for sea ice T_B simulation. For the simulation of the sea ice T_B we use the MIcrowave L-band LAyered Sea ice emission model (MILLAS) (Maaß et al., 2013). In Sect. 2.4 we show how we integrate the surface roughness statistics (PDF_α) with MILLAS using geometrical optics.

The three key results of this study, namely that (a) surface roughness reduces the polarisation difference, this change being most pronounced at incidence angles greater than 50^∘; (b) nadir T_B is little affected; and (c) comparison with the radiometer data and sensitivity study indicates that snow cover has a greater impact on the T_B than surface roughness, are subsequently presented and discussed in Sect. 3.

Section 4 summarises the results and discusses the implications of this study for current and future L-band missions.

2 Materials and methods

In this section we present the SMOSice 2014 campaign that is the key dataset of this study (Sect. 2.1). Section 2.2 presents the sea ice surface roughness measurements in the context of geometrical optics. Section 2.3 presents the sea ice emission model that we used.

Figure 1 The region of the SMOSice 2014 campaign. (a) Sea ice thickness on 24 March 2014, derived from SMOS. The SMOS sea ice thickness product with a resolution of 40 km is presented on a 15 km grid. An aggregation of thick ice (> 1 m) is visible along the Edgeøya's eastern coast. (b) The TerraSAR-X wide-swath mode (HH polarisation), with frames taken at 05:35 and 14:58 UTC. The aircraft tracks are marked in green (A) at 10:05–10:41 UTC and red (B) at 11:25–12:07 UTC.

2.1 SMOSice 2014 campaign

The SMOSice 2014 campaign took place between 21 and 27 March 2014 in the area between Edgeøya and Kong Karls Land, east of Svalbard. Hendricks et al. (2014) and Kaleschke et al. (2016) described the campaign extensively. In this study we analysed the data acquired during the flights on 24 and 26 March. From this point onwards, we focus solely on the parts relevant to the presented work.

In the period preceding the experiment from late January until early March the meteorological conditions in the region deviated strongly from the climatological means. The air temperature measured at Hopen Island meteorological station was on average 9 to 12 ^∘C higher than the climatological value for the period 1961–1990 (Strübing and Schwarz, 2014). Prevailing southerly winds pushed sea ice against the coasts of Nordaustlandet and into Hinlopen Strait, leaving a small strip of compacted ice along the coasts of Edgeøya. When sea ice returned with southerly drift in early March, the scene was set for the experiment. The thickest, most deformed ice was located in the western part of the studied region with a gradual decrease in thickness eastwards, where thin newly formed ice was dominant. This pattern can be observed in the SMOS sea ice thickness product displayed in Fig. 1. In this work, we focus on the data from the low-altitude flight at 70 m, because they are the data with highest spatial resolution of the airborne laser scanner (ALS) among all the flights. The analysis was further reduced to 24 March, due to the fact that the region covered on 26 March had a discontinuous ice cover and a large-scale swell was interfering with the surface elevation measurements. On 24 March, the Polar 5 research aircraft of the Alfred Wegener Institute (Bremerhaven, Germany) undertook measurement flights starting from the eastern coast of Edgeøya, along the lines marked in red and green in Fig. 1. The figure also shows TerraSAR-X wide-swath scenes, taken in the region. Flight A took place between 10:05 and 10:41 UTC, and flight B occurred from 11:25 to 12:07 UTC. The set of instruments mounted on the aircraft included an aerial camera to visually register the ice conditions, the Heitronics KT19.85 pyrometer for surface temperature measurements, the L-band radiometer EMIRAD-2, and the Airborne Laser Scanner (ALS) for high-resolution surface elevation measurements.

2.1.1 EMIRAD-2 radiometer

The EMIRAD-2 L-band radiometer (developed by DTUSpace) is a fully polarimetric system with advanced radio frequency interference (RFI) detection features (Søbjaerg et al., 2013). The setup mounted on the aircraft consists of two Potter horn antennae, one nadir pointing and one side looking at 45^∘ incidence angle. As the aircraft flew at the altitude of 70 m, the footprint of each antenna was 60 m × 90 m for the nadir-pointing antenna and 70 m × 90 m for the side-looking antenna. The receiver's sensitivity is 0.1 K for a 1 s integration window. Along with the L-band measurements, navigation data were collected to transform the polarimetric brightness temperature into the Earth reference frame (Hendricks et al., 2014). The EMIRAD-2 data were screened by evaluating kurtosis, polarimetric (Balling et al., 2012), and brightness temperature (T_B) anomalies. The screening showed RFI contamination for up to 30 % of the samples. After subtracting the mean value of the RFI-flagged data from the mean value of the full data, we found a 10 K difference for the nadir-looking horn, while the difference was non-significant for the side-looking horn. The analysis of the T_B during the “wing wags” calibration manoeuvers further revealed a 20 K offset relative to the nadir-looking vertical channel caused by a continuous wave signal from the camera that was mounted on the aeroplane to obtain visual images. The analysis concludes a purely additive characteristic that allowed for bias correction (Hendricks et al., 2014). In this study, we use the data preprocessed by the DTU team. The radiometer data were RFI-cleaned and bias-corrected and validated using aircraft wing wags and nose wags over open ocean (Hendricks et al., 2014).

2.1.2 Airborne laser scanner

In this study, the ALS (Riegel VQ-580 laser scanner) has two purposes: (1) to measure the surface elevation for subsequent estimation of the ice thickness and (2) to characterise the surface topography. The ALS near-infrared laser (wavelength 1064 nm) measures snow and ice elevation with the accuracy and precision of 0.0025 m. Across-track and along-track elevation measurements were obtained every 0.25 and 0.50 m, respectively. These sampling characteristics resulted from the combination of the flight altitude (70 m) and the setup of the ALS (pulse repetition rate of 50 kHz, cross-track range of ±30^∘). The field of view of the radiometer side-looking antenna was only partially covered by the ALS scans. Nonetheless, we assume that the roughness characteristics measured by the ALS are representative for both antenna fields of view. The data were calibrated and geo-referenced to the WGS84 datum. Further processing involved manual classification of tie points in leads in order to obtain local sea level and sea ice freeboard (Hendricks et al., 2014). The geo-referenced surface elevations are used to compute surface roughness statistics. The elevation data are interpolated to a regular grid of 0.5 m by 0.5 m to form a digital elevation model (DEM) of the sea ice surface. The DEM serves to derive surface slope orientation.

The estimate of sea ice thickness was built on the hydro-static equilibrium assumption. The data required to estimate sea ice thickness consists of (a) the densities of water and ice, (b) the snow load classically described by snow density and snow thickness, and (c) ALS's freeboard data.

The water, ice, and snow densities retained are 1027 kg m⁻³ (water) and 917 kg m⁻³ after Ricker et al. (2014) and 300 kg m⁻³ after Warren et al. (1999).

Snow thickness was meant to be provided by the onboard snow radar; however the equipment was still in the test phase at the time of the experiment. As a workaround, we followed Kaleschke et al. (2016) and used the approximations found in Yu and Rothrock (1996) and Mäkynen et al. (2013) that set the snow thickness to 10 % of the sea ice thickness.

2.2 Sea ice surface roughness

In this subsection we will analyse the data from the airborne laser scanner (ALS) that we presented in Sect. 2.1.2. We use the ALS data to measure the decimetre-scale surface roughness. The ALS is a laser instrument that measures the distance to the surface. If snow covers the ice the ALS will register the snow–air interface as the elevation. Therefore, the relief of the ice is modified by snow cover. During the SMOSice 2014 campaign the snow measurements were unavailable. We assume that snow cover is a plane-parallel layer over sea ice. This assumption does not account for snow dunes and drifts that may form on the ice. The implications of snow thickness on the radioactive transfer modelling are discussed in the section presenting the sensitivity analysis (Sect. 3.2).

In the context of radiation transfer, the surface roughness is characterised in relation to incident wavelength. The ALS along-track spatial sampling of 0.5 m is a few times larger than the L-band wavelength in sea ice (${\mathit{\lambda }}_{\text{ice}}=\mathit{\lambda }/\sqrt{{\mathit{ϵ}}_{\text{ice}}}\approx \mathrm{0.1}$ m), which makes it suitable to measure the large-scale roughness, the part of the roughness spectrum where GO can be used to approximate the path of radiation. The analysis of ALS elevation data is done in three steps.

In the first step, we identify the ice with different degree of surface roughness. For that purpose we divide the flight tracks into 1 s sections (approximately 70 m long), large enough to cover the entire nadir radiometer footprint, and we build a histogram of the standard deviations of surface heights computed for these sections. The number of bins in the histogram is set according to the formula N_bins=5log ₁₀(N_data), after Panofsky and Brier (1958). We chose standard deviation as a criterion for defining the roughness classes, as it is widely used to characterise surface roughness from elevation profiles. Also, unlike visual interpretation of the aerial photography of sea ice, it does not introduce personal biases. The resulting histogram in Fig. 2 shows the sea ice roughness classes as histograms bins. No sections within the lowest standard deviation of surface height were found. That is probably due to the fact that no refrozen lead of the scale of 70 m was found or the ALS laser signal was not reflected back from the surface resulting in missing data.

Figure 2 Histogram of the standard deviation of surface heights computed from 70 m flight strips; bins define the roughness classes of sea ice. Examples for the three roughness classes “smooth”, “medium rough”, and “rough” are marked in the colours blue, green, and red, respectively.

Figure 3 The f_R parameter illustrating the deviation from the uniform azimuth distribution along 1000 randomly selected samples. The average value of f_R is marked as a thick red line. The considered threshold of uniform distribution, 4.3 km, is marked by the blue dashed line.

In the second step, we interpolate the ALS elevation measurements to a regular 0.5 m grid in order to form a digital elevation model (DEM) of the sea ice surface. The sea ice surface in the DEM is represented as a set of triangular facets. Each facet orientation in the 3-D Cartesian space (for simplicity we assume the base vectors $\widehat{\mathit{x}},\widehat{\mathit{y}},\widehat{\mathit{z}}$ to be aligned with the aircraft principal axis, so $\widehat{\mathit{y}}$ points to the flight direction) is described by two angles: facet slope α ($\mathrm{0}\le \mathit{\alpha }<\mathit{\pi }$) and facet azimuthal direction γ ($-\mathit{\pi }\le \mathit{\gamma }<\mathit{\pi }$). Therefore, the ith facet local normal vector is described by

$$\begin{array}{}\text{(1)}& {\widehat{\mathit{n}}}_i=-\widehat{\mathit{x}}\sin \left({\mathit{\alpha }}_i\right)\cos \left({\mathit{\gamma }}_i\right)-\widehat{\mathit{y}}\sin \left({\mathit{\gamma }}_i\right)\sin \left({\mathit{\alpha }}_i\right)+\widehat{\mathit{z}}\cos \left({\mathit{\alpha }}_i\right).\end{array}$$

In the third step, we compute the normal vectors and their orientations for the individual facets. This is done for all roughness classes. We found that the azimuthal orientation angle γ does not show any preferred directions within any given roughness class. In the next subsections we present the analysis of the two angles characterising the facet: azimuthal direction and surface slope.

2.2.1 Facet azimuth orientation

In the previous section we used the DEM to calculate the vectors normal to the surface facets. In this subsection we analyse the orientations of facet azimuths. In order to evaluate the distribution of the facet azimuths, we define parameter f_R (Eq. 2). This parameter is calculated from a histogram of azimuth orientation. In Eq. (2) the denominator is the total number of measurements expressed as the number of angular bins multiplied by the mean number of samples per bin (total number of samples). The numerator is a sum of the differences between the mean number of samples per bin and the actual number of samples in each bin. There are 36 bins. The number of bins was determined with N_bins=5log ₁₀(N_data) considering the maximal number of samples in the 25 km flight track. The f_R parameter equals zero for the perfectly uniform distribution, in which case the number of counts in each bin (n_i) equals a mean number of counts (μ). The f_R parameter reaches its maximum value of $f_{R_{\text{max}}}=\mathrm{2}-\mathrm{4}/N_{\text{bins}}$ when the slopes are aligned, i.e. grouped in two bins.

$$\begin{array}{}\text{(2)}& \begin{array}{rl}& f_R={\displaystyle \frac{{\sum }_i^{N_{\text{bins}}}\left(\right|n_i-\mathit{\mu }\left|\right)}{N_{\text{bins}}\mathit{\mu }}}\\ & \mathit{\mu }={\displaystyle \frac{\mathrm{1}}{N_{\text{bins}}}}\sum _i^{N_{\text{bins}}}{n}_i\end{array}\end{array}$$

To evaluate f_R we selected 1000 random 15 km samples from the flight tracks. The analysis of the samples shows that the deviation from the uniform distribution decreases sharply with increasing distance over the first kilometre (Fig. 3). For distances along the flight track greater than 4.3 km the curve flattens at a value of f_R=0.05 in 90 % of the samples. We assume that at scales greater than 4.3 km (marked by the vertical dashed line in Fig. 3) slope orientations do not have a preferential direction beyond natural variability. This distance corresponds to an approximately 60 s section. In Fig. 3 the average value of f_R is marked as a thick red line. Several sample profiles are plotted as grey lines to illustrate the variability.

2.2.2 Facet slope angle

Section 2.2.1 looked at the azimuthal orientation of surface facets. This section focuses on the analysis of facet slopes. For all roughness classes we observe a similar probability density function (PDF) of surface slopes. The PDFs have a maximum at zero and a gradual decline in the likelihood of encountering the steeper slopes. Figure 4 shows the PDF_α on a logarithmic scale for the three distinct roughness classes: smooth 0.05 m < σ_z < 0.11 m (in blue), medium rough 0.28 m < σ_z < 0.34 m (in green), and rough 0.45 m < σ_z < 0.51 m (in red).

We decide to approximate the PDF of surface slopes with an exponential curve:

$$\begin{array}{}\text{(3)}& {\text{PDF}}_{\mathit{\alpha }}=C_{\text{norm}}\exp (-(\mathit{\alpha }/s_{\mathit{\alpha }}\left)\right),\end{array}$$

where C_norm is the normalisation constant and s_α is the “geometrical-slope roughness parameter”. Figure 4 presents the data and the exponential approximations. The log scale is very relevant because it becomes obvious that the chances of encountering steep slopes are getting slimmer the higher the slope angle. Consequently, it makes sense that the approximation functions misfit the observations at high slope angles as it is irrelevant to fit an approximation there. As s_α is the only parameter of the approximation function, it is descriptive of the surface roughness.

Figure 5 shows the relation between s_α and the standard deviation of surface heights corresponding to the roughness classes defined above. The error bars represent the uncertainty associated with each data point. Very rough ice has a large uncertainty because the number of samples was small (classes with σ_z>0.6 m accounted for less than 75 s of flight, out of a total 78 min). The quadratic relation is holding well for ice with up to 0.5 m standard deviation of surface heights. The equation of the fitted curve is $s_{\mathit{\alpha }}=\mathrm{51.61}{\mathit{\sigma }}_z^{\mathrm{2}}+\mathrm{1.50}{\mathit{\sigma }}_z+\mathrm{0.14}$.

Figure 4 Density of probability of surface slopes on a logarithmic scale for three roughness classes: smooth 0.05 m < σ_z < 0.11 m (in blue), medium rough 0.28 m < σ_z < 0.34 m (in green), and rough 0.45 m < σ_z < 0.51 m (in red). The exponential fits to the respective curves are marked as thin coloured dotted lines.

Figure 5 Surface roughness parameter s_α describing the probability distribution of surface slopes. Error bars are inversely proportional to the number of data points in each roughness class. The “smooth”, “medium rough”, and “rough” classes are marked in the colours blue, green, and red, respectively. The red dashed line marks the fitted curve.

2.3 Sea ice brightness temperature simulation

In this subsection we present the emission model for simulating the sea ice brightness temperature (T_B). We use the MIcrowave L-band LAyered Sea ice emission model described by Maaß et al. (2013), further referred to as MILLAS. This model is based on the radiative transfer model of Burke et al. (1979) (who used it for soils), with infinite half-space of seawater covered with layers of sea ice, snow, and a top semi-infinite layer of air. In contrast to the original model of Burke et al. (1979) and its usage by Maaß et al. (2013), the current version of MILLAS takes into account multiple reflections at the layer boundaries. The multiple reflections are expressed as subsequent terms of a geometric series. The summation over series stops when terms of the series contribute less than a given threshold to the total (in the following calculations the threshold was set to 0.001). MILLAS describes the brightness temperature above snow-covered sea ice as a function of temperature and permittivity of the layers. The water permittivity depends mainly on the water temperature and salinity (Klein and Swift, 1977). Ice permittivity can be approximately described as a function of brine volume fraction (Vant et al., 1978), which depends on ice salinity and the densities of the ice and brine (Cox and Weeks, 1982), which in turn depends mainly on ice temperature. We set the ice salinity to 4 g kg⁻¹ which is a mean value for first-year ice determined by Cox and Weeks (1974). The permittivity of dry snow can be estimated from its density and temperature (Tiuri et al., 1984). In the simulation, the ice and water salinity are kept constant (see Table 1). Furthermore, we assume that the system is in thermal equilibrium and that the water beneath the ice is at the freezing point. In this configuration, the T_B is simulated as a function of ice thickness (d_ice), snow thickness (d_snow), and surface temperature (T_surf). In our setup, the snow is assumed to be dry with a density of 300 kg m⁻³, which is the average snow density value for December Arctic measurements from 1954 to 1991 (Warren et al., 1999). The T_B simulations are only slightly sensitive to snow density; see Fig. 3 in Maaß et al. (2013). The permittivities of snow and ice are linked to their temperature. The temperature profiles within snow and ice are assumed to be continuous and linear. The values for the ice and snow thermal conductivity are taken from Yu and Rothrock (1996) and Untersteiner (1961). As the optimisation of the emission model lies beyond the scope of this work, we use the simplest setup variant of MILLAS consisting of a single layer of ice covered with a single layer of snow.

Table 1 Brightness temperature simulation setup of the MIcrowave L-band LAyered Sea ice emission model (MILLAS).

2.4 Simulation of brightness temperature of rough sea ice

In the previous sections, we described the sea ice surface as composed of facets with an orientation described by two angles: the slope α and the azimuthal direction γ. Subsequently, we analysed the ALS data to extract information about statistical distributions of slopes and their orientation. We concluded that the exponential function is suitable to describe the probability density function of surface slopes for a range of ice surfaces with different degrees of surface roughness.

In this section, we describe how we integrate the probability description of faceted sea ice surface with the MILLAS emission model. We will start by describing the coordinate system that we used in the T_B simulations. The relations between radiometer antenna-look direction ($\widehat{\mathit{r}}$) and the horizontal ($\widehat{\mathit{h}}$) and vertical ($\widehat{\mathit{v}}$) polarisation vectors are described in a Cartesian coordinate system. We show how these relations are represented in the coordinate system associated with a facet. The vectors defined in the tilted facet coordinate system are denoted with the subscript i. Subsequently, we will derive the equation that sums up the emissions originating from multiple facets.

We consider the Cartesian coordinate system ($\widehat{\mathit{x}}$, $\widehat{\mathit{y}}$, $\widehat{\mathit{z}}$) with the origin in the centre of the nadir sensor footprint. In this reference frame the radiometer antenna-look direction ($\widehat{\mathit{r}}$) is described as

$$\begin{array}{}\text{(4)}& \widehat{\mathit{r}}=\widehat{\mathit{x}}\sin \left({\mathit{\theta }}_{\mathrm{0}}\right)\cos \left({\mathit{\varphi }}_{\mathrm{0}}\right)+\widehat{\mathit{y}}\sin \left({\mathit{\theta }}_{\mathrm{0}}\right)\sin \left({\mathit{\varphi }}_{\mathrm{0}}\right)-\widehat{\mathit{z}}\cos \left({\mathit{\theta }}_{\mathrm{0}}\right),\end{array}$$

where the θ₀ is the antenna incidence angle and the ϕ₀ is the azimuthal direction of the antenna. In this particular case we set the reference system so ϕ₀=0 and $\widehat{\mathit{x}}$ is parallel to the ground. The antenna setting defines the directions of the horizontal ($\widehat{\mathit{h}}$) and vertical ($\widehat{\mathit{v}}$) polarisation vectors:

$$\begin{array}{}\text{(5)}& \begin{array}{rl}& \widehat{\mathit{h}}=-\widehat{\mathit{x}}\sin \left({\mathit{\varphi }}_{\mathrm{0}}\right)+\widehat{\mathit{y}}\cos \left({\mathit{\varphi }}_{\mathrm{0}}\right),\\ & \widehat{\mathit{v}}=-\widehat{\mathit{x}}\cos \left({\mathit{\theta }}_{\mathrm{0}}\right)\cos \left({\mathit{\varphi }}_{\mathrm{0}}\right)-\widehat{\mathit{y}}\cos \left({\mathit{\theta }}_{\mathrm{0}}\right)\sin \left({\mathit{\varphi }}_{\mathrm{0}}\right)-\widehat{\mathit{z}}\sin \left({\mathit{\theta }}_{\mathrm{0}}\right).\end{array}\end{array}$$

We are interested in finding a relationship between the radiation originating from a tilted facet and from a flat one. For that purpose we must consider the tilted coordinate system associated with the ith facet (variables associated with individual facets are denoted with the subscript i). The z coordinate in this tilted coordinate system is aligned with the facet normal vector ${\widehat{\mathit{n}}}_i$, followed by ${\widehat{\mathit{x}}}_i$ and ${\widehat{\mathit{y}}}_i$ calculated accordingly.

$$\begin{array}{}\text{(6)}& \begin{array}{rl}& {\widehat{\mathit{z}}}_i={\widehat{\mathit{n}}}_i\\ & {\widehat{\mathit{y}}}_i={\displaystyle \frac{{\widehat{\mathit{n}}}_i\times \widehat{\mathit{r}}}{|{\widehat{\mathit{n}}}_i\times \widehat{\mathit{r}}|}}\\ & {\widehat{\mathit{x}}}_i={\widehat{\mathit{y}}}_i\times {\widehat{\mathit{z}}}_i\end{array}\end{array}$$

Therefore, the local incidence angle θ_i is

$$\begin{array}{}\text{(7)}& {\mathit{\theta }}_{\text{i}}={\cos }^{-\mathrm{1}}(-\widehat{\mathit{r}}\cdot {\widehat{\mathit{n}}}_i),\end{array}$$

and the local horizontal and vertical polarisation vectors are

$$\begin{array}{}\text{(8)}& \begin{array}{rl}& {\widehat{\mathit{h}}}_i={\widehat{\mathit{y}}}_i,\\ & {\widehat{\mathit{v}}}_i=-{\widehat{\mathit{x}}}_i\cos \left({\mathit{\theta }}_{\text{i}}\right)-{\widehat{\mathit{z}}}_i\sin \left({\mathit{\theta }}_{\text{i}}\right).\end{array}\end{array}$$

The emissions from the facet at an angle θ_i and polarisation p are denoted with an asterisk: $T_{\text{B}}^{*}({\mathit{\theta }}_{\text{i}};p)$. In order to calculate the brightness temperatures observed in the global horizontal and vertical polarisation, it is necessary to account for the coordinates' rotation (Ulaby et al., 2014, p. 443):

$$\begin{array}{}\text{(9a)}& T_{\text{Bi}}({\mathit{\theta }}_{\text{i}};H)=(\widehat{\mathit{h}}\cdot {\widehat{\mathit{h}}}_i{)}^{\mathrm{2}}{T}_{\text{B}}^{*}({\mathit{\theta }}_{\text{i}};H)+(\widehat{\mathit{h}}\cdot {\widehat{\mathit{v}}}_i{)}^{\mathrm{2}}{T}_{\text{B}}^{*}({\mathit{\theta }}_{\text{i}};V),\end{array}$$

$$\begin{array}{}\text{(9b)}& T_{\text{Bi}}({\mathit{\theta }}_{\text{i}};V)=(\widehat{\mathit{v}}\cdot {\widehat{\mathit{h}}}_i{)}^{\mathrm{2}}{T}_{\text{B}}^{*}({\mathit{\theta }}_{\text{i}};H)+(\widehat{\mathit{v}}\cdot {\widehat{\mathit{v}}}_i{)}^{\mathrm{2}}{T}_{\text{B}}^{*}({\mathit{\theta }}_{\text{i}};V).\end{array}$$

We model the sea ice surface as a set of facets. Therefore the brightness temperature registered at the antenna aperture is a sum of contributions from N facets. We assume that each facet area A at the distance R is visible at the incidence angle θ_i and covers a patch of antenna field of view, equal to the solid angle Ω_i:

$$\begin{array}{}\text{(10)}& {\mathrm{\Omega }}_i={\displaystyle \frac{A\cos \left({\mathit{\theta }}_{\text{i}}\right)}{R^{\mathrm{2}}\cos \left({\mathit{\alpha }}_i\right)}}.\end{array}$$

The formula summing the contributions from N visible facets is

$$\begin{array}{}\text{(11)}& T_{\text{B}}({\mathit{\theta }}_{\mathrm{0}};p)={\displaystyle \frac{\mathrm{1}}{N\cos \left({\mathit{\theta }}_{\mathrm{0}}\right)}}\sum _{i=\mathrm{1}}^N{\mathit{\omega }}_i{T}_{\text{Bi}}({\mathit{\theta }}_{\text{i}};p){\mathrm{\Omega }}_i,\end{array}$$

with ω_i the antenna gain component.

At this stage of our study, we aim at modelling the effect of surface roughness on the T_B. As a first-order approximation we assume that the antenna gain is constant across the whole field of view and that the antenna is in a far field so the incidence angle θ₀ is assumed constant across the field of view. This assumption is more suitable for the space-borne system. The resulting formula is

$$\begin{array}{}\text{(12)}& T_{\text{B}}({\mathit{\theta }}_{\mathrm{0}};p)={\displaystyle \frac{\mathrm{1}}{N\cos \left({\mathit{\theta }}_{\mathrm{0}}\right)}}\sum _{i=\mathrm{1}}^N{T}_{\text{Bi}}({\mathit{\theta }}_{\text{i}};p)\sec \left({\mathit{\alpha }}_i\right)(-\widehat{\mathit{r}}\cdot {\widehat{\mathit{n}}}_i).\end{array}$$

This is the formula that we implement in the geometrical roughness model. Figure 6 presents the data flow in the geometrical roughness model. The model merges the MILLAS emission model (suitable for sea ice) with the geometrical characterisation of the faceted surface. In the presented setup the MILLAS emission model uses the sea ice surface temperature (T_surface), sea ice thickness (d_ice), and snow thickness (d_snow) as input variables. The GO part needs the cumulative probability distribution of surface slopes (CDF_α) and the antenna look direction ($\widehat{\mathit{r}}$). The orientation of N facets representing the sea ice surface is calculated with the inverse transform sampling (ITS) (Devroye, 2006). This method returns a random slope value from a given non-uniform distribution. The non-uniform distribution is described by a cumulative probability distribution, which in this case depends on the geometrical roughness parameter s_α. Similarly, the azimuthal angle γ_i is drawn from a uniform distribution. The result of this processing step is the set of N pairs of angles (α_i, γ_i) describing the orientation of N facets. In the next step, the local normal vector and the local incidence angle (θ_i) are calculated for each of the N facets. The θ_i is used to calculate the brightness temperature emitted from the ith facet with the emission model. Shadowing occurs when ${\mathit{\theta }}_{\text{i}}>\mathit{\pi }/\mathrm{2}$ and the radiation is emitted away from the antenna. In the current setup the double-bounce effects are not accounted for. $\widehat{\mathit{r}}$ and ${\widehat{\mathit{n}}}_i$ are used to calculate the local and global polarisation directions, as well as Ω_i. For the final step, the summed contributions from N facets are summed as in Eq. (12). The result is the brightness temperature of the surface observed under θ₀ and polarisation p.

Figure 6 Flow chart showing the processing chain within the geometrical roughness model. The model consists of two principal blocks: the emission model and the inverse transport sampling ITS module responsible for generating the facet orientation. The input parameters surface temperature T_surface, ice thickness d_ice, snow thickness d_snow, cumulative distribution of surface slopes CDF_α, and antenna direction $\widehat{\mathit{r}}$ are colour coded in red.

Figure 7 Effects of the large-scale surface roughness on the brightness temperature of sea ice, simulated with the geometrical roughness model. Vertical polarisation is in red, and horizontal polarisation is in blue. The black dashed lines mark the T_B for the specular surface. The surface roughness parameter s_α varies between 2 and 20; the thicker the line the higher the s_α. The inputs for the MILLAS emissivity model are kept constant: T_surface=260 K, d_ice=1.42 m, and d_snow=0.14 m.

The number of facets N impacts the quality of the simulation and the reproducibility of the results. In order to set the value of N, we looked at the standard deviation of 20 T_B simulation results for nadir and 45^∘. We decided that the standard deviation of T_B should be lower than 0.1 K, as this is the accuracy of the EMIRAD-2 radiometer for the 1 s integration time. This criterion is met for N greater than 10⁴, which we take as the N value for further simulations.

3 Results

The three main results of this study are as follows. (a) Surface roughness reduces the polarisation difference; this change is most pronounced at incidence angle greater than 50^∘. (b) Nadir T_B is little affected. (c) Comparison with the radiometer data and sensitivity study indicates that snow cover has a greater impact on the T_B than surface roughness.

In Sect. 3.1 we show the brightness temperature simulations for sea ice with different degrees of large-scale surface roughness. To interpret the results we conduct a sensitivity analysis (Sect. 3.2). The comparison of the simulated vs. measured T_B over 4.3 km flight track samples is shown in Sect. 3.3.

3.1 Brightness temperature simulations

In Sect. 2.2.2 we derived a parameterisation of the degree of surface roughness. We approximate the roughness by an exponential probability distribution function (PDF) of surface slopes. The shape of the PDF is fully described by the s_α parameter.

In our simulation s_α varies between 2 and 20 in accordance with the measurements of surface slopes conducted during SMOSice 2014. As the aim of this study is to characterise the effect of surface roughness on L-band T_B, in this section we keep the other parameters of the emissivity model constant: surface temperature and ice and snow thickness (T_surface=260 K, d_ice=1.42 m, d_snow=0.14 m). The brightness temperatures are calculated for every degree of incidence angle in the range 0–70^∘. Figure 7 shows the simulation results.

The effect of increasing surface roughness is two-fold. First, the overall near-nadir intensity is lowered by 2.6 K. Second, the polarisation difference decreases. For the highest value of the roughness parameter, at Brewster's angle the vertical polarisation decreases by 8 K and horizontal polarisation experiences a 4 K increase. The effect of roughness is more pronounced for larger values of roughness parameter s_α and is most visible at higher incidence angles (60^∘).

The polarisation mixing can be explained by the approach used in this study. The emissions from the facet in horizontal (${\widehat{\mathit{h}}}_i$) and vertical (${\widehat{\mathit{v}}}_i$) polarisations are partially mixed when expressed in the ($\widehat{\mathit{h}},\widehat{\mathit{v}}$) coordinate system (see Eq. 12).

The lowering of the intensity has two possible explanations. First is that our model does not take into account shadowing effects. When local incidence angle is greater than 90^∘, the facet is emitting away from the antenna. For the “near-nadir” angles (0–30^∘), the likelihood of shadowing is less than 1 % for the most rough ice (see Fig. 4).

The second explanation of this effect is associated with shape of the Fresnel emissivity curves. The polarisation difference for large incidence angles is larger than for the near-nadir ones. Therefore, the mean of the two polarisations ($T_{\text{B}}({\mathit{\theta }}_{\mathrm{0}};H)+T_{\text{B}}({\mathit{\theta }}_{\mathrm{0}};V))/\mathrm{2}$ i.e. total intensity) is fairly constant up to 30^∘ and than drops continuously ∼ 10 K by 60^∘. The trend continues for higher incidence angles. High values of s_α increase the likelihood of returning a large incidence angle in the inverse transform sampling. These large incidence angles contribute to the overall lowering of the T_B. The contribution of this mechanism is ≈ 2 K for T_B(0) in the case of the most rough ice. The 2 K estimate was obtained by integrating the drop in total intensity weighted by the PDF_α.

The above results are obtained with a Monte Carlo simulation. This method is a time-consuming approach. Therefore, we propose a parameterisation of the simulation results. The two effects, the polarisation mixing and the lowering of brightness temperature, can be expressed in a fashion similar to the HQ model proposed by Choudhury et al. (1979). Here we propose a formulation with two parameters H_α and Q_α.

$$\begin{array}{}\text{(13)}& T_{\text{B}}(\mathit{\theta };p)=\left[\right(\mathrm{1}-Q_{\mathit{\alpha }})\cdot T_{\text{B}}^{*}(\mathit{\theta };p)+Q_{\mathit{\alpha }}{T}_{\text{B}}^{*}(\mathit{\theta };q\left)\right]\cdot H_{\mathit{\alpha }},\end{array}$$

where p and q stand for the polarisation.

H_α accounts for the change in total intensity and Q_α for the polarisation mixing. The emissions from the specular surface are denoted with an asterisk: $T_{\text{B}}^{*}(\mathit{\theta };p,q)$. The proposed parameterisation approximates the results obtained with the Monte Carlo simulation with a root-mean-square difference of 0.45 K:

$$\begin{array}{}\text{(14)}& \begin{array}{rl}& H_{\mathit{\alpha }}=a_{\mathrm{1}}{s}_{\mathit{\alpha }}^{\mathrm{2}}+a_{\mathrm{0}}\\ & Q_{\mathit{\alpha }}=b_{\mathrm{1}}{s}_{\mathit{\alpha }}^{\mathrm{2}}+b_{\mathrm{0}}\end{array},\end{array}$$

where $a_{\mathrm{1}}=\mathrm{-}\mathrm{0.018}\times {\mathrm{10}}^{-\mathrm{3}}$, a₀=1 and $b_{\mathrm{1}}=\mathrm{0.532}\times {\mathrm{10}}^{-\mathrm{3}}$, b₀=0.

The emissions from the specular surface are an essential input for the geometrical roughness model used in this study. The exact shape of the simulated brightness temperature curves depends on the probability distribution of slopes as well as on the emission characteristics of the specular surface. In this paragraph, we will investigate how the shape of the $T_{\text{B}}^{*}(\mathit{\theta };p,q)$ influences the geometrical roughness model results. The shapes of the polarisation curves, i.e. the reflectivities for a given incidence angle, are described by the Fresnel equations. Equations that are determined by the permittivity of the medium (ϵ). (In this section we omit the question of penetration depth assuming that the emissions are coming from the isothermal ice layer of constant permittivity.) To investigate the impact of the varying ϵ, we take a range of permittivities specific to sea ice as calculated within the MILLAS model. In the present setup the sea ice permittivity depends on bulk ice salinity and ice temperature. We calculate the ϵ for a range of ice temperatures (250 K < T_ice < 271 K) and salinities (1 g kg⁻¹ $<S_{\text{ice}}<$ 12 g kg⁻¹).

The sea ice permittivities from the MILLAS model range between $\mathit{ϵ}=\mathrm{3.1}+\mathrm{0.05}i$ (for T_ice=271 K, S_ice=7 g kg⁻¹) and $\mathit{ϵ}=\mathrm{4.6}+\mathrm{0.8}i$ (for T_ice=253 K, S_ice=1 g kg⁻¹), where T_ice is the bulk ice temperature and S_ice is the bulk ice salinity. The curves corresponding to those values lie close together, indicating that the proposed parameterisation is suitable for all types of sea ice. The effect of permittivity on the polarisation mixing parameter (Q_α) is less pronounced. The dependence of the Q_α parameter on the roughness follows a similar quadratic curve regardless of the surface permittivity.

3.2 Sensitivity analysis

In this section, we investigate how sensitive the model is to its main variables: surface temperature, ice thickness, snow thickness, surface roughness, and the implicit assumption of 100 % sea ice concentration. This is a mandatory step toward the evaluation of the model (Sect. 3.3).

The two most important factors influencing the L-band brightness temperature over sea ice are the ice concentration and the ice thickness. We calculate the sensitivity of our model with regard to sea ice concentration by assuming a linear mixing of water and thick ice fractions within the radiometer footprint. The brightness temperature of seawater T_Bw (salinity of 33 g kg⁻¹, temperature 271.2 K) is approximately 110 K, and $T_{\text{Bi}}(\mathrm{0};p,q)$ of thick sea ice (T_surf=260 K, d_ice=1.5 m, bulk salinity of 3 g kg⁻¹) is 240 K. The resulting sensitivity to sea ice concentration is ≈ 1.5 K %⁻¹.

The sensitivity of the T_B to sea ice thickness over thin sea ice (d_ice<0.75 m) is fundamental for the sea ice thickness retrieval from L-band radiometry. T_B saturates when the sea ice thickness is significantly larger than the penetration depth of the L-band radiation. Therefore, our analysis focuses on sea ice thicker than 1 m, in order to single out the contributions of surface roughness.

Figure 8 Change in the brightness temperature as predicted by the MILLAS emission model as a function of surface temperature. The different line styles correspond to the different snow thickness assumptions. The calculation was carried out for the sea ice thickness of 1 m and surface roughness set to zero.

Table 2 Table with sensitivities of the brightness temperature at nadir, 45, and 60^∘ as simulated with geometrical optics surface roughness model. The input parameters are roughness parameter s_alpha, ice thickness d_ice, snow thickness d_snow, and surface temperature T_surf.

Table 2 contains the sensitivities of the geometrical roughness model to the input parameters: roughness parameter s_α, ice thickness d_ice, snow thickness d_snow, and surface temperature T_surf. Presented values are grouped into columns corresponding to the polarisation and three incidence angles: 0, 45, and 60^∘. The angles are chosen to reflect the antennae configuration during SMOSice 2014 with an addition of 60^∘, close to Brewster's angle where surface roughness effects are most pronounced.

The assumption about snow thickness has a considerable effect on the sea ice T_B (Maaß et al., 2013). For this reason the values of sensitivities are considered for a number of snow thickness values. Figure 8 shows the simulated L-band brightness temperatures at 0^∘ and 45^∘ as a function of the surface temperature. In this study we make considerable assumptions about snow thickness. To illustrate the uncertainty associated with the assumptions, the simulations are made for a range of snow thicknesses plotted in different line styles.

In the MILLAS model, ice permittivity is parameterised with ice temperature. The non-monotonic shape of the curves is caused by change in ice permittivity. Therefore, the sensitivities are calculated for two temperature ranges: lower (250–265 K) and higher (265–270 K). Table 2 contains the calculated sensitivities.

As far as the large-scale surface roughness is concerned, the sensitivity increases almost linearly for the values of s_α between 0 and 20^∘, which is the maximal value measured during the SMOSice 2014 campaign.

In order to interpret the results of the comparison between the simulation and measurements, it is necessary to evaluate the uncertainties associated with the input parameters: surface temperature, ice thickness, and snow thickness. In the following paragraphs, we use ”mean model sensitivity for the cold conditions” for the averaged absolute sensitivity of $T_{\text{B}}(\mathrm{0};H,V)$ and $T_{\text{B}}(\mathrm{45};H,V)$ at 250 K. We take the values for the lower temperature range as they match the conditions during the SMOSice 2014 campaign.

The uncertainty associated with surface temperature is estimated as the product of the sensitivity of the model times the measurement uncertainty. The surface temperature measurements carried out with the KT19.85 have an accuracy of 0.5 K. The mean surface temperature in the region covered by ice was 251.7 ± 3.5 K. We take the standard deviation of surface temperature measurements as the parameter uncertainty. We estimate the uncertainty associated with surface temperature to be 0.7 K.

The uncertainty associated with sea ice thickness is estimated as the product of the sensitivity of the model times the measurement uncertainty. The sea ice thickness measurements in this study are derived from the resampled ALS elevation data. The mean standard deviation of the resampled elevation measurements is 0.08 m. The assumption about the densities of snow, ice, and water combined with the assumption on the snow thickness of 1∕10 of ice thickness lead to the uncertainty of 0.4 m. Taking into account the mean model sensitivity for the cold conditions prevailing during the flights, we estimate that the uncertainty associated with sea ice thickness is 0.5 K.

The uncertainty associated with snow thickness is estimated as the product of the sensitivity of the model times the measurement uncertainty. Unfortunately, the snow thickness measurements are not available. The snow layer, although transparent for the L-band radiation, is not invisible. The refraction on the snow–ice and snow–air interfaces alters the local incidence angles. Snow cover also has an effect on the temperature profile within the ice. This indirectly affects the permittivity of sea ice. All these factors make an estimation of the uncertainty caused by snow thickness especially hard to quantify. We assume that snow thickness uncertainty is equal to the mean standard deviation of the resampled elevation measurements, that is 0.08 m. The mean model sensitivity to snow thickness for the cold conditions is 8.6 K m⁻¹. Therefore, we estimate the uncertainty associated with snow thickness to be 0.7 K.

An important factor which is not directly included in the model is the sea ice concentration. In the model we assume sea ice concentration to be 100% in order to single out the much smaller contribution of surface roughness. However, if a linear mixing model is applied, the sensitivity to sea ice concentration is −1.5 K %⁻¹. During the preprocessing of the airborne laser scanner (ALS) data we excluded the 70 m sections with more than 5 % missing values. The missing values are caused by the instrument setup (rotating mirror, edge of swath) or by the lack of return reflection from open water or thin ice. We estimate that the uncertainty associated with the sea ice concentration is up to 7.5 K.

To put the partial sensitivities into perspective, the expected changes in the T_B caused by the strongest surface roughness measured during the SMOSice 2014 campaign do not exceed −2.2 K for nadir and 1 K and −5.6 K for the horizontal and vertical polarisation of the 45^∘ antenna, respectively.

To conclude, the sensitivity analysis of the geometrical roughness model leads to the conclusion that the surface roughness effects will be hard to observe in the SMOSice 2014 flight data with the current emission model setup.

Figure 9 Scatter plots illustrating the comparisons between the EMIRAD-2 data and the T_B simulated without GO roughness included – specular. Results corresponding to the setup with snow are marked in green, and without snow is in blue.

Figure 10 Scatter plots illustrating the comparisons between the EMIRAD-2 data and the T_B simulated with GO roughness included. Results corresponding to the setup with snow are marked in green, and without snow is in blue.

3.3 Simulations vs. measurements

In this section, we compare the brightness temperature measured with the EMIRAD-2 radiometer with the brightness temperature simulations. The comparison is done on a 4.3 km section as to justify the assumption of the isotropic azimuth distribution. We want to determine the simulation setup that best reproduces the radiometer measurements and whether the inclusion of the surface roughness in the simulation brings significant improvement. The limitation of this approach is that we assume that the ice observed by the side-looking antenna and the ice below the flight path have the same properties. We consider the surface temperature, the sea ice thickness, and the surface roughness along the flight, and we use them to run the statistical roughness model, described in Sect. 2.4 with the MILLAS single ice layer setup as the brightness temperature module. In setups including snow layer, the snow thickness is set to be 10 % of the sea ice thickness. The calculation is done for 60 s averages, during which the aircraft covers the distance of approximately 4.3 km. At this scale we consider the surface slopes' orientation to be isotropic. The measured surface slopes are used to compute empirical CDF_α. Analysis at the 4.3 km scale is justified by the fact that the antennae gains are not known. Furthermore, the averaged sea ice properties are more representative for the space-borne radiometers.

For each radiometer channel we made four simulation setups: two without roughness (Flat no snow, Flat snow), and two with roughness included (Rough (GO) no snow, Rough (GO) snow). As for the performance metrics of the model setups, we use the coefficient of determination (r²), the root-mean-square error (RMSE), and bias. These metrics are widely used in the assessment of the performance of satellite measurements (Entekhabi et al., 2010). Table 3 contains the results of the comparison expressed in terms of r², RMSE, and bias. The corresponding scatter plots illustrating the comparison between measured and modelled brightness temperatures are presented in Figs. 9 and 10.

The values of r² for all “channel–simulation setup” combinations do not exceed 0.36. The simplified one-layer model managed to capture only 36 % of the signal variance even with surface roughness included. Furthermore, the inclusion of surface roughness brings little improvement to the statistics. In the case of vertical polarisation, where the model studies indicate the most sensitivity to roughness, the r² is even a little lower. The inclusion of a very crude snow thickness parameterisation is more successful in capturing the radiometer measurements' variability. All metrics show that the four model setups perform poorly in reproducing the EMIRAD-2 measurements.

The results of the comparison are also presented in the form of histograms of the difference between the measured and simulated T_B (Fig. 11). For all four antenna feeds the RMSE between simulated and measured T_B decreases whenever the setups include snow. However, the bias increases.

The high values of RMSE and bias show a general misfit of the model to the data. The sensitivity study of the model presented in Sect. 3.2 indicates that the effects of surface roughness are of comparable magnitude or smaller than uncertainties associated with the model input parameters. The side-looking vertical polarisation is predicted to be most affected by the surface roughness. However, the surface roughness is not measured directly in the side-looking antenna footprint. The scatter plots in Figs. 9 and 10 and the calculated statistics show that a simplified one-layer model struggles to capture the dynamics of the sea ice T_B regardless of the simulation setup.

Figure 11 Histograms of the differences between the EMIARAD2 measurements and simulation setups for four antenna feeds.

Table 3 Performance metrics for the different T_B simulation setups in terms of coefficient of determination r², RMSE (K), and bias (K). For EMIRAD-2 channels four model setups are tested: Flat no snow, Flat snow, GO no snow, and GO snow.

4 Discussion and conclusions

In this paper we address the knowledge gap concerning the influence of the decimetre-scale surface roughness on the L-band brightness temperature of sea ice. We used the airborne laser scanner (ALS) data to characterise the sea ice surface and to produce the digital elevation model (DEM) of the sea ice surface. From the DEM we derived the probability distribution of surface slopes (α) and their azimuthal orientation (γ). We found that the probability distribution function of α (PDF_α) can be described with an exponential function regardless of the degree of roughness of sea ice surface. The exponent parameter (s_α) is a quasi-quadratic function of the standard deviation of surface heights. In the second part of this work, we used the PDF_α in a Monte Carlo simulation of the emission from a faceted sea ice surface. The effect of surface roughness is marginal near the nadir, accounting for up to 2.6 K decreases in T_B. The polarisation curves around Brewster's angle are most affected. The vertical polarisation decreases by 8 K, and horizontal polarisation increases by 4 K for the roughest ice, compared with the specular sea ice surface. The effect of the large-scale surface roughness on polarisation curves is not linear, with the degree of the surface roughness described by s_α, meaning that the alteration of the T_B curves is strongest for the roughest surface. The overall change of emission due to the large-scale surface roughness can be expressed as a superposition of change in intensity (H_α) and an increase in polarisation mixing (Q_α). The change in intensity depends primarily on the surface permittivity, whereas the polarisation mixing shows little dependence on ϵ. This parameterisation is suitable for all types of sea ice. However, the sensitivity analysis of the simplified emission model demonstrates that the expected change in T_B is comparable in magnitude to the uncertainty associated with the model input parameters.

The results have implications for the current and future L-band missions. The operational SMOS sea ice thickness product relies on near-nadir T_B observations (0–30^∘). Therefore, the large-scale surface roughness will have little effect on the retrieval. The SMAP and CIMR missions, which operate at incidence angles of 40 and 55^∘, respectively, are more exposed to the surface roughness effects. The effect on the vertical polarisation is stronger than on the horizontal polarisation. Lastly, we compared the simulation of the brightness temperature (with and without surface roughness) with the radiometer measurements. Unfortunately, this showed that our model does not capture the brightness temperature variability at the scale of 4.3 km. The inclusion of surface roughness is less important than the inclusion of a crude snow thickness parameterisation. This is confirmed by the sensitivity analysis of the model. Another possible explanation is that the sea ice in the studied region was highly heterogeneous in terms of its thickness and snow cover. Furthermore, a simple two-layer emission model used in this study has its limitations in capturing the T_B variability. Better results might be obtained if a multi-layer model together with the snow thickness measurements is used. With such a setup the direct inclusion of sea ice facets orientation in the radiometer field of view will be a valuable option to improve the T_B simulation. However, this would require in situ measurements of sea ice permittivity, snow thickness, temperature, and roughness as well as detailed characterisation of the antenna gain. Thus, the authors' recommendation for future studies is to measure the microphysical snow and sea ice properties together with surface roughness directly in the radiometer's field of view.

Appendix A: The following abbreviations are used in this paper

ALS	airborne laser scanner
CDF	cumulative distribution function
DEM	digital elevation model
GO	geometrical optics
ITS	inverse transform sampling
MILLAS	MIcrowave L-band LAyered Sea ice emission model
PDF	probability distribution function
RFI	radio frequency interference
SAR	synthetic aperture radar
SMAP	soil moisture active and passive mission
SMOS	soil moisture and ocean salinity mission
TB	brightness temperature