Thermal inertia and roughness of the nucleus of comet 67P/Churyumov–Gerasimenko from MIRO and VIRTIS observations

Aims. Using data from the Rosetta mission to comet 67P/Churyumov–Gerasimenko, we evaluate the physical properties of the surface and subsurface of the nucleus and derive estimates for the thermal inertia (TI) and roughness in several regions on the largest lobe of the nucleus.
Methods. We have developed a thermal model to compute the temperature on the surface and in the uppermost subsurface layers of the nucleus. The model takes heat conduction, self-heating, and shadowing effects into account. To reproduce the brightness temperatures measured by the MIRO instrument, the thermal model is coupled to a radiative transfer model to derive the TI. To reproduce the spatially resolved infrared measurements of the VIRTIS instrument, the thermal model is coupled to a radiance model to derive the TI and surface roughness. These methods are applied to Rosetta data from September 2014.
Results. The resulting TI values from both instruments are broadly consistent with each other. From the millimetre channel on MIRO, we determine the TI in the subsurface to be <80 JK−1 m−2 s−0.5 for the Seth, Ash, and Aten regions. The submillimetre channel implies similar results but also suggests that higher values could be possible. A low TI is consistent with other MIRO measurements and in situ data from the MUPUS instrument at the final landing site of Philae. The VIRTIS results give a best-fitting value of 80 JK−1 m−2 s−0.5 and values in the range 40–160 JK−1 m−2 s−0.5 in the same areas. These observations also allow the subpixel scale surface roughness to be estimated and compared to images from the OSIRIS camera. The VIRTIS data imply that there is significant roughness on the infrared scale below the resolution of the available shape model and that, counter-intuitively, visually smooth terrain (centimetre scale) can be rough at small (micrometre–millimetre) scales, and visually rough terrain can be smooth at small scales.


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 na-ture 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 di erent 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 e ciently 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 −1 m −2 s −0.5 for Hartley 2 and lower than 45 JK −1 m −2 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 −1 m −2 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 a ect the surface properties of the comet. Choukroun et al. (2015) also calculated TI from MIRO data and obtained values of 10-40 JK −1 m −2 s −0.5 and 20-60 JK −1 m −2 s −0.5 for the submillimetre and millimetre channels, respectively. In addition, Gulkis et al. (2015) found low TIs of 10-50 JK −1 m −2 s −0.5 from MIRO observations as Rosetta rst 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 −1 m −2 s −0.5 at the nal landing site, Abydos (Spohn et al. 2015). This represents the best-tting 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 −1 m −2 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 −1 m −2 s −0.5 (Lowry et al. 2012).
For the asteroid Lutetia, which Rosetta passed during its ight 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 −1 m −2 s −0.5 ), which increased with depth, and a 50% fractional surface coverage with hemispherical craters was required to t 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 reected 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 pro le. We focussed our analyses on data obtained in September 2014. This choice was motivated by the opportunity to compare our nal 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 di erent regions.
In the next section, we describe the MIRO and VIRTIS instruments. Then, the rst step towards estimating the TI involves assessing the datasets to nd suitable data for comparison (Section 3.1). The thermal model applied to the data from both instruments and used to calculate the temperature gradients is described in Section 3.2, and the details of the MIRO radiative transfer model are given in Section 3.3. The VIRTIS modelling is described in Section 3.4. The results from MIRO and VIRTIS are given in Sections 4.1 and 4.2, respectively.

Instruments
The MIRO instrument was a microwave spectrometer consisting of a 30 cm o set parabolic re ector 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 nearsurface 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 seconds and was calibrated in situ every 34 minutes 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 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 de ned as the required temperature of a black body, which lls 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 1 0.94 was applied to the brightness temperatures from the submillimetre channel as a result of the di erence in beam e ciency 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 eld 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 of 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 µm), 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 re ected ux 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).

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 identi ed 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 horizontal resolution of approximately 20 m, and the SPICE soft-ware (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 (angle between 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 di erent regions of the nucleus but occasionally they both looked at the same location.
The two instruments acquired data in slightly di erent 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 m 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 VIR-TIS band. Since both instruments operated at di erent wavelengths and spatial resolutions, the roughness is treated di erently 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 ve of the best times; these times were selected because they all have wellconstrained MIRO and VIRTIS observations, good viewing geometries, and cover di erent areas of the nucleus. These times (speci ed in UT) are given in Table 1 along with the corresponding observed MIRO brightness temperatures. These observations provide data on ve areas located on di erent 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 the Ash 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.

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 pro le 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 onedimensional thermal model using the heat equation where ρ is the bulk density in units of 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 de ned as where I has SI units of [J K −1 m −2 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 −3 (Jorda et al. 2016) and assume the heat capacity to be 500 J kg −1 K −1 , estimated from Robie et al. (1970). For the surface boundary, the following condition is used to solve the heat equation where S , the solar ux at 1 AU, is equal to 1370 W m −2 , 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 , 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 rst 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, 2014, September 13, 2014, and September 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 di erent values of TI (see Eqs. (2 -4) and Groussin et al. (2013)): 5, 10, 20, 40, 80, 160 and 320 JK −1 m −2 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 sec 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 pro les 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 dT dt = 0. 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 hours to run a complete set of TI values at a given date.
As an output, we obtain, for a given date, the temperature pro le 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 −1 m −2 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 −1 m −2 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.

Radiative transfer model for MIRO data
The MIRO radiometer measures antenna temperatures at millimetre and submillimetre wavelengths, described in Section 2, from which brightness temperatures for the observed areas are derived. Using the thermal outputs from Section 3.2 as inputs to a radiative transfer model, simulated brightness temperatures (hereafter SBT) are compared to measured brightness temperatures (hereafter MBT) for the di erent 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 eld is applied, respectively. Their values are relative to that of a vacuum, which is equal to 8.85×10 −12 F m −1 . 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 pro le of the subsurface and the permittivity D. Marshall et al.: Thermal inertia and roughness of 67P/CG 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 e ective spectral brightness intensity BT ν,e f f and the surface emissivity, Υ, is as follows where ν is the frequency of the observation, θ e , the facet emission angle de ned 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 re ectivity coe cients. The e ective 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 e ective spectral brightness intensity computed for a given facet, assuming a scatterfree homogeneous layer under conditions of local thermodynamic equilibrium, is expressed as follows (Ulaby et al. 2014) 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 eld can penetrate into the material.
Regarding the MIRO observations, we expect the submillimetre receiver to be sensitive to the thermal radiation integrated over the top subsurface of the layer and the millimetre Article number, page 5 of 15 A&A proofs: manuscript no. 33104corr 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).

Radiance model for VIRTIS data
With the VIRTIS data, we focus our study on only two speci c dates in Table 1: September 12 and September 15, 2014. Fig. 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 re ected solar radiation. The observed areas cover the regions of Aten, Babi, Khepry, and Imhotep, as de ned 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 signi cant 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 rst convert the modelled surface temperature into radiance, I(λ, T ), by combining re ected solar spectral radiance and emitted IR radiation as follows where α and are the constant geometric albedo and emissivity (respectively 0.05 and 0.95), F the solar ux modulated by a re ectance function (the Lommel-Seeliger Law), and B the black body radiation at wavelength λ and temperature T . We neglected the directionality of epsilon, which is a ected 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 where h, c, k are Planck's constant, the speed of light, and the Boltzmann constant, respectively. The value R S un is the sun radius, r h the heliocentric distance. The values µ 0 and µ are the cosines of incidence and emission angle.
Finally, the emitted IR radiation is the Planck function at T and λ, is We applied our model to all facets covered by the VIRTIS data set, except those in shadows for which a more complex photometric model should be considered, accounting for multiple scattering. Shadows were calculated with a custom ray tracing engine 1 widely used to simulate Rosetta images, and calibrated against observations.

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;Groussin et al. 2013;Keihm et al. 2012;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 in uenced by the roughness, and so the energy received is a function of the relative position of the observer and 1 shapeViewer: www.comet-toolbox.com 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 rst 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 nal 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 sections 2 and 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 e ects. In this approach, we rst 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 com-pared to OSIRIS images of the same region, as they have typically much higher spatial resolution than VIRTIS data. Wenally 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 inputs would 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 di erent geometries. However, since this is not always available in our dataset, we choose a less constrained approach, as outlined above. Fig. 4, 5, 6, 7, and 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 to 320 JK −1 m −2 s −0.5 , for the ve dates given in Table 2. As stated earlier, we constrain the TI when SBT intersects MBT for electrical skin depths less than ∼3.0 m and ∼1.0 m for the millimetre and submillimetre cases, respectively.

MIRO results
For example, in Fig. 4 in the millimetre channel, the SBT with TI = 5 JK −1 m −2 s −0.5 intersects the MBT but TI = 10 JK −1 m −2 s −0.5 does not. We hence derive that the TI in the millimetre channel for September 1 is <10 JK −1 m −2 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 insu cient, 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 di erent observations performed with the millimetre receiver are consistent with a TI less than 80 JK −1 m −2 s −0.5 , while the value for September 1, 2014 implies that the TI should be less than 10 JK −1 m −2 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 −1 m −2 s −0.5 , but the other dates imply higher values of 160 and 320 JK −1 m −2 s −0.5 are also possible for their observations. The results are compiled in the Table  2.

VIRTIS results
Disentangling TI from roughness e ects is a complex problem for which we do not always have enough constraints. To reiterate, we attempt to t rst the TI as it relates directly to the physical properties of the material. We then add the surface roughness as a second parameter e ectively describing which subpixel area contributes to the ux that reaches the detector.
Therefore, the rst 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 −1 m −2 s −0.5 ). We found that the modelled ux is almost always below the measured ux (see Fig. 9). As we increase TI, the average ux 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 t is achieved for a TI larger than 80 JK −1 m −2 s −0.5 . We therefore use this TI value as a lower limit in our analysis. Higher values lead to an equally good t, suggesting that TI is not the main reason for the dispersion we observe in Fig. 9, although the t is slightly worse for TI values of 160 and 320 JK −1 m −2 s −0.5 in the case of September 15, 4.75 µm.
Having selected TI=80 JK −1 m −2 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 de ne the corrected ux for a rough terrain as I r = c r × I, with c r a function of incidence and emission angles de ned 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 in uence of opening angle is less important than that of the crater density for the wavelengths considered. Therefore, we x the opening angle to a value of 180 (which leads to maximal c r ) and nd the crater density that minimizes the quantity I r − I m : the di erence between the calculated radiance, I r , and the VIRTIS measured radiance, I m . The crater density is allowed to vary between 0 and π 2 √ 3 , the maximum density for packing equal circles on a two-dimensional plane. This approach e ectively gives us a lower bound for the crater density on each facet in the region of interest.
We note that in principle, one should nd the best parameters by varying both TI and roughness in parallel when getting the minimum least-squares t. However, this is a large parameter space and running the complete inversion goes beyond the scope of this work. Fig. 10 shows an example of the corrected radiances for two di erent dates. The most remarkable e ect 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 ner spatial resolution for both parameters would certainly improve the ts, 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 minutes in this work, during which the viewing geometry changes as the comet rotates and the spacecraft moves. This also a ects 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 nd 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 di erent textures. The northern part is rougher with many metre-size 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 rst, 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 . This region has a very 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-size scale, display at 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 at surfaces . Further study is required to link the measured roughness to these morphological processes.

Thermal inertia
All of the results obtained from the MIRO millimetre observations imply that the TI has a low value, <80 JK −1 m −2 s −0.5 . The result from September 1 even implies that it should be Fig. 9. Left: 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: 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 t, 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. Constraining the TI from the subsurface temperature observations is sometimes di cult, however. For the measurements on September 12, 2014 and September 15, 2014, the submillimetre observation cannot constrain the TIso well (5-320 JK −1 m −2 s −0.5 and >20 JK −1 m −2 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 −1 m −2 s −0.5 .
In comparison, VIRTIS work implies a best-tting value of 80 JK −1 m −2 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 di erence. 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 Figure 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 −1 m −2 s −0.5 (Leyrat et al. 2016). Neglecting spatial and depth variations, a bulk TI value in the range 30-50 JK −1 m −2 s −0.5 would t 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.  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) can display high thermal roughness, indicating that the spatial scale of this parameter is far below the resolution of our images.
Taken together, the Rosetta measurements are consistent with a low TI (≤170 JK −1 m −2 s −0.5 ), and there is notable overlap in the error bars of these results.

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 de ne the surface roughness properly, our approach shows that one can use IR radiance measurements and a simpli ed roughness