© Y. A. Ilyushin and P. Hartogh 2020

1 Introduction

The Galilean icy moons (Greenberg 2010) differ from most objects of the Solar System in several respects. In addition to remarkably high reflectance in the optical (Buratti 1985) and microwave spectrum (Black et al. 2001a,b), they exhibit anomalous polarization signatures of radar echoes (Black et al. 2001b; Campbell 2012), which are not typical for most planetary surfaces. For this reason, many studies address the issues of their structure, composition, and thermal regime.

Geological data imply high resurfacing rates, tidal cracking of the icy crust (Aglyamov et al. 2017), and modification of the surface by high-energy particles, that is, ionic and meteoritic bombardment. Due to the lock resonances of the icy moons, the surfaces are modified asymmetrically. Therefore, the optical properties of the leading and trailing hemispheres of the Galilean satellites differ significantly, and their radiative properties in other spectral bands should also be different. Disk-integrated and disk-resolved measurements in different wavelength ranges (Orton et al. 1996; Trumbo et al. 2018; de Kleer et al. 2019) provided information on the surface temperature of the moons and its variation over the diurnal cycle. Various optical models, based on single scattering and also accounting for the multiple scattering, have been used for the interpretation of the available optical data (Buratti 1985, 1995). The key radiative properties of the surfaces of the icy moons, such as the single scattering albedo and scattering phase functions, have been retrieved or constrained (Golombek & Buratti 1985; Verbiscer et al. 1990; Domingue et al. 1991; Buratti 1995). Radar echoes coming from the Galilean satellites, which also indicated high reflectivity and anomalous polarization signatures (Ostro & Shoemaker 1990; Black et al. 2001a), have also been interpreted and explained with different theoretical models.

As an explanation of the observed depolarization ratios (Ostro et al. 1980; Black et al. 2001a), double reflection from crater walls (Ostro & Pettengill 1978; Campbell 2012; Von Eshleman 1986), refraction scattering in the bulk ice (Hagfors et al. 1985, 1997; Baron et al. 2003) and coherent backscattering enhancement (Hapke 1990; Black et al. 2001b) have been discussed. Pronounced multiple scattering in the medium, which is necessary for effective coherent backscattering enhancement (Ilyushin 2012), implies deep penetration of radio waves into the scattering medium and a low volume absorption coefficient of the bulk ice.

Theoretical studies on the subsurface radar location at low frequencies (Bruzzone et al. 2011; Ilyushin 2014; Romero-Wolf et al. 2015; Hartogh & Ilyushin 2016; Aglyamov et al. 2017) are also quite optimistic about the ice penetrability for the electromagnetic waves and therefore its purity, which is required for low absorption.

Passive microwave radiometry has long been used in lunar studies to explore the structure of the lunar surface and its thermal regime. During the 1970s, the role of subsurface scattering of microwave radiation in the lunar regolith was revealed (Keihm 1982). Various models of subsurface scattering have been used for theoretical simulations, including statistical variations in the dielectric profile with an effective correlation length (Fisher & Staelin 1977), and Rayleigh (England 1975) and Mie (Keihm 1982; Tsang & Kong 1977) scattering theory. Effects of the nonsphericity of the scatterers have been investigated by (Greenberg et al. 1971; Purcell & Pennypacker 1973). The effects of vertical stratification of the regolith have also been extensively studied (Keihm & Cutts 1981).

The Submillimeter Wave Instrument (SWI) is a spectrometer and radiometer with two channels working in two frequency bands, namely 530–625 GHz and 1080–1275 GHz (Hartogh 2020). It will measure and record the radio brightness temperatures of the icy surface at these frequencies. These parameters and their temporal dependencies can be used to identify the physical properties of the icy surface, such as the thermal inertia, impurity content, and the size of the irregularities and small structures inside the crust. The JUICE spacecraft trajectory will include several fly-bys of the Galilean moons Callisto, Ganymede, and Europa before it finally orbits Ganymede, first in an elliptical orbit and then in two circular orbits at 5000 km and 500 km altitudes at the end of the mission.

The SWI does not have polarimetric capabilities. However, the plane of polarization of the receiver rotates when the spacecraft changes its orientation. Therefore, observing some particular location on the surface from different directions and orientations,the SWI will measure not only the total intensity of the microwave radio emission, but also two more of its Stokes parameters. This will provide the possibility to study the global structure of the icy crust and tectonic and resurfacing processes.

The previous instrument of this type, the Microwave Instrument on the Rosetta Orbiter (MIRO; Gulkis et al. 2007), which recently flew on the Rosetta space mission, was also a two-channel dual-frequency radiometer. Among others, it provided a large data set of surface temperatures of the 67/P Churyumov-Gerasimenko cometary nucleus. These data were used to retrieve key physical parameters of the upper layer of the cometary nucleus, providing meaningful insight into the thermal regime of the cometary crust and physical processes taking place in it (Choukroun et al. 2015; Schloerb et al. 2015). In the crust of the Jovian icy moons, volume scattering should be expected to be more significant in comparison to the lunar regolith and dusty crust of the cometary nuclei, because there is probably lower absorption.

For correct interpretation of the microwave radiometry observational data, relevant radiative transfer models that adequately treat transport of the electromagnetic radiation in millimeter and submillimeter microwave bands should be elaborated together with the approaches for retrieval of the physical quantities of interest. The development of such models and approaches is the primary motivation of this paper.

The secondary objective of this study is to investigate the thermal radiation polarization, and its role and importance for the retrievals. The role of polarization in the radiative transfer theory itself is now relatively well understood. The neglect of polarization in radiance calculations not only induces systematic errors (Adams & Kattawar 1970; Ilyushin 2012; Emde & Mayer 2018) but also prevents qualitative analysis of the physical mechanisms of the formation of the radiative field (Ilyushin 2019). However, the information value of the polarization in radiative transfer inverse problems is still questionable, and for any particular retrieval a special assessment should be made. This assessment is also one of the goals of this study.

The remainder of the paper is organized as follows. In Sects. 2 and 3, thermal and radiative transfer models of the investigated medium are presented. Formulation of the polarized radiative transfer equation, boundary conditions for the polarized radiance at the refractive boundary, the matrix algebraic technique of the radiative transfer equation solution, and derivation of the radiative properties of the medium are provided in the Sects. 3.1–3.4, respectively. Details of the Bayesian retrieval approach adopted in this study are provided in Sect. 4. Those who are sufficiently familiar with these subjects can skip these sections, simply using them as a reference, and can proceed directly to the results (Sect. 5). The results are presented separately for the monodisperse and polydisperse medium models (Sects. 5.1 and 5.2, respectively). Conclusions along with some final remarks are drawn in the final section. For convenience, a list of notations is provided at the end of the paper.

2 Thermal regime of the surface

The traditional approach to simulation of the thermal regime of atmosphereless celestial bodies under insolation assumes that all the incoming radiation is absorbed immediately on the surface. If optically thick particles constitute the bulk of grained regolith, solar radiation is absorbed by a layer with a thickness of the order of several particle diameters. If the particles are transparent, solar energy penetrates beneath the surface and is absorbed at depths which are no longer negligible. If the regolith particles are also opaque in the thermal IR band, heat must be transported to the surface by conduction. If the temperature T(t, z) at depth is higher than the surface temperature, this is the so-called solid-state greenhouse effect (Brown & Matson 1987; Urquhart & Jakosky 1996; Hapke 1996a,b).

However, estimated values of the effective penetration depth of the sunlight for the Galilean icy moons are rather small compared to the diurnal thermal skin depth, and are equal to zero in some model parameter fits (Urquhart & Jakosky 1996). For this reason, we neglect penetration of solar radiation beneath the surface. We also do not simulate thermal balance at the ice to the vacuum interface and set the temperature T(t, z) at the surface z = 0 immediately in the boundary condition, approximating it by its first Fourier series terms, $$T(t,0)=T_0+\Delta T\sin (\omega t),$$(1)

where T₀ and ΔT are the mean temperature and daily cycle amplitude, respectively, and ω is the angular frequency of the moon rotation.

The temperature in the subsurface layers of the icy regolith obeys the heat-conduction equation: $$\rho C\frac{\partial T(t,z)}{\partial t}=\kappa \frac{{\partial }^{2}T(t,z)}{\partial z^2},$$(2)

where ρ is the effective density of the subsurface, C is its specificheat, and κ is its effective thermal conductivity. The formal solution of Eq. (2), which satisfies the boundary condition (Eq. (1)) $$T(t,z)=T_0+\Delta T\sin \left({\varphi }_{\text{H}}-\frac{z}{l_{\text{T}}}\right)\exp \left(-\frac{z}{l_{\text{T}}}\right),$$(3)

where ${\varphi }_{\text{H}}=\frac{2\pi (\text{LT}-12)}{24}$ is the solar hour angle, LT is local time in hours, and $$l_{\text{T}}=\sqrt{\frac{2\kappa }{\rho C\omega }},$$(4)

is the diurnal thermal skin depth. The mean temperature T₀ and diurnal variation ΔT can be retrieved from the IR observations, for example, for Ganymede (Orton et al. 1996). Seasonal variations of the temperature are much weaker, because of the small eccentricity of the orbit of Jupiter, and are much slower because of its long rotation period. We therefore neglect these variations in our study.

3 Microwave thermal radiation simulation

The regolithis a random medium consisting of more or less densely packed individual particles, where the microwave radiation is emitted, scattered, and absorbed. Radiative transfer theory is typically applied for simulations of such radiation fields.

As the regoliths of the Jovian icy moons are thought to be loose and porous, radiative parameters of the medium (volume absorption and extinction coefficients, and scattering phase functions) can be evaluated from individual particle-scattering properties, as was done for example by Black et al. (2001b). If it is necessary to account for close electromagnetic interaction between the particles in a densely packed medium, various approaches can be applied, for example multiplying the differential scattering cross section by the static structure factor (Mishchenko 1992).

In this study we do not investigate the possibility of the immediate retrieval of physical parameters of the medium, as was done by Black et al. (2001b). Instead, we simulate the thermal emission of the media with different lengths of extinction of thermal radiation l_T, single scattering albedos Λ, and sizes of scatterers. The physical parameters of the medium, for example the dielectric permittivity of scattering particles and their size distribution, can be retrieved or constrained from the known radiative parameters of the medium. As SWI is sensitive to the polarization of the radiation, as noted above, a fully vectorial treatment of the microwave radiative transfer problem is performed.

3.1 Polarized radiative transfer equation

In this study we simulate the transfer of polarized microwave radiation in the homogeneous scattering and lossy medium. The only source of radiation is own thermal emission of the medium, while the temperature depends only on the depth beneath the surface. We also assume a perfectly flat medium boundary, that is, we neglect the surface roughness. To simulate the radiative transfer, we numerically solve the vectorial radiative transfer equation (VRTE) in a 1D vertically stratified medium $${\mu }_z\frac{\partial }{\partial \tau }I(\tau ,\Omega )=-I(\tau ,\Omega )+\frac{\Lambda }{4\pi }{\displaystyle {\oint }_{4\pi }\widehat{x}}(\Omega ,{\Omega }^{\prime })I(\tau ,{\Omega }^{\prime })\text{d}{\Omega }^{\prime }+f(\tau ,\Omega ),$$(5)

where I = {I, Q, U, V } is the vector of Stokes parameters, τ = z∕l is the optical depth, z is the vertical coordinate in the medium, l is the meanfree path length of the submillimeter radiation in the medium, $\widehat{x}(\Omega ,{\Omega }^{\prime })$ is the (Mueller) phase scattering matrix, Ω = {μ_x, μ_y, μ_z} is the unit direction vector, μ_x = sin θ cos ϕ, μ_y = sin θsin ϕ, μ_z = cos θ are the direction cosines, f(τ, Ω) is the source function, and θ and ϕ are the zenith and azimuth angle, respectively (see Fig. 1). In an electrically isotropic medium, the thermal source function is also isotropic, that is, f(τ, Ω) = f(τ). Further, in a flat layered medium with isotropic thermal sources of radiation, the angular polarized radiance distribution I (τ, Ω) is rotationally symmetric and does not depend on the azimuth angle ϕ. Thus the scattering operator reduces to a one-fold integral: $$\mu \frac{\partial }{\partial \tau }I(\tau ,\mu )=-I(\tau ,\mu )+\Lambda {\displaystyle {\int }_{-1}^1\widehat{x}}(\mu ,{\mu }^{\prime })I(\tau ,{\mu }^{\prime })\text{d}{\mu }^{\prime }+f(\tau ),$$(6)

where μ ≡ μ_z. The kernel of the scattering integral is $$\widehat{x}(\mu ,{\mu }^{\prime })=\frac{1}{4\pi }{\displaystyle {\int }_0^{2\pi }\widehat{x}}(\arccos \mu ,\varphi ,\arccos \mu ,{\varphi }^{\prime })\text{d}{\varphi }^{\prime }.$$(7)

Equation (6) is solved with the discrete ordinates (DO) technique $${\mu }_i\frac{\partial }{\partial \tau }{I}_i(\tau )=-I_i(\tau )+{\displaystyle \sum _{l,j}\widehat{x}}({\mu }_i,{\mu }_j)a_j{I}_j(\tau )+f_i(\tau ),$$(8)

where μ_i are directional cosines (discrete ordinates) according to the Gauss quadrature formula (Chandrasekhar 1950), a_j are the corresponding quadrature weights, and I_i(τ) is the polarized radiance in the discrete directions μ_i.

If the sources f(τ) are isotropic and nonpolarized, the third and fourth Stokes parameters of the radiance field (U and V) vanish due to the rotational symmetry of the VRTE solution, and so the resulting scattering matrix is effectively a 2 × 2 matrix and the VRTE solution is the angular distribution of the two-component Stokes vector with the I, Q components (Drusch & Crewell 2005).

In addition, the thermal source function of the elementary volume of the medium is ε(1 − Λ)B(T), where ε is the volume extinction coefficient, $$B(T)=\frac{2h{\nu }^3}{c^2}\frac{1}{{\text{e}}^{\frac{h\nu }{k_{\text{B}}T}}-1}$$(9)

is the Planck black body radiation function (Janssen 1993; Drusch & Crewell 2005), h and k_B are the Planck and Boltzmann constants, respectively, and ν is the frequency. If the quantum energy at the radiation frequency is small compared to the temperature ℏω ≪ kT, the Planck function is almost proportional to the temperature (the so-called Rayleigh-Jeans approximations): $$B(T)\approx \frac{2{\nu }^2{k}_{\text{B}}T}{c^2}=\frac{2k_{\text{B}}T}{{\lambda }^2},$$(10)

and so all the Stokes parameters can be expressed directly in the radio brightness temperature units (λ is the wavelength). At the working frequencies of the SWI instrument and temperatures of Ganymede, the Rayleigh-Jeans approximation (Eq. (10)) does not work very well, and the deviation from the exact expression (Eq. (9)) can reach tens of per cent. However, according to Janssen (1993) this approximation can still be used as a definition of equivalent radio brightness temperature T_b in order to express the intensity I and all the other Stokes parameters in temperature units (degrees of Kelvin): $$T_{\text{b}}\equiv \frac{{\lambda }^2}{2k_{\text{B}}}I.$$(11)

To evaluate the source function of the VRTE, we use the exact Planck function B(T) (Eq. (9)) of the physical temperature (Eq. (3)). Making use of the standard Taylor expansion, we linearize this expression with respect to ΔT and express it in radio brightness temperature units and finally obtain $$T_{\text{b}}({\varphi }_H,z)=T_{\text{b0}}+\Delta T_{\text{b}}\sin \left({\varphi }_{\text{H}}-\frac{z}{l_{\text{T}}}\right)\exp \left(-\frac{z}{l_{\text{T}}}\right),$$(12)

where $$T_{\text{b0}}=\frac{h\nu }{k_{\text{B}}}\frac{1}{{\text{e}}^{\frac{h\nu }{k_{\text{B}}{T}_0}}-1},$$(13) $$\Delta T_{\text{b}}=\frac{{\text{e}}^{\frac{h\nu }{k_{\text{B}}{T}_0}}}{{\text{e}}^{\frac{h\nu }{k_{\text{B}}{T}_0}}-1}\frac{h^2{\nu }^2}{k_B^2{T}_0^2}\Delta T.$$(14)

The thermal radiation source function f(τ) in (Eq. (8)) respectively reduces to $$f_I=(1-\Lambda )T_{\text{b}}({\varphi }_H,z),$$(15a) $$f_Q=0,$$(15b)

in the equations for I and Q Stokes parameters, respectively.

Fig. 1 General sketch of the radiometric experiment.

3.2 Boundary conditions

Solving the radiative transfer equation in the lossy semi-infinite medium with the temperature asymptotically approaching a constant value at large depth, we should require that the solution be bounded at infinity. To fulfil that, we have to exclude exponentially growing terms from the solution. In practice, we solve Eq. (5) in the layer of the medium with sufficient thickness. Due to extinction and absorption, outgoing radiance at the top boundary layer becomes almost insensitive to the boundary condition at the opposite boundary. At depths considerably greater than the mean free path length l and thermal skin depth l_T, the temperatureof the medium is almost constant and uniform, and the radiance distribution is isotropic. To implement this into a numerical scheme, one can simulate a reflective boundary surface, that is, a specularly reflecting horizontal plane.

The top boundary of the medium is a perfectly flat interface of the medium to free space. Boundary conditions for I and Q at this refractive boundary are expressed through the transmission and reflection matrices (Kattawar & Adams 1989) $$T_{\text{MA}}=\left(\text{}\text{}\begin{array}{cc}\alpha +\eta & \alpha -\eta \\ \alpha -\eta & \alpha +\eta \end{array}\text{}\text{}\right),$$(16a) $$R_{\text{MA}}=\left(\text{}\text{}\begin{array}{cc}{\alpha }^{\prime }+{\eta }^{\prime }& {\alpha }^{\prime }-{\eta }^{\prime }\\ {\alpha }^{\prime }-{\eta }^{\prime }& {\alpha }^{\prime }+{\eta }^{\prime }\end{array}\text{}\text{}\right),$$(16b)

where $$\alpha =\frac{1}{2}{\left(\frac{\tan ({\theta }_i-{\theta }_t)}{\tan ({\theta }_i+{\theta }_t)}\right)}^2,$$(17a) $$\eta =\frac{1}{2}{\left(\frac{\sin ({\theta }_i-{\theta }_t)}{\sin ({\theta }_i+{\theta }_t)}\right)}^2,$$(17b) $${\alpha }^{\prime }=\frac{1}{2}{\left(\frac{2\sin {\theta }_t\cos {\theta }_i}{\sin ({\theta }_i+{\theta }_t)\cos ({\theta }_i-{\theta }_t)}\right)}^2,$$(18a) $${\eta }^{\prime }=\frac{1}{2}{\left(\frac{2\sin {\theta }_t\cos {\theta }_i}{\sin ({\theta }_i+{\theta }_t)}\right)}^2,$$(18b)

where θ_i and θ_t refer to the incident and transmitted angles and are related by Snell’s law $$\sin {\theta }_i=n\sin {\theta }_t,$$(19)

where n is the refractive index of the medium relative to the free space. For the θ_t > θ_crit = sin⁻¹n⁻¹, total internal reflection occurs, so that T_MA = 0, R_MA = 1.

3.3 A matrix algebraic solution of the VRTE

To solve thepolarized radiative transfer equation in the scattering medium with uniform scattering properties, we follow the well developed matrix algebraic approach, which is outlined for example by Dave (1975) and Karp et al. (1980). In the matrix notation, $$\frac{\partial }{\partial \tau }\widehat{I}=\widehat{B}\widehat{I}+{\mu }_i^{-1}{\widehat{f}}_i,$$(20) $$\frac{\partial }{\partial \tau }\exp (-\widehat{B}\tau )\widehat{I}=\exp (-\widehat{B}\tau ){\mu }_i^{-1}{\widehat{f}}_i,$$(21)

where $\widehat{I}$ – is the vector of the unknown polarized radiances, $\widehat{B}$ – matrix of the linear equation coefficients, the solution is expressed through the matrix exponent function $$\exp (-\widehat{B}\tau )\widehat{I}(\tau )-\widehat{I}(0)={\displaystyle {\int }_0^{\tau }\exp }(-\widehat{B}\tau ){\mu }_i^{-1}{\widehat{f}}_i\text{d}\tau .$$(22)

The matrix exponent operator $\exp (-\widehat{B}\tau )$ can be diagonalized in the following manner: $$\exp (-\widehat{B}\tau )=Û\exp (-\widehat{\Gamma }\tau ){Û}^{-1},$$(23)

where $\widehat{\Gamma }$ is the diagonal matrix, and Û is an orthogonal operator. For the source function (Eq. (15)) with the vertical temperature profile (Eq. (3)) the right hand side of the Eq. (22) is represented in the form $${\displaystyle {\int }_0^{\tau }\exp }(-\widehat{B}\tau ){\mu }_i^{-1}{\widehat{f}}_i\text{d}\tau =Û{\displaystyle {\int }_0^{\tau }\exp }(-\widehat{\Gamma }\tau )\text{d}\tau {Û}^{-1}{\mu }_i^{-1}{\widehat{f}}_i,$$(24)

where ${\widehat{f}}_i$ is the right-hand part of Eq. (8). The temperature depth profile T(t, z) is a scalar factor permutative with the Û matrix, and so integration can be performed before the matrix operations. For the temperature profile (Eq. (3)) integration yields $$\begin{array}{l}{\displaystyle {\int }_0^{\tau }T}\exp (-\widehat{\Gamma }\tau )\text{d}\tau =T_{b0}\text{Diag}\frac{1-\exp (-{\Gamma }_i\tau )}{{\Gamma }_i}\\ +\Delta T_{\text{bDiag}}\frac{{\text{e}}^{-\left({\Gamma }_i+\frac{1}{{\lambda }_{\text{T}}}\right){\tau }_0}{\lambda }_T\left(\cos \left(\frac{{\tau }_0}{{\lambda }_T}-{\varphi }_{\text{H}}\right)+({\lambda }_T{\Gamma }_i+1)\sin \left(\frac{{\tau }_0}{{\lambda }_T}-{\varphi }_{\text{H}}\right)+{\text{e}}^{\left({\Gamma }_i+\frac{1}{{\lambda }_T}\right){\tau }_0}(-\cos ({\varphi }_{\text{H}})+{\lambda }_T{\Gamma }_i\sin ({\varphi }_{\text{H}})+\sin ({\varphi }_{\text{H}}))\right)}{{\lambda }_T{\Gamma }_i({\lambda }_T{\Gamma }_i+2)+2},\end{array}$$(25)

where Diag(⋅) is the diagonal matrix, and Γ_i are the corresponding diagonal elements of the matrix $\widehat{\Gamma }$, and λ_T ≡ l_T∕l is the normalized thermal skin depth. The system of linear algebraic Eq. (22) with the right hand side (Eq. (24)) is appended with the proper boundary conditions on the boundary surfaces of the slab, and scaled with the scaling transform operation (Karp et al. 1980). A perfectly conservative medium, which needs a special scaling procedure, is not expected to be the case in our problem. In addition, such a medium would not radiate at all, and so it would not be of interest to us. Therefore, it is enough to multiply both sides of Eq. (22) by the matrix $\widehat{S}{Û}^{-1}$, where Ŝ is the diagonal matrix with the diagonal coefficients (Budak et al. 2004; Ilyushin 2011) $$S_i=\{\begin{array}{cc}\exp ({\Gamma }_i{\tau }_0),& {\Gamma }_i<0;\\ 1,& {\Gamma }_i>0.\end{array}$$(26)

After the scaling transform, the system (Eq. (22)) is solved for the unknown outgoing radiative fluxes L_i escaping from the medium with the standard matrix algebra techniques. Necessary computations have been implemented with MatLab.

3.4 The radiative properties of the medium

We simulatethe medium as a matrix of lossy ice containing randomly distributed spherical voids (bubbles) of constant radius a. Thus, the diffraction parameter of the bubbles x = ka is one of the key parameters of the medium, which are to be retrieved from the observational data. The real permittivity of ice is approximately ε ≈ 3.15, and refractivity is therefore n ≈ 3.15^1∕2 ≈ 1.77. Losses in the medium are characterized by the single scattering albedo Λ. The scattering matrices $\widehat{x}(\Omega ,{\Omega }^{\prime })$ of the bubbles of different sizes have been calculated using the T-matrix algorithm (Moroz 2005) and then integrated over the azimuth ϕ numerically (Ilyushin 2019). In the T-matrix evaluation of individual scattering properties the loss in ice was disregarded, that is, the loss was only accounted for in the single scattering albedo Λ.

If the scatterer sizes are randomly distributed, the scattering matrices $\widehat{x}(\Omega ,{\Omega }^{\prime })$ have to be properly weighted and averaged. Black et al. (2001b) assume a power size distribution law, $$n(a)\text{d}a=K{a}^{-\beta }\text{d}a,$$(27)

in the whole range of scatterer size a from zero to infinity. Physical considerations require that the total scattering cross section of the unit volume $$\sigma ={\displaystyle {\int }_0^{\infty }\pi }{a}^2{Q}_{\text{sca}}(a)n(a)\text{d}a,$$(28)

where Q_sca is the scattering efficiency (Bohren & Huffman 1998), and the total volume fraction v occupied by the scatterers, $$v={\displaystyle {\int }_0^{\infty }\frac{4}{3}}\pi a^{3}n(a)\text{d}a,$$(29)

converges. Recalling that Q_sca ∝ a⁴ and Q_sca ≈ 2 for small and large a, respectively, convergence of the scattering cross section (Eq. (28)) implies that 3 < β < 7. The requirement of convergence of the total volume fraction integral (Eq. (29)) at the lower and upper limits results in contradictory demands β < 4 and β >4, respectively. It is worth mentioning that the β values retrieved by Black et al. (2001b) from the experimental data all lie in the range 3 < β < 4. Following Black et al. (2001b), we ignore the volume fraction divergence and investigate the whole range 3 < β < 7.

4 Bayesian retrieval approach

The regular almost spherical shape of the icy moons, contrary to the nucleus of the Churyumov-Gerasimenko (67/P) comet with its irregular shape (Choukroun et al. 2015), should allowus to collect a large amount of data for the thermal radiance of the surface for different observation angles and local solar times. If the available amount of the data is enough, physical parameters of the medium can be retrieved with one of the well-developed procedures, such as for example the Bayesian one (Rodgers 2000).

4.1 Maximum likelihood retrieval technique

The Bayesian approach suggests a formal, well-developed tech- nique for investigation of how observational information improves our a priori knowledge of the object. The result is a posteriori probability distribution, which relates to the objects or hypotheses. One can say that the Bayesian procedure of the initial a priori knowledge is a sample of an adaptive system with the ability to self-learn.The Bayesian theorem, or the Bayesian rule, is widely applied to a broad range of statistical problems. It provides large freedom of choice of a priori information from all of the available data by the proper selection of an a priori probability distribution. If the likelihood function accounts for all the observational data, a posteriori probability incorporates all the available information, both a priori known and obtained from the observation. In the classical approach, the a priori information is typically accounted for in the form of exact constraints on the unknown model parameters. These constraints, which are usually expressed as inequalities, are inconvenient for practical implementation and often do not yield satisfactory results. On the contrary, in the Bayesian approach our a priori knowledge of the object is expressed through a probability density. This approach often appears more convenient and flexible than the classical one, and various pieces of heterogeneous a priori information can be assimilated in a common uniform formal procedure.

For a nonlinear problem, such as the inversion of radiometry data, the maximum a posteriori approach is primarily used. A posteriori probability density for a nonlinear forward model with Gaussian measurement error and Gaussian a priori estimate is (Rodgers 2000) $$\begin{array}{ll}-2\text{ln}P(x|y)=\hfill & {[y-F(x)]}^T{S}_{\epsilon }^{-1}[y-F(x)]\hfill \\ \hfill & +{[x-x_a]}^T{S}_a^{-1}[x-x_a]+c_3,\hfill \end{array}$$(30)

where y is the vector of observed data values, x is the vector of the unknown model parameters, x_a is an a priori estimate of these parameters, F(x) is a forward model of the problem, S_a and S_ε are matrices of variance of the a priori estimate and measurement errors, respectively, and c₃ is a normalization constant to normalize the total a posteriori probability P(x |y) to unity.

The maximum probability state $\widehat{x}$ is the stationary point of the probability density function (Eq. (30)). To find it, we equate the derivative of Eq. (30) to zero. If the forward model F(x) is linear, the resulting equation is a system of linear algebraic equations and can be solved by standard linear algebra routines. For nonlinear forward models, special algorithms can be applied, such as for example Newtonian iterations (Rodgers 2000).

4.2 Bayesian estimate dispersion and Cramer-Rao bounds

To assess the quality of our Bayesian estimates of the medium parameters, we make use of the Cramer-Rao bounds (CRB). This is an inequality which sets a lower bound on the variance of the estimator. It has been demonstrated to be an informative and efficient tool to characterize minimum uncertainties of observational data inversion in various research fields (Gu et al. 2018; Holmes et al. 2013; McMahon et al. 2011). For a bibliographic and historical review of the subject, as well as a brief tutorial, the reader is referred to Gonsalves (1976), for example, and books and articles cited therein.

For the given physical model of the measurement process and statistical model of the noise, CRBs can be efficiently evaluated for any set of the unknown model parameters and any configuration of the measuring equipment. They directly provide the minimum estimation variance which can be possibly achieved by any unbiased estimator of any unknown model parameter. Here, CRBs are evaluated numerically using the microwave radiative transfer model described above, and a conventional Gaussian probabilistic model of the noise.

Let us denote x as the vector of unknown medium parameters, which are the diffractional parameter of the voids in ice x = ka, the single scattering albedo Λ, and the normalized thermal skin depth λ_T = l_T∕l. According to Helstrom (1968), the variance of any unbiased estimator $\widehat{x}(I)$ of any element x_i of x is bounded as $$\langle {({\widehat{x}}_i(y)-x_i)}^2\rangle \ge \text{CRB}{(x)}_{i,i},$$(31)

where CRB(x) is a square matrix whose dimensions are the same as the length of x. The CRB of x_i estimation corresponds to the ith diagonal element of CRB(x). Furthermore, $$\text{CRB}(x)=I_{\text{F}}^{-1}(x)$$(32)

is the inverse of the Fisher information matrix I_F(x) (Helstrom 1968; Rodgers 2000) with elements $$I_F{(x)}_{i,j}=\langle \frac{\partial \ln P(y|x)}{\partial x_i}\frac{\partial \ln P(y|x)}{\partial x_j}\rangle ,$$(33)

where P(y|x) is the likelihood of observing y given x. In the case of the multivariate Gaussian distribution of the observed quantities y with mean ⟨y⟩ and covariance matrix Γ, Eq. (33) can be simplified. If the covariance matrix Γ is independent of unknown parameters, (Eq. (33)) becomes $$I_{\text{F}}{(x)}_{i,j}=\frac{\partial \langle y\rangle }{\partial x_i}{\Gamma }^{-1}\frac{\partial \langle y\rangle }{\partial x_j},$$(34)

where the angular brackets denote the average. When using the known forward models described above, the evaluation of the Fisher information matrix I_F (x) only requires calculation of derivatives ∂⟨y⟩∕∂x_i, as well as knowledge of the instrumental errors Γ on the measurements. The derivatives from Eq. (34) can be calculated numerically from finite differences. The matrix I_F (x) is then inverted to obtain the CRB matrix CRB$(x)=I_{\text{F}}{(x)}^{-1}$. Thus, the minimum standard deviation that can possibly be achieved by any unbiased estimator ${\widehat{x}}_i(y)$ can be directly evaluated. In our study the vector ofobserved quantities y is the set of the Stokes parameters of thermal emission of the surface, registered at different observing angles and local times. The vector of unknown parameters ⟨x⟩ consists of three parameters of the model adopted for the simulations (the single scattering albedo Λ, the normalized thermal skin depth λ_T, and the size parameter of the Mie theory ka of size distribution exponent β). In our numerical study, the observed data vector is also simulated numerically with the same forward problem model y = F(x₀), where x₀ is also a set of some values of themodel parameters, which we call “true” values. To study the importance of polarization measurements, we include the second Stokes parameter Q in the vector y or exclude them from it.

5 Results and discussion

The polarized radiative transfer Eq. (8) was solved numerically with the matrix algebra technique as described above with the Gauss quadrature rule of sixteenth-order of accuracy. An optical depth of the medium of τ₀ = 30 was assumed for simulations, which is well above both the extinction and transport length of the medium, as well as the longest thermal wave length considered here.

According to the data provided by Orton et al. (1996), in our simulations we assume the values T_b0 = 100 K and ΔT_b = 20 K, respectively. Minimal and maximal temperatures correspond to local times LT = 6 h and LT = 18 h, which is sunrise and sunset. This mean radio brightness temperature T_b0 according to Eq. (13) corresponds to mean physical temperatures of T₀ = 113 K and T₀ = 126 K at frequencies of ν = 600 GHz and ν = 1200 GHz, respectively. The temperature deviation ΔT_b is much less sensitive to the approximations made and is about 20 K for both these frequencies. By simple calculation, it can easily be verified that these temperatures and their variations agree well with the estimated value of the thermal inertia of the Ganymede subsurface reported by Orton et al. (1996) ($\sqrt{\kappa \rho C}\approx 70{\text{Jm}}^{-2}{\text{K}}^{-2}{\text{s}}^{-1/2}$).

Fig. 2 Diurnal cycle of the angular distribution of the outgoing radiance, first Stokes parameter I (intensity). ka = 1, Λ = 0.1, λT = 1. Local times are shown by the numerical labels near each curve.

Fig. 3 Diurnal cycle of the angular distribution of the outgoing radiance, first Stokes parameter I (intensity). ka = 1, Λ = 0.9, λT = 1. Local times are shown by the numerical labels near each curve.

5.1 Monodisperse medium model (fixed size of the scatterers)

Numerically simulated angular dependencies of the outgoing intensity and polarization (I and Q Stokes parameters) for the monodisperse scattering medium model (single value of the scatterer size a) are shown in Figs. 2–7. Only three of 16 discrete ordinates correspond to directions of outgoing rays. All the others experience total internal reflection.

Discretized values of I and Q are plotted in the figures with dots. The curves are polynomial fits.

It is clearly seen that the angular distribution of the radiated intensity is sensitive to our model parameters, namely x = ka, Λ and λ_T = l_T∕l. In particular, the amplitude of the daily variation increases with the increase of normalized thermal skin depth λ_T. Maximal and minimal values of the radiated intensity roughly correspond to the surface temperature maximum and minimum, respectively. The polarization of the outgoing radiance is notable. The polarization degree Q∕I is considerablyhigher than typical values for thermal radiation of for example terrestrial atmospheric precipitation (Kutuza & Smirnov 1980; Kutuza et al. 1998; Evtushenko et al. 2002; Ilyushin & Kutuza 2016, 2017). This should be attributed to the polarizing properties of the refractive boundary which themselves are due to different transmission and reflection coefficients for s- and p-polarized radiation. Thus, for reliable measurement of both first and second Stokes parameters with a radiometer sensitive to linear polarization, one should properly account for mutual orientation of the instrument polarization plane and local normal to the boundary surface.

Bayesian likelihood functions (Eq. (30)) are calculated for the simulated sets of observed data for the whole 24 h rotation period (one observation per hour). Observational data vectors y (x) in our numerical simulations were formed from polarized radiances in the direction μ_i available from numerical solution for every hour of the 24h LT rotation period. For the sake of simplicity, observation angles correspond to the discrete ordinates of the used quadrature formula (Gauss quadrature of sixteenth order), properly corrected for the refraction at the boundary. A priori knowledge of x_a was assumed to be completely absent, that is, $S_a^{-1}=0$. Measurement error correlation matrix S_ε corresponds to the 1 K rms error of the radio brightness temperature measurement. Plots of logarithms of the Bayesian likelihood functions are presented in color graphs for two cases, when one (I) and two (I, Q) Stokes parameters are accounted for Figs. 8 and 9, respectively. In these figures, the Bayesian likelihood function is shown versus three parameters of our model: Λ, x = ka, and λ_T = l_T∕l. A global minimum of the deviation (Eq. (30)) is clearly seen, corresponding to a maximum of the likelihood function and uniquely determining the estimated values of three unknown parameters.

One can see that knowledge of the second Stokes parameter Q almost does not change the graph likelihood function, and therefore does not improve the quality of our retrieval. This means that a radiometer without polarimetric capabilities, like the SWI heterodyne instrument, can work well and can provide good retrievals. However, the instrument can be sensitive to the polarization state of radiation, that is, in can measure some linear superposition of the Stokes parameters, rather than pure intensity I. This featureof the instrument should be measured during the calibration tests and taken into account in the retrieval algorithm (Bayesian or any other). In other words, it does not make sense to measure the Stokes parameters separately and independently, but the sensitivity of the instrument to the polarized radiation should be accounted for. Here we considered the case of an idealized polarimetric radiometer measuring pure I and Q Stokes parameters. Nevertheless, due to the relative unimportance of the polarization state revealed by us, our assessment of the retrievals quality should be largely valid also for the case of the radiometer sensitive to the polarization state.

Cramer-Rao bounds determine the theoretical limits of accuracy potentially achievable for any unbiased estimator of any unknown parameters of the forward model F(x) (Eq. (31)). Assuming a radio brightness temperature measurement error of δT = 1 K, we evaluate the CRBs with the formulae (32) and (34) in a straightforward manner. The results of our evaluation of the Cramer-Rao bounds of estimates of the three model parameters x = ka, Λ, and λ_T are shown in Figs. 10–15.

As we can see, the CRB plots show that the model parameters (x, Λ, λ_T) can be retrieved or constrained from the SWI radiometry data but not all at the same time and not in the whole rangeof possible variation of their values. First, we can see that the CRB of the diffraction parameter x = ka is only low enough for satisfactory constraints for high albedos Λ and moderate x. This is intuitively expected, because to retrieve physical properties of the scatterers (air bubbles) they should be “clearly seen”, that is, the role of the scattering in the formation of outgoing radiance must be considerable. The range of x values where the CRB is low is roughly bounded by the value x = 4. This range coincides with the one where the scattering asymmetry parameter g with the increase of x monotonously grows from zero to its maximal value (see Fig. 16). The scattering asymmetry parameter g (Bohren & Huffman 1998) is the mean scattering cosine, which is the weighted average of the cosine of the scattering angle with the phase scattering function $$g\equiv \frac{{\displaystyle {\int }_0^{2\pi }{\displaystyle {\int }_{-1}^1{x}_{11}}}(\mu ,\varphi )\mu \text{d}\mu \text{d}\varphi }{{\displaystyle {\int }_0^{2\pi }{\displaystyle {\int }_{-1}^1{x}_{11}}}(\mu ,\varphi )\text{d}\mu \text{d}\varphi },$$(35)

where x₁₁ is the upper left diagonal element of the phase scattering matrix $\widehat{x}(\Omega ,{\Omega }^{\prime })$, which is the phase scattering function in the scalar radiative transfer approximation. This parameter largely characterizes the shape of the scattering pattern of the infinitesimal volume of the medium, in particular, it is a common measure of its elongation. The results of the CRB analysis shown in Figs. 12 and 15 expectedly indicate the fact that the Bayesian retrieval of the size parameter x works reasonably well in that range of variation of this parameter, where the shape of the phase scattering function changes significantly. In addition, it is worth noting the low values of CRBs of all three model parameters for all values of Λ and λ_T and low x.

Fig. 4 Diurnal cycle of the angular distribution of the outgoing radiance, first Stokes parameter I (intensity). ka = 1, Λ = 0.9, λT = 5. Local times are shown by the numerical labels near each curve.

Fig. 5 Diurnal cycle of the angular distribution of the outgoing radiance polarization. The 2nd Stokes parameter Q. ka = 1, Λ = 0.1, λT = 1. Local times are shown by the numerical labels near each curve.

Fig. 6 Diurnal cycle of the angular distribution of the outgoing radiance polarization. The 2nd Stokes parameter Q. ka = 1, Λ = 0.9, λT = 1. Local times are shown by the numerical labels near each curve.

Fig. 7 Diurnal cycle of the angular distribution of the outgoing radiance polarization. The 2nd Stokes parameter Q. ka = 1, Λ = 0.9, λT = 5. Local times are shown by the numerical labels near each curve.

Fig. 8 Likelihood function ln P(x|y) of the Bayesian procedure of the parameters retrieval of the unknown medium parameters x = ka, Λ, λT. The second Stokes parameter Q is ignored. The vector y of the observed quantities simulated for ka = 6, Λ = 0.6, λT = 3 model parameter values.

Fig. 9 Likelihood function ln P(x|y) of theBayesian procedure of parameter retrieval of the unknown medium parameters x = ka, Λ, and λT. The second Stokes parameter Q is included. The vector y of the observed quantities is simulated for ka = 6, Λ = 0.6, λT = 3 model parameter values.

Fig. 10 Cramer-Rao boundaries for the thermal skin depth λT vs. x = ka and Λ. Here, λT = 4. The second Stokes parameter Q is ignored.

Fig. 11 Cramer-Rao boundaries for the single scattering albedo Λ vs. x = ka and Λ. Here, λT = 4. The second Stokes parameter Q is ignored.

Fig. 12 Cramer-Rao boundaries for the wave parameter x = ka vs. x = ka and Λ. Here, λT = 4. The second Stokes parameter Q is ignored.

Fig. 13 Cramer-Rao boundaries for the thermal skin depth λT vs. x = ka and Λ. Here, λT = 4. The second Stokes parameter Q is included.

Fig. 14 Cramer-Rao boundaries for the single scattering albedo Λ vs. x = ka and Λ. Here, λT = 4. The second Stokes parameter Q is included.

Fig. 15 Cramer-Rao boundaries for the wave parameter x = ka vs. x = ka and Λ. Here, λT = 4. The second Stokes parameter Q is included.

Fig. 16 Scattering asymmetry parameter g of the bubbles in ice plotted against the diffraction parameter x = ka (monodisperse model, solid curve) and against the size distribution parameter β (polydisperse model, dashed curve).

5.2 Polydisperse medium model

Simulations for the polydisperse medium model show that some parameters of this model can also be constrained. As monodisperse model results showed that accounting for the second Stokes parameter Q does not significantly impact the Bayesian retrieval, for the polydisperse model we do not address this question and incorporate Q in all model retrievals, thus incorporating the second Stokes parameter Q measurements in all the simulations. Figure 17 shows the Bayesian likelihood function of simultaneous retrieval of the polydisperse medium parameters (size distribution parameter β, single scattering albedo Λ, thermal skin depth λ_T) for the true values β = 4.6, Λ =0.6, and λ_T = 3. The maximum of the likelihood function is clearly seen.

Figures 18–20 show the Cramer-Rao bounds for the polydisperse medium parameters β, Λ, and λ_T. It can be seen that β and Λ can be best constrained for high values of thesingle scattering albedo Λ. This is not unexpected, because these parameters characterize scattering properties of the medium, which are most pronounced when the albedo is high. For these values of the albedo, the CRBs for size distribution parameter β are reasonably low in the whole range considered here for this parameter (3 ≤ β ≤ 7). This is due to monotonous decrease of the scattering asymmetry parameter g with the growth of β throughout this range (see Fig. 16), contrary to the monodisperse model discussed above. One can reveal from Figs. 16 and 18 that the highest values of the derivative d g∕dβ correspond to the lowest values of the CRB bound. The thermal skin depth λ_T can always be constrained with an accuracy better than Δλ_T = 1, which qualitatively reveals that the thermal skin depth is on the order of the free path length of the radiation in the medium.

Thus our results demonstrate the principal possibility and basic limitations of retrieval for constraining the thermal radiative transfer model parameters from submillimeter wave radiometry with the SWI instrument.

Fig. 17 Likelihood function ln P(x|y) of the Bayesian procedure of the parameters retrieval of the unknown polydisperse medium parameters β, Λ, λT. The second Stokes parameter Q is included. The vector y of the observed quantities simulated for β = 4.6, Λ = 0.6, λT = 3 model parameter values.

Fig. 18 Cramer-Rao boundaries for the size distribution parameter β of the polydisperse model vs. β and Λ. Here, λT = 1. The second Stokes parameter Q is included.

Fig. 19 Cramer-Rao boundaries for the single scattering albedo Λ of the polydisperse model vs. β and Λ. Here, λT = 1. The second Stokes parameter Q is included.

Fig. 20 Cramer-Rao boundaries for the thermal skin depth λT of the polydisperse model vs. β and Λ. Here, λT = 1. The second Stokes parameter Q is included.

6 Conclusions and remarks

We developed a model of submillimeter wave propagation in the icy crust of the Galilean moons of Jupiter. Furthermore, we performed numerical simulations of the radiative transfer of polarized thermal radiation of the submillimeter spectral band in lossy ice with bubbles. We show that taking the polarization state of the thermal radiation into account does not significantly improve the retrievals. However, polarization sensitivity of the instrument has to be properly accounted for.

The possibility of Bayesian retrieval of the model parameters from radiometric data was also investigated. We show that some combinations of the model parameters (air-bubble size, single scattering albedo, thermal skin depth) can be effectively constrained from the microwave radiometry data.

Acknowledgements

Support from the Max Planck Institute of Solar System Research is greatly appreciated. The authors thank the reviewer and the Editor for valuable comments and constructive suggestions.

Appendix A List of notations

a – individual scattering radius in meters

$\widehat{B}$ – discretized VRTE matrix

C – specific heat of the ice

CRB – Cramer-Rao bounds

F(x) – forward model of the Bayesian retrieval

f – VRTE source function

${\widehat{f}}_i$ – vector of the VRTE source function values

g – scattering asymmetry parameter

h – Planck constant

I – Stokes parameters vector

I, Q, U, V – Stokes parameters

I_F – Fisher information matrix

$\widehat{I}$ – Stokes parameter vector

k – wave vector, m⁻¹

k_B – Boltzmann constant

K – normalizing constant of the scatterer size distribution

LT – local time in hours

l_T – thermal skin depth in meters

n – volume concentration of the scatterers

P – probability density

Q_sca – scattering efficiency

R_MA – reflection matrix of the polarized radiance at the refracting boundary

S_a – a priori estimate covariances matrix

S_ε – measurement errors matrix

Ŝ – scaling transform matrix

T_b – radio brightness temperature in Kelvin

T_b0 – mean radio brightness temperature in Kelvin

T_MA – transmission matrix of the polarized radiance at the refracting boundary

T – temperature in Kelvin

T₀ – mean temperature of the ice in Kelvin

Û – orthogonal operator diagonalizing the VRTE matrix

v – volume fraction in cubic meters

x = ka – Mie theory particle size parameter

$\widehat{x}(\cdot ,\cdot )$ – phase scattering matrix

x – vector of the model parameters

$\widehat{x}$ – maximum likelihood estimate of the model parameters

x_a – a priori known model parameters vector

y – vector of the observed data values

z – depth in meters

α, η, α′, η′ – auxiliary quantities

β – scatterer size distribution parameter

$\widehat{\Gamma }$ – diagonal matrix of the VRTE matrix eigenvalues

Γ_i – VRTE matrix eigenvalues

ΔT – daily variation of the temperature of the ice in Kelvin

ΔT_b – daily variation of the radio brightness temperature in Kelvin

${\widehat{\mu }}_i$ – z direction cosine vector

μ_x, μ_y, μ_z m μ – directional cosines

ρ – density

κ – thermal conductivity

ϕ – azimuth angle

ϕ_H – local solar hour angle

λ_T – normalized thermal skin depth

Λ – single scattering albedo

ε – volume extinction coefficient, m⁻¹

ε – dielectric permittivity

θ – zenith angle

θ_i – incident angle

θ_t – transmitted angle

τ – the optical thickness

τ₀ – total opticalthickness of the medium

σ – total scattering cross section in square metres

Ω – the unit vector of the direction

References

All Figures

	Fig. 1 General sketch of the radiometric experiment.
In the text

	Fig. 2 Diurnal cycle of the angular distribution of the outgoing radiance, first Stokes parameter I (intensity). ka = 1, Λ = 0.1, λT = 1. Local times are shown by the numerical labels near each curve.
In the text

	Fig. 3 Diurnal cycle of the angular distribution of the outgoing radiance, first Stokes parameter I (intensity). ka = 1, Λ = 0.9, λT = 1. Local times are shown by the numerical labels near each curve.
In the text

	Fig. 4 Diurnal cycle of the angular distribution of the outgoing radiance, first Stokes parameter I (intensity). ka = 1, Λ = 0.9, λT = 5. Local times are shown by the numerical labels near each curve.
In the text

	Fig. 5 Diurnal cycle of the angular distribution of the outgoing radiance polarization. The 2nd Stokes parameter Q. ka = 1, Λ = 0.1, λT = 1. Local times are shown by the numerical labels near each curve.
In the text

	Fig. 6 Diurnal cycle of the angular distribution of the outgoing radiance polarization. The 2nd Stokes parameter Q. ka = 1, Λ = 0.9, λT = 1. Local times are shown by the numerical labels near each curve.
In the text

	Fig. 7 Diurnal cycle of the angular distribution of the outgoing radiance polarization. The 2nd Stokes parameter Q. ka = 1, Λ = 0.9, λT = 5. Local times are shown by the numerical labels near each curve.
In the text

	Fig. 8 Likelihood function ln P(x|y) of the Bayesian procedure of the parameters retrieval of the unknown medium parameters x = ka, Λ, λT. The second Stokes parameter Q is ignored. The vector y of the observed quantities simulated for ka = 6, Λ = 0.6, λT = 3 model parameter values.
In the text

	Fig. 9 Likelihood function ln P(x|y) of theBayesian procedure of parameter retrieval of the unknown medium parameters x = ka, Λ, and λT. The second Stokes parameter Q is included. The vector y of the observed quantities is simulated for ka = 6, Λ = 0.6, λT = 3 model parameter values.
In the text

	Fig. 10 Cramer-Rao boundaries for the thermal skin depth λT vs. x = ka and Λ. Here, λT = 4. The second Stokes parameter Q is ignored.
In the text

	Fig. 11 Cramer-Rao boundaries for the single scattering albedo Λ vs. x = ka and Λ. Here, λT = 4. The second Stokes parameter Q is ignored.
In the text

	Fig. 12 Cramer-Rao boundaries for the wave parameter x = ka vs. x = ka and Λ. Here, λT = 4. The second Stokes parameter Q is ignored.
In the text

	Fig. 13 Cramer-Rao boundaries for the thermal skin depth λT vs. x = ka and Λ. Here, λT = 4. The second Stokes parameter Q is included.
In the text

	Fig. 14 Cramer-Rao boundaries for the single scattering albedo Λ vs. x = ka and Λ. Here, λT = 4. The second Stokes parameter Q is included.
In the text

	Fig. 15 Cramer-Rao boundaries for the wave parameter x = ka vs. x = ka and Λ. Here, λT = 4. The second Stokes parameter Q is included.
In the text

	Fig. 16 Scattering asymmetry parameter g of the bubbles in ice plotted against the diffraction parameter x = ka (monodisperse model, solid curve) and against the size distribution parameter β (polydisperse model, dashed curve).
In the text

	Fig. 17 Likelihood function ln P(x|y) of the Bayesian procedure of the parameters retrieval of the unknown polydisperse medium parameters β, Λ, λT. The second Stokes parameter Q is included. The vector y of the observed quantities simulated for β = 4.6, Λ = 0.6, λT = 3 model parameter values.
In the text

	Fig. 18 Cramer-Rao boundaries for the size distribution parameter β of the polydisperse model vs. β and Λ. Here, λT = 1. The second Stokes parameter Q is included.
In the text

	Fig. 19 Cramer-Rao boundaries for the single scattering albedo Λ of the polydisperse model vs. β and Λ. Here, λT = 1. The second Stokes parameter Q is included.
In the text

	Fig. 20 Cramer-Rao boundaries for the thermal skin depth λT of the polydisperse model vs. β and Λ. Here, λT = 1. The second Stokes parameter Q is included.
In the text