© ESO 2018

1 Introduction

Comets are considered to be remnants of the pristine material which formed during the early solar system, and as such, they may hold clues to their local formation environment and the evolution of the solar system. One way to assess the pristine nature of a comet is to measure the thermal properties of the surface layer because these control the energy transport through the surface and thus how much material is still unprocessed by solar insolation since its formation.

Energy transport through the surface of the nucleus is a very complex process (e.g. Gortsas et al. 2011). It is mainly initiated by the incoming solar insolation and the dust-ice matrix of the surface allows for the propagation of energy into the subsurface layers. The dust and ice can emit radiation and the volatile ice can also sublimate from the cometary material when heated. Since the volatiles are thought to be distributed throughout the surface and subsurface layers (De Sanctis et al. 2015), sublimation of ices could arise from different depths within the nucleus as heat is transported internally. Dust can then also be expelled from the surface as a result of ice sublimation. In addition, the inhomogeneous surface topography leads to shadowing and self heating (Keller et al. 2015).

Thermal inertia (hereafter TI) and surface roughness are two of the key mediators for the propagation of energy into a comet nucleus. The TI describes how efficiently the solar energy propagates into the surface and represents the ability of the surface to respond to the temperature forcing provided by the solar illumination. Roughness can locally enhance or reduce the amount of energy received and emitted at a given location. These two properties can be determined from the spatial and temporal surface temperature.

The Rosetta mission was launched in 2004 and arrived at its target destination, the comet 67P/Churyumov–Gerasimenko (hereafter 67P), in August 2014. The study of 67P represents an opportunity to better understand the physical properties of the nucleus as inferred from the energy transport in a cometary nucleus. For two years, the spacecraft observed the behaviour of the comet as it approached and receded from the sun. The mission ended in September 2016 when the spacecraft descended onto the surface. There were 11 instruments onboard Rosetta, and we used data obtained by two remote sensing units to determine the thermal properties of the surface and subsurface: MIRO (Gulkis et al. 2007, Microwave Instrument for the Rosetta Orbiter) and VIRTIS (Coradini et al. 2007, Visible and InfraRed Thermal Imaging Spectrometer). These are described in the next section.

Before the Rosetta mission, the Deep Impact spacecraft measured the thermal properties of comets 9P/Tempel 1 and 103P/Hartley 2 (Groussin et al. 2013). The TIs for these two comets were found to be lower than 250 JK⁻¹ m⁻² s^−0.5 for Hartley 2 and lower than 45 JK⁻¹ m⁻² s^−0.5 for Tempel 1 (3σ upper limit). The regions with exposed water ice had temperatures over 100 K higher than the sublimation temperature (~200 K), indicating the presence of ice-free dust at the subpixel scale. In addition, work by Davidsson et al. (2015) demonstrated that disc-integrated thermal emission from comets depends upon the type of roughness on the surface as well as the TI.

Regarding 67P, in a previous work by Schloerb et al. (2015), the TI in the Imhotep and Ash regions was calculated from MIRO data. These authors derived a low TI between 10–30 JK⁻¹ m⁻² s^−0.5 and found that the MIRO millimetre (hereafter mm) emission arises from a depth of approximately 4 cm, whilst the submillimetre (hereafter smm) emission comes from approximately 1 cm. Their work also shows residuals that are correlated to surface features, suggesting that regional heterogeneity may affect the surface properties of the comet. Choukroun et al. (2015) also calculated TI from MIRO data and obtained values of 10–40 and 20–60 JK⁻¹ m⁻² s^−0.5 for the submillimetre and millimetre channels, respectively. In addition, Gulkis et al. (2015) found low TIs of 10–50 JK⁻¹ m⁻² s^−0.5 from MIRO observations as Rosetta first approached 67P.

From the Philae lander, MUPUS (Spohn et al. 2007, Multipurpose Sensors for Surface and Sub-Surface Science) measured a TI of 85 ± 35 JK⁻¹ m⁻² s^−0.5 at the final landing site, Abydos (Spohn et al. 2015). This represents the best-fitting value from their thermal models to the infrared radiometer data from the MUPUS thermal mapper. This is higher than the MIRO measurements but this could be attributed to heterogeneities in the surface layer, such as porosity, temperature, dust-to-ice ratio, composition, dust distribution, and layer thickness (Spohn et al. 2015). The MUPUS measurement is also larger than unresolved observations made by the Spitzer infrared space observatory, which gives a low TI, 20 JK⁻¹ m⁻² s^−0.5 (Lamy et al. 2008), and is more consistent with the results obtained from MIRO measurements. The results from thermophysical modelling applied to light curves also gave a low TI value equal to 15 JK⁻¹ m⁻² s^−0.5 (Lowry et al. 2012).

For the asteroid Lutetia, which Rosetta passed during its flight to 67P, Keihm et al. (2012) used observations obtained by MIRO and VIRTIS to determine the TI and roughness. These authors found that a low TI (<30 JK⁻¹ m⁻² s^−0.5), which increased with depth, and a 50% fractional surface coverage with hemispherical craters was required to fit the model to the observations.

We aim to determine thermal properties of the surface layer of the nucleus of comet 67P, such as the TI and roughness, by analysing data from the VIRTIS and MIRO experiments. Both of these datasets provide measurements of the emitted and reflected radiance of the cometary nucleus. In the case of VIRTIS, the emitted radiance is related to the surface temperature (uppermost few tens of microns), whereas for MIRO, which operates at longer wavelengths, the emitted radiance is related to the subsurface temperature profile. We focussed our analyses on data obtained in September 2014. This choice was motivated by the opportunity to compare our final results with previous publications from the MIRO science team, such as by Schloerb et al. (2015), Gulkis et al. (2015), and Choukroun et al. (2015), who analysed MIRO data obtained during the same time period although from different regions.

In the next section, we describe the MIRO and VIRTIS instruments. Then, the first step towards estimating the TI involves assessing the datasets to find suitable data for comparison (Sect. 3.1). The thermal model applied to the data from both instruments and used to calculate the temperature gradients is described in Sect. 3.2, and the details of the MIRO radiative transfer model are given in Sect. 3.3. The VIRTIS modelling is described in Sect. 3.4. The results from MIRO and VIRTIS are given in Sects. 4.1 and 4.2, respectively.

2 Instruments

The MIRO instrument was a microwave spectrometer consisting of a 30 cm offset parabolic reflector telescope; a millimetre heterodyne receiver, operating at a centre-band frequency of 188 GHz (1.6 mm wavelength); and a submillimetre heterodyne receiver, operating at 562 GHz (0.5 mm wavelength; Gulkis et al. 2007). The frequency bands of the millimetre and submillimetre receivers had angular resolutions of 23.8 arcmin (6.9 mrad) and 7.5 arcmin (2.2 mrad), respectively, and each contained a single broadband channel for the measurement of near-surface temperatures. Additionally, a Chirp Transform Spectrometer (Hartogh & Hartmann 1990, CTS) was connected to the submillimetre receiver for the observation of absorption and emission lines of water, carbon monoxide, methanol, and ammonia. The work presented here focusses on the temperatures measured by the two broadband channels.

The instrument observed the source for 30 s and was calibrated in situ every 34 min by observing onboard hot and cold calibration targets. In September 2014, these typically had temperatures of 18 °C and –47 °C. The output from the receivers was proportional to the observed power and converted to an antenna temperature, T_A, by (1)

where P is the observed power density and k is the Boltzmann constant. In this work, we converted the measured antenna temperatures from MIRO to a brightness temperature, T_B, which is defined as the required temperature of a black body, which fills the antenna beam pattern to produce the observed power. A Planck function was used to calculate T_B. We assumed that each beam could be approximated as a Gaussian beam, and a multiplicative correction factor of was applied to the brightness temperatures from the submillimetre channel as a result of the difference in beam efficiency between the two channels, as described by Schloerb et al. (2015).

The VIRTIS instrument was a hyperspectral imager that consisted of a high-spectral-resolution point spectrometer (VIRTIS-H, 1.88–5.0 μm) and two mapping channels (VIRTIS-M-VIS, 0.22–1.05 μm; VIRTIS-M-IR, 0.95–5.1 μm; Coradini et al. 2007). We used data acquired by the IR mapping channel VIRTIS-M-IR (hereafter abbreviated as VIRTIS), the cryocooler of which failed in May 2015, thereby ending its science output. Its instantaneous field of view was 250 μrad × 250 μrad. At each basic acquisition, VIRTIS recorded a frame, a spectrally resolved spatial line (432 spectral bands times 256 spatial samples; 9.5 nm spectral sampling). A series of consecutive frames forms a data cube, a spectrally resolved two-dimensional image. For the cubes considered here, a basic acquisition typically had an exposure duration in the order of seconds, whereas the acquisition of a cube takes on the order of tens of minutes up to hours, a time span over which the observation and illumination geometry can change considerably from comet rotation and spacecraft motion.

The VIRTIS infrared spectrometer was mostly sensitive to the surface temperature (uppermost few tens of micrometre), rather than the centimetre subsurface as in the case of MIRO. The longest wavelengths that VIRTIS could measure generally provide the best diagnostics since the contribution from reflected flux is less important at these wavelengths. We always see the rising wing on the short wavelength side (the Wien’s part) of the thermal emission peak, starting from about 3.0–3.6 μm (Filacchione et al. 2016).

3 Methods

3.1 Observational overlap

In order to better determine the thermal properties of the nucleus in some regions, we decided to combine MIRO and VIRTIS data, providing more observational constraints than with a single instrument and enabling us to investigate possible heterogeneities between surface and subsurface thermal properties. We therefore identified when MIRO and VIRTIS were observing the same location on the comet at the same time. The data were averaged into one-minute time intervals and using a digital shape model of the nucleus, we could determine the footprints on the nucleus of the two instruments at a given time. We used a decimated version of shape model “cg-dlr_spg-shap7-v1.0”, which has 124938 facets (Preusker et al. 2017) with a horizontalresolution of approximately 20 m, and the SPICE software (Acton 1996). This way, we could determine the averaged viewing geometry of each observation from the positions of the spacecraft, the comet, and the sun, and then calculate the viewing angle (anglebetween facet normal vector and the vector connecting the facet and the spacecraft) and beam size of the submillimetre and millimetre beams, as well as the footprint of the observation and the facet illumination. Most of the time MIRO and VIRTIS observed different regions of the nucleus but occasionally they both looked at the same location.

The two instruments acquired data in slightly different ways. The MIRO millimetre and submillimetre beams probed a single circular area (assuming a Gaussian beam shape), giving two subsurface temperatures, one from each beam. At 30 km from the nucleus the millimetre and submillimetre beams had approximate diameters of 210 and 64 m, respectively. The millimetre beam hence observed between 55 and 367 facets of the decimated shape model for the observations given in Table 1, and between 3 and 26 facets in the submillimetre beam. The VIRTIS instrument, on the other hand, was a slit spectrometer, observing a narrow swath across the surface at a single exposure. At 30 km, it instantaneously observed an area of approximately 8 m × 2000 m and could resolve the temperature of the 256 facets placed along the swath. We therefore looked for observations in which the MIRO beam intersects the VIRTIS band. Since both instruments operated at different wavelengths and spatial resolutions, the roughness istreated differently in each case. The MIRO instrument takes large-scale topographical roughness at the scale of the shape model (metres) into account when inferring the TI, and VIRTIS models subfacet roughness, since it has a higher spatial resolution.

In September 2014, observations from the MIRO and VIRTIS instruments overlapped spatially and temporally ~150 times, and most of these occurred between September 1–2 and September 12–15. For this work, we chose to focus on five of the best times; these times were selected because they all have well-constrained MIRO and VIRTIS observations, good viewing geometries, and cover different areas of the nucleus. These times (specified in UT) are given in Table 1 along with the corresponding observed MIRO brightness temperatures. These observations provide data on five areas located on different geomorphological areas (El-Maarry 2015). On September 2 and 15, 2014, the instrument probed areas located in the Aten region, while on September 1 and 13, it observed areas located in the Seth and Ash regions, respectively. On September 12, MIRO and VIRTIS probed an area located at the border between theAsh and Aten regions. The measurements correspond to the pre-landing mapping phase shortly after arrival at the comet. During this epoch Rosetta was orbiting 67P at a distance of 30 km, at 3.38 AU from the sun.

3.2 Thermal model

To interpret the MIRO and VIRTIS data, we need to calculate the temperature of the surface and subsurface in response to insolation. The thermophysical temperature profile of a comet nucleus can be determined by solving the heat conduction equation with the assumption that the heat induced by the solar radiation is transported into depth and emitted into space (Groussin et al. 2013). This can also be used to sublimate volatiles as considered by Capria et al. (2017), but we neglected this part in this work, since water sublimation is negligible at 3.38 AU. The subsurface temperatures can be approximated with a one-dimensional thermal model using the heat equation (2)

where ρ is the bulk density in units of (kg m⁻³), c_p is the specific heat (J kg⁻¹ K⁻¹), and κ is thermal conductivity (W m⁻¹ K⁻¹). The partial derivatives and are for the changes of temperature with time, t, and depth, z, respectively. In Eq. (2), constant surface properties are assumed that do not depend on depth or temperature or time, and there is no internal heat source. The TI is defined as (3)

where I has SI units of (J K⁻¹ m⁻² s^−0.5). This is the quantity of interest, and with these equations, the temperature is found as only a function of TI, assuming Bond albedo and emissivity are known. We take the density to be 532 kg m⁻³ (Jorda et al. 2016) and assume the heat capacity to be 500 J kg⁻¹ K⁻¹, estimated from Robie et al. (1970).

For the surface boundary, the following condition is used to solve the heat equation (4)

where S, the solar flux at 1 AU, is equal to 1370 W m⁻², A_H is the Bond albedo, taken to be 0.0108, from a standard geometric value of 0.04 and phase integral value of 0.27 (Lamy et al. 2007), ϵ is the infrared emissivity, taken to be 0.95 (Birkebak 1972), σ is the Stefan– Boltzmann constant, r_h is the heliocentric distance, and i is the angle of incidence of the sun to the local surface.

As a first step, the illumination geometry and projected shadows are calculated for each facet of the shape model over one complete nucleus rotation. This is carried out for three dates for which we have overlapping VIRTIS and MIRO data: September 2, 13, and 15, 2014.

In the second step, we calculate the temperature of each facet with a thermal model that takes into account insolation, emission, self-heating, and heat conductivity with different values of TI (see Eqs. (2)–(4) and Groussin et al. 2013): 5, 10, 20, 40, 80, 160 and 320 JK⁻¹ m⁻² s^−0.5. In this work, we test the TI at these discrete values. In the conclusions, rather than discuss error bars, we consider a range of acceptable values denoted by the tested values. The thermal model is run over several nucleus rotations for each date (between 8 and 22 rotations depending on the thermal inertia) with a time step of 12.4 s and up to a depth of 10 diurnal skin depths (between 2.2 cm to 1.4 m depending on the TI). Running to these depths ensures convergence of the temperature profiles up to ten wavelengths of the MIRO millimetre channel (1.6 cm), and by checking this convergence, we negated the need for a full seasonal model for the uppermost layers of interest. The temperature at the lowest depth was assumed to be 30 K, and a depth is used such that . To account for the large number of facets (124938 facets) and the physical processes such as heat conductivity and self-heating, the code has been optimized and parallelized. On a computer with four CPU cores, it takes approximately 48 h to run a complete set of TI values at a given date.

As an output, we obtain, for a given date, the temperature profile for each facet of the shape model from the surface to a depth of 10 MIRO mm wavelengths or more. Figure 1 shows an example of temperature as a function of depth, for a TI of 80 JK⁻¹ m⁻² s^−0.5, on September 15, 2014 (UT 00:00). As expected for the dayside, the temperature decreases with depth. Figure 2 shows an example of surface temperature, for varying values of the TI in the range 5–320 JK⁻¹ m⁻² s^−0.5, on September 15, 2014 (UT 00:00). The maximum temperature decreases with increasing TI (since more energy penetrates into the nucleus by heat conductivity), as does the diurnal thermal amplitude.

3.3 Radiative transfer model for MIRO data

The MIRO radiometer measures antenna temperatures at millimetre and submillimetre wavelengths, described in Sect. 2, from which brightness temperatures for the observed areas are derived. Using the thermal outputs from Sect. 3.2 as inputs to a radiative transfer model, simulated brightness temperatures (hereafter SBT) are compared to measured brightness temperatures (hereafter MBT) for the different MIRO obervations. The TI of the material layer is constrained if the SBT curve for a given TI intersects the MBT curve. The SBT varies as a function of the TI and the relative complex permittivity (hereafter permittivity) of the material taken into account in the radiative transfer.

The permittivity, ϵ, consists of a real part, ϵ′, and an imaginary part, ϵ^″, which describe the polarisability of a material and the capacity of motion of free charges in the material when an electric fieldis applied, respectively. Their values are relative to that of a vacuum, which is equal to 8.85 × 10⁻¹² F m⁻¹. The real part of the permittivity is often referred to as the dielectric constant. The permittivity depends on the frequency, related to polarization mechanisms, composition, temperature, and bulk density of the material (Campbell & Ulrichs 1969; Mattei et al. 2014; Brouet et al. 2016).

For a given observed area and at a given wavelength of observation, the SBT is determined from the inverse Planck function taking into account the total spectral brightness intensity received from all facets covered by the MIRO beam. The spectral brightness intensity of a given facet varies as a function of the temperature profile of the subsurface and the permittivity of the material. It also depends on the emission angle, which refers to the orientation of the facet with respect to the main beam of the MIRO antenna. In this way, roughness at the scale of the shape model (of the order of metres) is treated in the interpretation of the MIRO measurements. A weighting function is applied to calculate the spectral brightness intensity of each facet. This function is based on an assumed Gaussian beam pattern for the MIRO instrument (Gulkis et al. 2007), where the centre of the beam contributes more strongly than the edges.

For one facet, the spectral brightness temperature intensity BT_ν as a function of the effective spectral brightness intensity BT_ν,eff and the surface emissivity, Υ, is as follows (5)

where ν is the frequency of the observation, θ_e, the facet emission angle defined as the angle between the facet normal and the main beam of MIRO antenna during the observation, and θ_t, the transmission angle in the subsurface layers of the facet. The transmission across a specular boundary is considered in order to compute the emissivity of the facet surface using the Fresnel reflectivity coefficients. The effective spectral brightness intensity is computed integrating the contribution of each subsurface layer, weighted by a radiative transfer function, which expresses the extinction of the spectral brightness intensity as it propagates towards the surface. The effective spectral brightness intensity computed for a given facet, assuming a scatter-free homogeneous layer under conditions of local thermodynamic equilibrium, is expressed as follows (Ulaby et al. 2014) (6)

where B(ν, T(z)), the Planck function for a given frequency of observation and temperature, T(z), which varies with depth and is derived from the thermal model, and δ_el, the electrical skin depth of the layer at a given wavelength. The value H is the thickness of a layer equalling ten thermal skin depths. The electrical skin depth is inversely proportional to the imaginary part of the permittivity and proportional to the wavelength and the square root of the dielectric constant. It characterizes how deep an electromagnetic field can penetrate into the material.

Regarding the MIRO observations, we expect the submillimetre receiver to be sensitive to the thermal radiation integrated over the topsubsurface of the layer and the millimetre receiver to be sensitive to the thermal radiation integrated over a greater depth. The electrical skin depth is constrained by a permittivity model based on experimental results obtained with lunar regolith samples reported in Olhoeft & Strangway (1975). The temperature dependence of the dielectric constant is considered to be linear (Calla & Rathore 2012; Brouet et al. 2015).

Fig. 1 Temperature as a function of depth (on the surface and 1 cm, 2 cm, and 4 cm below the surface) for a TI of 80 JK−1 m−2 s−0.5, for September 15, 2014 (UT 00:00).

Fig. 2 Temperature on the surface for various values of the thermal inertia in the range 5–320 JK−1 m−2 s−0.5, for September 15, 2014 (UT 00:00).

Fig. 3 Radiance measurements from VIRTIS projected onto the shape model on 12 Sep. 2014 (left) and 15 Sep. 2014 (right).

3.4 Radiance model for VIRTIS data

With the VIRTIS data, we focus our study on only two specific dates in Table 1: September 12 and September 15, 2014.

Figure 3 shows the shape model of 67P rendered from the observing view of the spacecraft at those two epochs, and the radiance measured at 4.95 μm by VIRTIS around that time. This wavelength was chosen because it is the longest in the VIRTIS range that is not directly at the detector edge and is outside certain spectral features that may be real or else caused by calibration issues (Filacchione et al. 2016). The measured radiance is a combination of directional thermal emission and reflected solar radiation. The observed areas cover the regions of Aten, Babi, Khepry, and Imhotep, as defined by El-Maarry (2015). They are located in the northern hemisphere, at latitudes between 0 and +60 degrees and longitudes between 60 and 150 degrees. Because the data set covers several morphological regions, we expect significant variations in roughness values. Facets shared by the two datasets cover mainly the northern part of Imhotep and a small fraction of Aten.

In order to compare data and model, we first convert the modelled surface temperature into radiance, I(λ, T), by combining reflected solar spectral radiance and emitted IR radiation as follows (7)

where α and ϵ are the constant geometric albedo and emissivity (respectively 0.05 and 0.95), F the solar flux modulated by a reflectance function (the Lommel-Seeliger Law), and B the black body radiation at wavelength λ and temperature T. We neglected the directionality of epsilon, which is affected by spatially unresolved surface roughness.

For a given facet on the shape model, the solar spectral radiance is given by Planck’s law, heliocentric distance, and a photometric model. We chose the Lommel–Seeliger (Hapke 2012) correction, which depends only on incidence and emission angles. This equation is written as (8)

where h, c, k are Planck’s constant, the speed of light, and the Boltzmann constant, respectively. The value R_Sun is the sun radius, r_h the heliocentric distance. The values μ₀ and μ are the cosines of incidence and emission angle.

Finally, the emitted IR radiation is the Planck function at T and λ, is (9)

We applied our model to all facets covered by the VIRTIS data set, except those in shadows for which a more complex photometricmodel should be considered, accounting for multiple scattering. Shadows were calculated with a custom ray tracing engine¹ widely used to simulate Rosetta images, and calibrated against observations.

3.5 Importance of roughness

In addition to TI and emissivity, it has been noted by many authors that roughness at spatial scales unresolved by the instrument plays an important role when determining the surface temperature of an airless body with a low TI, at many length scales (Kuehrt et al. 1992; Keihm et al. 2012; Groussin et al. 2013; Davidsson et al. 2015). First, the micro-topography (on centimetre scale) of the surface leads to re-radiation and multiple scattering of the incoming solar and thermal radiation, which can locally enhance or decrease the amount of energy received by forming heat traps or sinks. The temperature therefore varies on subpixel scales and we consequently observe a thermal radiation that cannot simply be represented by a single Planck function. One example is the case of craters, in which the shadowed wall can be heated up by radiation from the illuminated wall. Second, the illuminated fraction of the surface seen by the detector is influenced by the roughness, and so the energy received is a function of the relative position of the observer and the energy source. Both issues limit the accuracy of temperature retrieval by remote sensing and must be accounted for. It should be noted that simulating accurately the roughness at the subpixel resolution of the OSIRIS (Keller et al. 2007, Optical, Spectroscopic, and Infrared Remote Imaging System) camera over a large surface is computationally far too expensive to be achieved in most cases. In practice, the surface temperature is first approximated with a classical Lambert surface. In a second step, the roughness is modelled by multiplying the Lambert temperature with a correction factor that is a function of wavelength, incidence angle, emission angle, projected phase angle, and a statistical description of the subpixel topography.

The subpixel topography can be represented in many ways, including fractal surface, coherent noise, and craters. Davidsson et al. (2015) showed that, for non-directional terrains (i.e. with facets orientated in all directions), it does not matter much in the final result which description is used. We adopted the approach proposed by Kuehrt et al. (1992), representing topographic variations with a set of spherically shaped depressions. The model is scale independent, and therefore our only free parameters are the area fraction covered by the craters and opening angle of the spherical segments. The opening angle refers to the centre of the corresponding sphere so an angle of 180° creates a half-spherical depression.

In this study, we want to assess whether the roughness of surface of 67P can be estimated reliably from the combination of infrared measurements by the Rosetta VIRTIS instrument and a thermal model of the surface, described in Sects. 2and 3.2, respectively. To do so, we assume that after accounting for the TI, any further discrepancy between model and data can be corrected by a factor that results solely from roughness effects. In this approach, we first derive a best TI for the observed area and then allow roughness to vary spatially in our model. We know all the geometric parameters for any facet of the shape model at any time, so we can use the roughness model to relate this correction factor to the area covered by the spherical segment depressions. This roughness estimation is then compared to OSIRIS images of the same region, as they have typically much higher spatial resolution than VIRTIS data. We finally assess how the spatial distribution of our retrieved thermal roughness (micrometre–millimetre scale) compares to the visual roughness (centimetre scale) in the images.

We therefore have three unknowns: TI and two roughness parameters (crater density and opening angle). Ideally, combining three inputswould provide good constraints on the parameters: for instance, the radiance in at least three wavelength channels for a given facet, or radiance measurements at the same wavelength but acquired with different geometries. However, since this is not always available in our dataset, we choose a less constrained approach, as outlined above.

4 Results

4.1 MIRO results

Figures 4–8 show the SBT as functions of the electrical skin depth for each of the thermal inertia values used in the thermal model, in the range 5–320 JK⁻¹ m⁻² s^−0.5, for the five dates given in Table 2. As stated earlier, we constrain the TI when SBT intersects MBT for electrical skin depths less than ~3.0 and ~1.0 m for the millimetre and submillimetre cases, respectively.

For example, in Fig. 4 in the millimetre channel, the SBT with TI = 5 JK⁻¹ m⁻² s^−0.5 intersects the MBT but TI = 10 JK⁻¹ m⁻² s^−0.5 does not. We hence derive that the TI in the millimetre channel for September 1 is <10 JK⁻¹ m⁻² s^−0.5, as given in Table 2.

We choose not to further constrain the TI value by interpolating between the curves. Linear interpolation seems to be insufficient, so instead of devising a complex interpolation scheme, we quote the upper and lower bounds (where possible) from the discrete TI values, which initially went into the thermal model.

The MBT of the different observations performed with the millimetre receiver are consistent with a TI less than 80 JK⁻¹ m⁻² s^−0.5, while the value for September 1, 2014 implies that the TI should be less than 10 JK⁻¹ m⁻² s^−0.5. The results from the submillimetre receiver on September 1 and 2, 2014 are also consistent with a TI less than 80 JK⁻¹ m⁻² s^−0.5, but the other dates imply higher values of 160 and 320 JK⁻¹ m⁻² s^−0.5 are also possible for their observations. The results are compiled in the Table 2.

4.2 VIRTIS results

Disentangling TI from roughness effects is a complex problem for which we do not always have enough constraints. To reiterate, we attempt to fit first the TI as it relates directly to the physical properties of the material. We then add the surface roughness as a second parameter effectively describing which subpixel area contributes to the flux that reaches the detector.

Therefore, the first step in our inversion of roughness parameters is to determine which thermal inertia brings our model the closest to the measured radiance. We calculated the radiance for the seven values of TI investigated by the thermal model described above (TI = 5, 10, 20, 40, 80, 160, 320 JK⁻¹ m⁻² s^−0.5). We found that the modelled flux is almost always below the measured flux (see Fig. 9). As we increase TI, the average flux decreases. This is because higher TI means that facets must be illuminated for a longer time before reaching their maximum temperature. Our simulations show that a good fit is achieved for a TI larger than 80 JK⁻¹ m⁻² s^−0.5. We therefore use this TI value as a lower limit in our analysis. Higher values lead to an equally good fit, suggesting that TI is not the main reason for the dispersion we observe in Fig. 9, although the fit is slightly worse for TI values of 160 and 320 JK⁻¹ m⁻² s^−0.5 in the case of September 15, 4.75 μm.

Having selected TI = 80 JK⁻¹ m⁻² s^−0.5 as our reference TI, we then postulate that the remaining discrepancy between model and observation is caused by the surface roughness. Following the approach of Kuehrt et al. (1992), we define the corrected flux for a rough terrain as I_r = c_r × I, with c_r a function of incidence and emission angles defined by the shape model facet normals, projected phase, crater density, crater opening angle, and wavelength. The projected phase, or azimuth, is the angle between the observer direction and solar direction vectors projected onto the facet surface. In this approach, the only free parameters are the crater density and the crater opening angles. Our test cases show that the influence of opening angle is less important than that of the crater density for the wavelengths considered. Therefore, we fix the opening angle to a value of 180 (which leads to maximal c_r) and find the crater density that minimizes the quantity I_r − I_m: the difference between the calculated radiance, I_r, and the VIRTIS measured radiance, I_m. The crater density is allowed to vary between 0 and , the maximum density for packing equal circles on a two-dimensional plane. This approach effectively gives us a lower bound for the crater density on each facet in the region of interest.

We note that in principle, one should find the best parameters by varying both TI and roughness in parallel when getting the minimum least-squares fit. However, this is a large parameter space and running the complete inversion goes beyond the scope of this work.

Figure 10 shows an example of the corrected radiances for two different dates. The most remarkable effect is the strongly reduced scattering of our data points with respect to Fig. 9, indicating a very good agreement with the observations. The residual discrepancies can be explained by the fact that we have used a relatively coarse shape model (124938 facets, 20 m resolution) and the same value of TI for all facets. A finer spatial resolution for both parameters would certainly improve the fits, as would a more advanced roughness model and spatially varying emissivity. In addition, it should be noted that VIRTIS measurements are not instantaneous, but require up to 20 min in this work, during which the viewing geometry changes as the comet rotates and the spacecraft moves. This also affects our inversion.

Figures 11 and 12 illustrate our results; the retrieved crater density is projected on the shape model. One can see that the distribution of densities is not random, but rather organized in regions whose distribution resembles the morphological regions derived from optical data alone. We find a similar crater density for regions with a similar morphology, although this retrieved parameter may appear sometimes counter-intuitive. For instance, when looking at morphology alone, the Imhotep plain appears to be separated in two areas with different textures. The northern part is rougher with many metre-sized boulders, while the southern area is very smooth and has topographic variations that are basically undetectable at the best OSIRIS resolution (<50 cm px⁻¹). Interestingly, our roughness inversion also discriminates between these two subregions, thus helping to validate our method, but the area with highest crater density is the subregion that appeared to be the smoothest in OSIRIS images.

Although this may seem surprising at first, it is actually explained by the fact that the roughness that plays a role for IR measurements is simply at a much lower scale than topographic variations. Taking this into account, it seems reasonable to postulate that consolidated terrain is likely to show uniform textures on a centimetre scale, while a visually smooth terrain may be covered in pebbles and dust aggregates that form a very rough surface at small scales.

This result is consistent with several observations from the Rosetta lander Philae. Images acquired by the ROLIS camera during the descent to the Agilkia region on November 12, 2014 show that a region that appeared as a smooth dusty plain in OSIRIS images is actually covered by centimetre-sized pebbles when observed at higher resolution (Mottola et al. 2015). This region has avery granular appearance as seen by ROLIS (Rosetta Lander Imaging System) and may be a result of physical processes occurring on the comet such as airfall of small material or talus deposits. The same may be true for the Imhotep dust plain, as revealed by our thermal analysis. On the other hand, CIVA (Comet Infrared and Visible Analyser) images obtained at Abydos after landing show that the consolidated terrain, although very rough on a metre-sizescale, display flat and smooth facets at centimetre to decimetre scale (Bibring et al. 2015). Flat and smooth morphologies could be a result of a variety of process, including the disaggregation of boulders by complete erosion to create rounded structures, or the fragmentation of hard consolidated material to reveal flat surfaces(Mottola et al. 2015). Further study is required to link the measured roughness to these morphological processes.

Fig. 4 Simulated brightness temperatures (SBT; coloured lines), as functions of the electrical skin depth, for various values of TI in the range 5−320 JK−1 m−2 s−0.5 (SI units). The black dotted lines represent the MBT from the millimetre and submillimetre receivers for September 1, 2014 at 23:48. The millimetre SBT curves with TI values less than 10 JK−1 m−2 s−0.5 intersect the millimetre MBT line; and the submillimetre SBT curves with TI values less than 80 JK−1 m−2 s−0.5 intersect thesubmillimetre MBT line.

Fig. 5 Same as Fig. 4 but for September 2, 2014. The results are given in Table 2.

Fig. 6 Same as Fig. 4 but for September 12, 2014. The results are given in Table 2.

Fig. 7 Same as Fig. 4 but for September 13, 2014. The results are given in Table 2.

Fig. 8 Same as Fig. 4 but for September 15, 2014. The results are given in Table 2.

5 Discussion and conclusions

5.1 Thermal inertia

All of the results obtained from the MIRO millimetre observations imply that the TI has a low value, <80 JK⁻¹ m⁻² s^−0.5. The result from September 1 even implies that it should be <10 JK⁻¹ m⁻² s^−0.5 in the Seth region. For this date, the submillimetre channel gives a higher upper bound of 80 JK⁻¹ m⁻² s^−0.5, similar to other results.

Constraining the TI from the subsurface temperature observations is sometimes difficult, however. For the measurements on September 12, 2014 and September 15, 2014, the submillimetre observation cannot constrain the TIso well (5–320 JK⁻¹ m⁻² s^−0.5 and >20 JK⁻¹ m⁻² s^−0.5, respectively) based on the assumed dielectric model. The results from the millimetre observations on these two dates are better constrained however, giving an upper bound of 80 JK⁻¹ m⁻² s^−0.5.

In comparison, VIRTIS work implies a best-fitting value of 80 JK⁻¹ m⁻² s^−0.5 across the observed Aten, Babi, Khepry, and Imhotep regions. The VIRTIS value was calculated over a much larger area than the MIRO values as the viewing area of VIRTIS is larger than the MIRO beam area. However, this does not seem to make a difference. With the VIRTIS data restricted to a subset containing only the facets viewed by the MIRO beam, the VIRTIS result remains the same; a similar thermal inertia value is predicted.

In Fig. 13, we compare our measurements for the TI from this work with previous observations. Our measurements for the TI from MIRO are consistent with past MIRO measurements of the TI (Gulkis et al. 2015; Schloerb et al. 2015; Choukroun et al. 2015) and there is also some overlap with the lander observation (Spohn et al. 2015). The measured VIRTIS TI is in good agreement with the MUPUS value, which was also determined from an infrared instrument: MUPUS-TM operates at 5–25 μm and VIRTIS observes at 2–5 μm. In addition, further VIRTIS calculations have yielded a TI in the range 25–170 JK⁻¹ m⁻² s^−0.5 (Leyrat et al. 2016). Neglecting spatial and depth variations, a bulk TI value in the range 30–50 JK⁻¹ m⁻² s^−0.5 would fit all of the measurements from Rosetta instruments, except for the MIRO millimetre measurement made here on September 1. This value is similar to the early TI estimates from Lamy et al. (2008) and Lowry et al. (2012), which are much lower than most of the Rosetta instrument measurements.

Taken together, the Rosetta measurements are consistent with a low TI (≤170 JK⁻¹ m⁻² s^−0.5), and there is notable overlap in the error bars of these results.

Fig. 9 Left panel: example of modelled radiance minus observed radiance at 4.95 μm on September 15, 2014 for TI = 5 JK−1 m−2 s−0.5. Increasing the TI makes the plot more symmetrical around zero. Right panel: mean value of model minus observation as a function of TI for two epochs and two wavelengths. We observe that increasing the TI beyond 80 JK−1 m−2 s−0.5 does improve the fit, and slightly worsens it for the case September 15, 4.75 μm. We therefore use this TI value as a lower limit in our analysis.

Fig. 10 Difference between modelled rough radiances and VIRTIS measurements for September 12, 2014 (left panel) and September 15, 2014 (right panel). Although we still observe some scattering, the dispersion is far smaller than when not considering roughness, as in Fig. 9.

5.2 Roughness

In the MIRO analysis, roughness is considered on the scale of the horizontal shape model resolution because the angle between each facet normal of the shape model and the direction to the spacecraft is taken into account when computing the SBT. We considered adding subfacet roughness into the MIRO analysis but this would be computationally intensive. In addition, the uncertainty on the TI is larger than the expected uncertainty caused by subfacet roughness because of the valid range of values considered for the electrical skin depth.

The VIRTIS analysis showed that the roughness could be considered on subpixel scales that are below the resolution of the OSIRIS images. While work remains to be done to define the surface roughness properly, our approach shows that one can use IR radiance measurements and a simplified roughness model to separate the surface into regions of different properties. This study raises the prospect of consistently mapping the full surface at various epochs, in particular in areas like Imhotep in which significant changes have been detected at large and small scales over the course of the Rosetta mission (Groussin et al. 2015). We conclude that the roughness below the resolution of the images plays a key role in interpreting thermal infrared observations, and that one cannot use the visible images to estimate the surface roughness at this scale, since the small-scale (IR wavelength) roughness is not well correlated with the roughness on the large scale (visible wavelength).

This work demonstrates the importance of comparative studies of results from different instruments. By combining observationsfrom the MIRO, VIRTIS, and OSIRIS instruments on Rosetta, we have gained new insights into the small-scale roughness of the surface and found more evidence for a low TI of the comet nucleus.

Fig. 11 Left panel: regional map of the comet from El-Maarry (2015) with morphological units. Right panel: roughness parameter, i.e. the density of the half-spherical depressions projected onto the shape model. The density value is represented by the colour bar on the right.

Fig. 12 Half-spherical depression density overlaid on OSIRIS images, 3D projected onto the shape model with the same colour scheme as in Fig.11. The images are centred on approximately 130° longitude and 10° latitude. One can see that areas with similar morphologies share similar roughness. Visually smooth terrains (at the scale of OSIRIS data, i.e. 75 cm px−1) can display high thermal roughness, indicating that the spatial scale of this parameter is far below the resolution of our images.

Fig. 13 Summary of measured TIs (units: J K−1 m−2 s−0.5) from a variety of sources. Measurements in black represent previous published works and coloured bars represent this work. From Left to right, these measurements represent the TI from: (a) Spitzer observations (Lamy et al. 2008); (b) a thermophysical model applied to light curve observations (Lowry et al. 2012); (c) MIRO observations by Gulkis et al. (2015); (d) MIRO observations by Schloerb et al. (2015); (e) MIRO observations by Choukroun et al. (2015), with TI from the millimetre channel in black and TI from the submillimetre channel in grey; (f) MUPUS observations, with average represented in white (Spohn et al. 2015); (g) VIRTIS observations (Leyrat et al. 2016); (h) MIRO observations from this work in Seth region, with TI from the millimetre channel in dark green and TI from the submillimetre channel in light green (horizontal lines represent upper or lower bounds); (i) MIRO observations from this work in Aten region; (j) MIRO observations from this work in Ash region; (k) MIRO observations from this work in Ash region; (l) MIRO observations from this work in Aten region; and (m) VIRTIS observations from this work; the best-fitting value is indicated in white.

Acknowledgements

This project has received funding from the European Union’s Horizon 2020 research and innovation programme under grant agreement No 686709. This work was supported by the Swiss State Secretariat for Education, Research and Innovation (SERI) under contract number 16.0008-2. The opinions expressed and arguments employed herein do not necessarily reflect the official view of the Swiss Government. We also thank Lei Pan and Paul Springer, who work to maintain and improve the JPL database.

References

All Tables

All Figures

	Fig. 1 Temperature as a function of depth (on the surface and 1 cm, 2 cm, and 4 cm below the surface) for a TI of 80 JK−1 m−2 s−0.5, for September 15, 2014 (UT 00:00).
In the text

	Fig. 2 Temperature on the surface for various values of the thermal inertia in the range 5–320 JK−1 m−2 s−0.5, for September 15, 2014 (UT 00:00).
In the text

	Fig. 3 Radiance measurements from VIRTIS projected onto the shape model on 12 Sep. 2014 (left) and 15 Sep. 2014 (right).
In the text

Fig. 4 Simulated brightness temperatures (SBT; coloured lines), as functions of the electrical skin depth, for various values of TI in the range 5−320 JK−1 m−2 s−0.5 (SI units). The black dotted lines represent the MBT from the millimetre and submillimetre receivers for September 1, 2014 at 23:48. The millimetre SBT curves with TI values less than 10 JK−1 m−2 s−0.5 intersect the millimetre MBT line; and the submillimetre SBT curves with TI values less than 80 JK−1 m−2 s−0.5 intersect thesubmillimetre MBT line.

In the text

	Fig. 5 Same as Fig. 4 but for September 2, 2014. The results are given in Table 2.
In the text

	Fig. 6 Same as Fig. 4 but for September 12, 2014. The results are given in Table 2.
In the text

	Fig. 7 Same as Fig. 4 but for September 13, 2014. The results are given in Table 2.
In the text

	Fig. 8 Same as Fig. 4 but for September 15, 2014. The results are given in Table 2.
In the text

Fig. 9 Left panel: example of modelled radiance minus observed radiance at 4.95 μm on September 15, 2014 for TI = 5 JK−1 m−2 s−0.5. Increasing the TI makes the plot more symmetrical around zero. Right panel: mean value of model minus observation as a function of TI for two epochs and two wavelengths. We observe that increasing the TI beyond 80 JK−1 m−2 s−0.5 does improve the fit, and slightly worsens it for the case September 15, 4.75 μm. We therefore use this TI value as a lower limit in our analysis.

In the text

	Fig. 10 Difference between modelled rough radiances and VIRTIS measurements for September 12, 2014 (left panel) and September 15, 2014 (right panel). Although we still observe some scattering, the dispersion is far smaller than when not considering roughness, as in Fig. 9.
In the text

	Fig. 11 Left panel: regional map of the comet from El-Maarry (2015) with morphological units. Right panel: roughness parameter, i.e. the density of the half-spherical depressions projected onto the shape model. The density value is represented by the colour bar on the right.
In the text

Fig. 12 Half-spherical depression density overlaid on OSIRIS images, 3D projected onto the shape model with the same colour scheme as in Fig.11. The images are centred on approximately 130° longitude and 10° latitude. One can see that areas with similar morphologies share similar roughness. Visually smooth terrains (at the scale of OSIRIS data, i.e. 75 cm px−1) can display high thermal roughness, indicating that the spatial scale of this parameter is far below the resolution of our images.

In the text

Fig. 13 Summary of measured TIs (units: J K−1 m−2 s−0.5) from a variety of sources. Measurements in black represent previous published works and coloured bars represent this work. From Left to right, these measurements represent the TI from: (a) Spitzer observations (Lamy et al. 2008); (b) a thermophysical model applied to light curve observations (Lowry et al. 2012); (c) MIRO observations by Gulkis et al. (2015); (d) MIRO observations by Schloerb et al. (2015); (e) MIRO observations by Choukroun et al. (2015), with TI from the millimetre channel in black and TI from the submillimetre channel in grey; (f) MUPUS observations, with average represented in white (Spohn et al. 2015); (g) VIRTIS observations (Leyrat et al. 2016); (h) MIRO observations from this work in Seth region, with TI from the millimetre channel in dark green and TI from the submillimetre channel in light green (horizontal lines represent upper or lower bounds); (i) MIRO observations from this work in Aten region; (j) MIRO observations from this work in Ash region; (k) MIRO observations from this work in Ash region; (l) MIRO observations from this work in Aten region; and (m) VIRTIS observations from this work; the best-fitting value is indicated in white.

In the text