Submillimeter Wave Instrument radiometry of the Jovian icy moons Numerical simulation of the microwave thermal radiative transfer and Bayesian retrieval of the physical properties

Context. We address the issue of remote sensing of the surfaces of Galilean icy moons. We investigate the prospects for retrieval of the physical parameters of the surface of the Jovian icy moons from submillimeter wave radiometry data. Aims. We show that the model parameters could not be completely retrieved from the polarized radiometry data, but some of their combinations can be effectively constrained. Methods. The polarized radiative transfer in lossy porous ice was numerically simulated. A Bayesian maximum likelihood retrieval algorithm was developed and tested on the simulated data in a wide range of variation of the model parameters. The uncertainty of the retrievals was evaluated with the Cramer-Rao bounds. We established the combinations of model parameters that can be effectively constrained from the measured data. Results. We reveal that the effective scatterer size can be reliably constrained for a range of values where the scattering asymmetry parameter uniquely depends on the wave parameter, and for relatively high values of the single scattering albedo, for which the scattering in the medium is signiﬁcant. Similarly, the domains of reliable retrieval of the single scattering albedo and thermal skin depth are established.


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(Buratti , 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(Hagfors et al. , 1997Baron 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 twochannel 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.

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, where T 0 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: where ρ is the effective density of the subsurface, C is its specific heat, and κ is its effective thermal conductivity. The formal solution of Eq. (2), which satisfies the boundary condition (Eq. (1)) is the solar hour angle, LT is local time in hours, and is the diurnal thermal skin depth. The mean temperature T 0 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.

Microwave thermal radiation simulation
The regolith is 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.

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 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 mean free path length of the submillimeter radiation in the medium,x(Ω, Ω ) is the (Mueller) phase scattering matrix, fractured ice with bubbles and dust submillim eter wave radiomet er o r b it in g s p a c e c r a f t a n t e n n a p a t t e r n z θ Fig. 1. General sketch of the radiometric experiment.
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: where µ ≡ µ z . The kernel of the scattering integral iŝ Equation (6) is solved with the discrete ordinates (DO) technique 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 twocomponent 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 A24, page 3 of 13 A&A 644, A24 (2020) extinction coefficient, 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): 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): 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 where The thermal radiation source function f (τ) in (Eq. (8)) respectively reduces to in the equations for I and Q Stokes parameters, respectively.

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 temperature of 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) where where θ i and θ t refer to the incident and transmitted angles and are related by Snell's law where n is the refractive index of the medium relative to the free space. For the θ t > θ crit = sin −1 n −1 , total internal reflection occurs, so that T MA = 0, R MA = 1.

A matrix algebraic solution of the VRTE
To solve the polarized 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, whereÎ -is the vector of the unknown polarized radiances, B -matrix of the linear equation coefficients, the solution is expressed through the matrix exponent function The matrix exponent operator exp(−Bτ) can be diagonalized in the following manner: A24, page 4 of 13 Y. A. Ilyushin and P. Hartogh: Submillimeter Wave Instrument radiometry of the Jovian icy moons whereΓ 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 wheref 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 See above Eq. (25) where Diag(·) is the diagonal matrix, and Γ i are the corresponding diagonal elements of the matrixΓ, 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ŜÛ −1 , whereŜ is the diagonal matrix with the diagonal coefficients (Budak et al. 2004;Ilyushin 2011) 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.

The radiative properties of the medium
We simulate the 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 matricesx(Ω, Ω ) 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 matricesx(Ω, Ω ) have to be properly weighted and averaged. Black et al. (2001b) assume a power size distribution law, 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 where Q sca is the scattering efficiency (Bohren & Huffman 1998), and the total volume fraction v occupied by the scatterers, converges. Recalling that Q sca ∝ a 4 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.

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 allow us 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).

Maximum likelihood retrieval technique
The Bayesian approach suggests a formal, well-developed technique 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 A24, page 5 of 13 A&A 644, A24 (2020) 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 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 3 is a normalization constant to normalize the total a posteriori probability P(x|y) to unity. The maximum probability statex 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).

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 estimatorx(I) of any element x i of x is bounded as 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, is the inverse of the Fisher information matrix I F (x) (Helstrom 1968;Rodgers 2000) with elements 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 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 F (x) −1 . Thus, the minimum standard deviation that can possibly be achieved by any unbiased estimatorx i (y) can be directly evaluated. In our study the vector of observed 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 0 ), where x 0 is also a set of some values of the model 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.

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 τ 0 = 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 0 = 113 K and T 0 = 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 π/4 π/3 π/2 0 π−θ 18 π/6 π/4 π/3 π/2 0 π−θ 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)

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 π/4 π/3 π/2 0 π−θ π/4 π/3 π/2 0 π−θ π/6 π/4 π/3 π/2 0 π−θ π/4 π/3 π/2 0 π−θ and minimum, respectively. The polarization of the outgoing radiance is notable. The polarization degree Q/I is considerably higher 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. 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 ppolarized 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 −1 a = 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 feature of 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 range of 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 where x 11 is the upper left diagonal element of the phase scattering matrixx(Ω, Ω ), 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.

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 the single 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 dg/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.

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.