© The Authors 2025

1 Introduction

The Martian Moons Exploration mission MMX (JAXA) will be launched in October 2026 towards the Martian system. After its arrival in 2027, the probe will closely study the system, specifically looking at its two moons: Phobos and Deimos (Kuramoto et al. 2022). The aim of this sample return mission is to provide new insights that would help to lift the mystery around the origins of the Martian moons.

Phobos and Deimos share strong spectral similarities with D-type asteroids (Thomas 1979), and hence seem to have been captured by Mars’ gravity. However, their orbital parameters (low inclination and eccentricity) suggest a formation resulting from a large impact on Mars (Craddock 2011).

To support the MMX mission, the onboard MMX InfraRed Spectrometer (MIRS) will perform observations of the surfaces of Phobos and Deimos at a resolution of up to few meters/pixel and will record the reflected flux (Barucci et al. 2021) in the [0.9–3.6] μm range. In parallel with helping in the selection of the sampling site on Phobos, MIRS will provide spectral data to build a global map of the moons at a high spectral and spatial resolution and will allow a close analysis of the surface composition. Such precise and global analysis may hold the keys to unravelling the origins of Phobos and Deimos (Barucci et al. 2021) as it will provide new insights into their surface composition and evolution. Specifically, MIRS is designed to search for hydrated minerals and organic materials by looking at the structural hydroxyl absorption features at 2.7 μm (and their associated overtones at 1.4 and 1.8 μm) and the aliphatic–aromatic organic materials near 3.5 μm that would support the capture hypothesis. In the giant impact scenario, the precursor materials would be heated to peak temperatures exceeding 1000 K (Hyodo et al. 2017), likely resulting in the complete loss of water.

The spectral range of MIRS is also compatible with the search of anhydrous mineral signatures (0.9–1 and 2 μm) that would instead support the impact hypothesis (Craddock 2011).

Because of the low near-infrared (NIR) albedo at the Phobos and Deimos surfaces (Fornasier et al. 2024), at a heliocentric distance of 1.5 AU, the solar flux is sufficient to raise the surface temperature in the range ([130–353] K, Giuranna et al. (2011)). Therefore, its thermal emission increases above the reflected part of the spectra in a significant spectral range of MIRS data (Fraeman et al. 2014; David et al. 2024). This contribution will increase the measured flux starting from 2.5 μm (Fraeman et al. 2012, 2014) and change the shape of any measured spectra, potentially hiding absorption features, and will therefore greatly complicate the analysis and the spectral data interpretation as absorption bands associated with hydrated minerals and organics are found around ≃2.7 μm and ≃3.2 μm, respectively. During the warmest hour for a given surface, the higher measured flux after 2.5 μm will increase the spectral slope specifically in this region. Moreover, with surface temperatures varying across the Phobos surface, the spectral location and shape of this thermal tail will change within the spectrum. Therefore, a correction for the thermal emission contribution is mandatory (Clark et al. 2011) to retrieve spectral parameters specifically the spectral slope and potential absorption features, and thus interpret spectroscopic data. David et al. (2024) has already proposed a singular iterative method to retrieve the temperature from a measured spectrum based on Clark et al. (2011) and its ‘thermal tail’.

In this work, we present a different but complementary approach based on geometry, thermal modelling, and an estimation of the thermal properties of the surface to compute the surface temperature at a global scale on the moons. Although it depends strongly of the knowledge of the surface, an issue that will be fixed for the MMX mission by independent observations from the other instruments on MMX, the method would be insensitive to measurements artefacts and will be a support for more detailed analysis of the future MIRS spectra.

Fig. 1 Simple schematic of the model pipeline. The surface incoming flux is first computed in order to derive the surface temperature based on the ephemeris, the Mars albedo mosaic from Christensen et al. (2001), the shape model of the targeted moon, and its albedo. The self-heating is computed iteratively after a first run of the thermal model until the convergence is reached, reinjecting the previous estimation to compute the new surface temperature. The thermal model takes the grain size Dg, the surface porosity p and density ρ, and the non-porous thermal conductivity km and specific heat Cp as inputs. Eventually the surface thermal emission is computed for a rough surface knowing the single scattering albedo (SSA) ϖ; the surface Hapke’s average slope angle $\overline{\theta }$; the opposition effect amplitude B0 and angular width parameter hs; the phase function asymmetry parameter b; and the geometry of observation regrouping i, e, and ϕ for the incidence, observation, and phase angle, respectively. This thermal emission is to be used as a reference to remove the emission component in the spectra.

2 Models

We established a simple pipeline to compute the surface emission in order to remove its contribution in the NIR spectra. We detail this pipeline below (see Fig. 1).

2.1 Estimating the absorbed flux

The first step in determining the surface temperature of the Areolian moons is to estimate the total flux reaching the surface and being absorbed by the surface. Therefore, we adapted a simulation already existing for the Cronian system (Le Gall et al. 2023; Bonnefoy et al. 2020; Leyrat et al. 2015, 2011; Ferrari & Leyrat 2006) describing how light propagates into the Martian system, accounting for the reflections between Mars and its moons. For each instant of the simulation (every 20 minutes on a seasonal cycle) the flux emitted by the Sun at the distance of Mars is computed. Thanks to the NAIF SPICE toolkit (Acton 1996; Acton et al. 2018) the geometry of the problem (i.e. the distance between objects, illumination, and reflection angles) is also computed, with respect to the position of the bodies at each time of the simulation. We mapped the moons’ surfaces with a 3 × 3°.pixel⁻¹ (≃600 m ×600 m at the Phobos equator) grid projected onto the shape model (Willner et al. 2014) of the moons. Using the shape model allows the right incidence angle for the illumination to be computed at each coordinate of the moon and is a way to take into account the topography, and thus the associated large-scale shadowing.

The absorbed flux F_abs is computed at each time of the simulation as $$\begin{array}{rl}{F}_{abs}=(1-a_{moon}^{VIS}).[& F_{\odot }^{dir}.cos({\theta }_{i,\odot })+F_{\odot }^{ref}.cos({\theta }_{i,♂})]\\ & +(1-a_{moon}^{IR}).(F_{♂}^{emi}).cos({\theta }_{i,♂}),\end{array}$$(1)

θ_i,source being the incidence angle to the local normal of the surface from the light source (Sun or Mars). It takes into account the solar flux $F_{\odot }^{dir}$ at the moon heliocentric distance, the received solar flux on Mars (Mars-shine) reflected towards the moon $F_{\odot }^{ref}$, and the flux emitted by Mars at its equilibrium temperature $F_{♂}^{emi}$, and the Bond albedos in the visible and infrared $a_{moon}^{[Vis/IR]}$ (Fornasier et al. 2024).

Direct and reflected solar fluxes strongly depend on the heliocentric distance and the eclipses. The moon orbital period defines the simulation time step at which the absorbed flux is computed. To sample correctly the eclipses (3 per 25 h for Phobos), we set the time step to ≃22 min. The sub-Mars hemisphere ([315, 45]° E) is always irradiated by Mars thermal emission. Mars dust storms influence $F_{\odot }^{ref}$ and $F_{♂}^{emi}$ as its surface visible albedo increases by 15% (Christensen 1988), thereby increasing the reflected light. However, the surface is no longer visible in this condition, and the Mars surface thermal emission is set to 0: $$F_{\odot }^{dir}={\sigma }_{SB}.T_{\odot }^4.{\left(\frac{R_{\odot }}{d_{moon}}\right)}^2,$$(2) $$F_{\odot }^{ref}={\sigma }_{SB}.T_{\odot }^4.{\left(\frac{R_{\odot }}{d_{♂}}\right)}^2.a_{♂}^{VIS}.2\pi .{\left(\frac{R_{♂}}{h_{moon}}\right)}^2,$$(3) $$F_{\delta }^{emi}={\sigma }_{SB}.T_{♂,eq}^4.(1-a_{♂}^{IR}).2\pi .{\left(\frac{R_{♂}}{h_{moon}}\right)}^2.$$(4)

Here T_⊙ and T_♂,eq are respectively the temperatures of the Sun and Mars, d_moon and d_♂ are the moon and Mars heliocentric distances, h_moon is the moon altitude over the Mars surface, R_⊙ and R_♂ are the Sun and Mars radii, and a_♂ is the Mars albedo. This latest is retrieved from the Mars Global Surveyor (MGS) Thermal Emission Spectrometer (TES) global mosaic (Christensen et al. 2001).

However, we did not consider any surface albedo variation of the moons since only a low resolution albedo map exists for Deimos (Thomas et al. 1996) and is only partial for Phobos (Simonelli et al. 1998; Fornasier et al. 2024). Albedos were therefore fixed to $a_{Phobos}^{VIS}=0.0161$ (averaged between green and black filter from Fornasier et al. (2024)) for the visible part and $a_{Phobos}^{IR}=0.02$ ($a_{Phobos}^{IR}=1-e_{Phobos}^{IR}$ from Giuranna et al. (2011)) for the infrared part of the spectra.

The flux in regions where incoming light is blocked by large-scale topographic features represented in the shape model was assumed to be zero. However, occlusion due to unresolved surface roughness was not accounted for, which may result in an overestimation of flux near the terminator. Nevertheless, this would not be a significant issue for MIRS because the measurements are likely to be limited to 60° of incidence.

Fig. 2 Simple drawing gathering all flux considered in the simulation. Solar flux (direct and reflected by Mars) is indicated in yellow. Thermal emission (Mars and self-heating) is represented in red.

2.2 Self-heating

In order to keep the model simple and to avoid too long computation times, we did not consider the self-illumination (reflections from the moon surface towards other parts of the moon) because of the low surface albedo (<5%). However, the self-heating consisting on emitted light at the surface on convex shapes cannot be neglected and is estimated.

We thus developed a simple model that computes the self-heating at the equilibrium temperature T_eq (taken as the average temperature on the year after several computations of the thermal model), in which we reinject the self-heating from the previous iteration. The incident flux on a map-pixel i coming from all elements illuminating it is $$F_{SH}(i)=(1-a_{moon}^{IR}).\sum _{j\ne i}\frac{\nu (i,j).(1-a_{moon}^{IR}).{\sigma }_{SB}.T_{eq}^4.cos({\theta }_{i,j}).s_j.A_i}{d_{i,j}^2},$$(5)

where ν(i, j) is either 1 or 0 if the surface elements are viewing each other, θ_i,j represents the elements normal angular separation, s_j is the surface of the i_th pixel, A_i is the surface area of the i_th pixel, and d_i,j is the estimated distance between the two elements computed on the Phobos ellipsoid.

To ensure a good convergence of the self-heating effect on the temperature, we computed the self-heating several times as a test until the average flux in the region converged. At each iteration, we reused the previous estimation of the self-heating in the computation of the new surface temperature, which was then reinjected for the new estimation of the self-heating.

At the end of the simulation, the flux absorbed by the moons was computed over an entire revolution around the Sun, which was required for thermal modelling.

2.3 Computing the surface temperature

2.3.1 A 1D thermal transfer

The surface temperature was computed over a full seasonal cycle using a thermal model developed for the thermal emission of particles inside Saturn’s ring by Ferrari & Leyrat (2006), which was then adapted to asteroids and icy satellites (Leyrat et al. 2011, 2012; Le Gall et al. 2014; Leyrat et al. 2015; Bonnefoy et al. 2020; Le Gall et al. 2023) and now for the Areolian moons.

The model solves numerically the time-dependent heat equation in 1D and considers heat transfer by conduction and radiation on the vertical axis only (from the surface towards the sub-surface) using a semi-implicit algorithm (the Crank-Nicholson method).

At the surface, the imposed limit condition is a radiative equilibrium (Eq. (6)) between the outgoing (thermal emission towards space) and incoming fluxes (solar direct, Mars-shine, Mars thermal emission, and self-heating at the surface). The temperature is considered constant at large depths Eq. (7) (deeper than ten times the seasonal thermal skin depth δ_th,s; adiabatic condition) and because of their size and the low energy dissipated in the moons from the Mars gravitational tides, we do not consider any endogenous heat sources (radioactive decay, tidal heating): $$F_{abs}=(1-a_{moon}^{IR}){\sigma }_{SB}.T^4-{\kappa }_{eff}\frac{\mathrm{\partial }T}{\mathrm{\partial }r},$$(6) $${\frac{\mathrm{\partial }T}{\mathrm{\partial }r}|}_{r>10{\delta }_{th,s}}=0.$$(7)

Moreover, the model only considers constant thermal properties in the sub-surface, which is the main limit to this model. Because the model considers only heat transfer along the vertical axis, the maximum resolution at which the surface temperature can be estimated is approximately 10 seasonal skin depths (δ_th,s), which given their thermal properties, lies around 10 m on Martian moons.

The two main parameters controlling the heat transfer in the sub-surface are the effective thermal conductivity κ_{ef f} and specific heat capacity C_p of the material at the surface. κ_{ef f} is strongly linked to the texture of the material (namely the grain size D_g and the porosity p) as the heat transfer relies on either conduction or radiation, depending on the texture.

2.3.2 A model for the surface thermal conductivity κ_eff

The thermal conductivity of the surface κ_{ef f} is determined following Eq. (6) from Sakatani et al. (2017), assembling the radiative and conductive terms. This relation between the texture and the thermal conductivity of the medium has been tested under vacuum conditions with glass beads as a model material for powders, for porosity ranging from 30 to 86.7% (Sakatani et al. 2017). We chose this model as it had been used in the case of microgravity for the thermal study of the Ryugu asteroid (Sakatani et al. 2018; Wada et al. 2018), and validated experimentally for high porosity (up to 86.7%) and a wide range of grain sizes up to 1000 μm and down to a few micrometres. This model assembles the two mechanisms of heat transfer: solid conductivity κ_solid and radiation κ_rad.

The solid conductivity κ_solid is determined as $${\kappa }_{solid}=\frac{4}{{\pi }^2}{\kappa }_m(1-p)C\xi \frac{r_c}{D_g}.$$(8)

Here κ_m represents the non-porous material thermal conductivity, p is the porosity, and C is the coordination number from Suzuki et al. (1980) $$C=\frac{2.812(1-p{)}^{-1/3}}{f^2(1+f^2)},$$(9)

where f = 0.07318 + 2.193p − 3.357p² + 3.194p³.

In Eq. (8), r_c describes the contact radius and is defined by $$r_c={\left[\frac{3(1+{\nu }^2)}{4E}(F+\frac{3}{2}\pi \gamma D_g+\sqrt{3\pi \gamma D_{g}F+\left(\frac{3}{2}\pi \gamma D_g\right)})\right]}^{1/3},$$(10)

where ν is Poisson’s ratio, E the Young modulus, F the external compressive force, and γ the surface energy (Johnson et al. 1971). The compressional stress expressed as $F=\frac{2\pi D_g^2}{\sqrt{6}(1-p)}\times {\rho }_0(1-p)gz$ (Sakatani et al. 2017) has been neglected since its contribution in the first centimetres would be small and the thermal model does not consider variations in depth.

The radiative contribution to the thermal conductivity is expressed as $${\kappa }_{rad}=\frac{4ϵ}{2-ϵ}{\sigma }_{SB}L{T}^3,$$(11)

where ε, σ_{S B}, and T are respectively the grain emissivity (here taken as equal to the surface emissivity), the Stefan-Boltzmann constant, and the temperature; and L is a factor describing the geometrical clearance (mean free path for radiation) in the surface, and is expressed as $L=2\zeta {\left(\frac{p}{1-p}\right)}^{1/3}{D}_g$, as described in Piqueux & Christensen (2009), ζ being a scaling factor linked to the particle size as $\zeta =0.68+\frac{7.6\times {10}^{-5}}{D_g}$ in Sakatani et al. (2018) and Wada et al. (2018).

2.4 Estimating the flux emitted from the rough surface

Following the method in Davidsson et al. (2015) and Davidsson & Rickman (2014), the thermal emission of a rough surface can be assimilated to the superposition of the thermal emission of a large number of facets having their own temperature due to their relative orientation to both the light source and the observer. In the case of the self-heating, this large number of facets also increase the multiple-scattering at the surface, which increases the flux budget for each facet. In addition, at large incidence angles, the facets project shadows onto their neighbours, decreasing the absorbed flux budgets and thus the temperature. Finally, the emission of the rough surface is also dependent on the geometry since facets will be hidden by their neighbours due to the roughness.

Following this approach, we simulated the emission of an unresolved roughness as the sum of the emission of N_sub sub-element within our map-pixels. Thus, we computed a large number of virtual sub-elements (N_sub = 100) within the map-pixel, each tilted at an angle θ_sub with respect to the map-pixel normal. This tilting angle θ_sub is randomly drawn from the realistic facet angle distribution computed from Hapke’s average slope parameter $\overline{\theta }$, derived from the High Resolution Stereo Camera (HRSC) observation of Phobos (Fornasier et al. 2024). This angle modulates the received flux F.cos(θ_sub,i) in each sub-element and ensures a realistic heating due to the unresolved roughness.

To model the shadows projected by some sub-elements onto the others due to the roughness in the map-pixel, we added a term to shut down the illumination of a fraction of the sub-elements at each time step. The proportion of sub-elements in the shadows is estimated using the shadowing function S (i, e, ϕ) (Equation (12.28) from Hapke 2012). S (i, e, ϕ) depends on the geometry of illumination and observation (incidence angle i, emergence angle e, and phase angle ϕ). We computed S (i, e, ϕ) individually for each map-pixel at each time step in the simulation. The emergence angle was set to e = 0, as our goal was to estimate the number of non-illuminated sub-elements within a map-pixel only due to the unresolved roughness. However, the incidence angle i = i(t, ((λ^o, δ^o))) of a given map-pixel depends on its coordinates (λ^o, δ^o) and the local time t; the phase angle ϕ is set to ϕ = i since the shadowing function is computed in the nadir direction.

Then, for each sub-element, the temperature is computed using the thermal model described earlier.

Although it is possible to do so, we do not vary the surface texture (grain size, porosity) and composition properties from one sub-element to another; instead, we assume that the surface is homogeneous within a map-pixel.

Hence the emission of each sub-element is computed as $$I_i=H(\varpi )B\left(T_i\right)B_{SHOE}(e_i,h_S),$$(12)

where H(ϖ) is the Chandrasekhar H function evaluated at the single scattering albedo (SSA) ϖ, B(T_i) represents the black-body emission of the sub-element i at the temperature T_i, B_{S HOE}(e_i, h_S) is Hapke’s shadow hiding opposition effect (SHOE) function evaluated at the emission angle e_i of the sub-element i and at the angular-width parameter h_S.

Then the map-pixel emission is simply constructed as the sum of each sub-element weighted by the density of probability of the surface slope distribution a(i) giving $I_{th}={\sum }_i^N{I}_{i}a(i)$.

A strong difference with the work from Davidsson et al. (2015) and Davidsson & Rickman (2014) lies in the absence, in our case, of taking the projection of shadows, which will have an important impact on temperature at large incidence angles. Additionally, we take into account the shadowing in emission in a simpler way using the Hapke’s SHOE function.

To take into account the reflective part in the spectrum, we compute the reflectance spectrum of the surface using the Hapke’s parameters derived from the observations of Phobos by HRSC.

3 Removing the surface thermal emission from spectra

3.1 Creating synthetic data

As a test of our model and in order to remove the contribution of the surface thermal emission in the spectra, we created a set of synthetic data similarly to David et al. (2024).

For each map-pixel, we simulated a sub-scale roughness similar to the one described above (see Section 2.4) to simulate the emission of a rough surface. We used the output of a thermal model to compute the surface temperature of N = 100 sub-elements, with an arbitrary value for the grain size D_g, the porosity p, the non-porous material specific heat C_p, and the thermal conductivity κ_m, but consistent with the Phobos. Then, we computed the map-pixel thermal emission as a superposition of the N black body (see Section 2.4).

To this black-body spectrum I_th, we added a reflective component from Hapke’s model using parameters derived from the HRSC observation (Fornasier et al. 2024). The reflectance r_Hapke computed from Hapke’s model was then multiplied by the solar spectrum S_⊙ at Mars, and added to the thermal emission spectrum. The Hapke’s parameters that were required to compute the reflectance (namely the SHOE opposition amplitude B₀, the angular half-width h_s, the asymmetry parameter g, and the mean slope $\overline{\theta }$) were derived from the HRSC observation. Only the SSA ϖ was from the CRISM data (Fraeman et al. 2014). The new spectrum was then converted into I/F by dividing again by S_⊙ and multiplying by π.

Some Gaussian noise G_noise was then added to the spectrum to simulate the future MIRS observation S/N, the noise is taken as 1% of the I/F value at each wavelength.

The final synthetic spectrum was computed as $$I/F=\pi .\frac{S_{\odot }.r_{Hapke}+I_{th}}{S_{\odot }}+G_{noise}.$$(13)

3.2 Test with CRISM data

To test this correction on the CRISM data, we used the cube FRT00002992_03 acquired on October 23, 2007, around 1h18 UTC. These data contain only the I/F and no correction for the observation geometry was performed.

In this observation, we defined a region of interest (RoI) close to the sub-Mars hemisphere on the Phobos surface. The choice of this RoI relies on the apparent smoothness of the area (no craters) and the constant phase angle ϕ ≃ 37°, emission angle e ≃ 5°, and incidence angle i ≃ 30°.

We used our pipeline (our flux simulation on October 23, 2007, the thermophysical model, and the thermal emission model) to compute the I/F in this RoI at the specific time of the observation. We computed the surface temperature using parameters derived from HRSC photometric observations (Fornasier et al. 2024), the non-porous material specific heat C_p = 760, and the thermal conductivity κ_m = 1.6 (to simulate chondritic composition (Opeil et al. 2010)), but letting the grain size D_g and the porosity p be free parameters. Once again, the Hapke model was used to compute the reflective counterpart of the spectrum.

4 Results

4.1 Flux variations at the surface of Phobos

Fig. 3 presents an output of the flux model absorbed by the Phobos surface on October 23, 2007 (at 01^h17^min UTC).

At the date of the simulation, the sub-solar point was lying at −20°N : 326°E, east of Stickney crater Fig. 3. The maximum flux absorbed by the Phobos surface at this date reaches 550 W m⁻² on Stickney’s east rim and west inner wall.

Over the year, the instantaneous absorbed flux by the Phobos surface reaches a maximum of 718 W m⁻² at the equator from external sources (the Sun and Mars). Mars-shine and Mars thermal emission contributes a maximum of 77 W m⁻² on the sub-Mars hemisphere to this illumination. The self-illumination in the near-infrared by Phobos itself however adds at most an additional 34 W m⁻² specifically on convex surfaces.

Fig. 3 Simulation of the fluxes absorbed by the surface of Phobos on October 23, 2007. Stickney is the large crater on the right (310°E). This simulation shows the total absorbed flux on Phobos (direct and reflected solar flux, Mars thermal emission and self-illumination). The rectangle in black is the region of interest (RoI) where we plotted the absorbed flux over the complete revolution in Fig. 4. The black pentagon represents the position of the sub-solar point at the time of the simulation.

Fig. 4 In the RoI (black square in Fig. 3), the absorbed flux reaches 600 W.m−2, while on Phobos the maximum is 713 W m−2. A minimum of 35 W m−2 is constantly illuminating the Mars sub-point due to the Mars thermal emission.

4.2 Surface temperature variations

A temperature global mosaic computed for Phobos October 23, 2007, at the surface with our model is presented in Fig. 5a. We set the grain size to D_g = 500 μm; the porosity to p = 50%, according to the porosity derived for Phobos in Fornasier et al. (2024); the bulk thermal conductivity to κ_bulk = 1.6 W m⁻¹ K⁻¹ to be consistent with chondritic material, in the range [1, 2] W m⁻¹ K⁻¹ (Opeil et al. 2010; Yomogida & Matsui 1983) measured at 200 K; and the specific heat to C_p = 760 J kg⁻¹ K⁻¹ (Kuzmin et al. 2003). Although Phobos surface temperatures, specifically close to the sub-solar point, are expected to be higher, eventually reaching more than 300 K, Opeil et al. (2010) stated that thermal conductivity and specific heat are somewhat constant for this temperature range. Computed at Phobos equilibrium temperature (T = 267.2 K), this gives an effective thermal conductivity of κ_{ef f} = 2.710⁻² W m⁻¹ K⁻¹, consistent with values from Kuzmin et al. (2003) where conductivity is estimated to be in the range [7.5 × 10⁻³; 2.3 × 10⁻²] W m⁻¹ K⁻¹. The bulk density for Phobos was set to 1873 kg cm⁻³ (Pieters et al. 2014; Rosenblatt 2011). Overall, the thermal inertia computed from these values and used in this work is $I_{th}=\sqrt{{\kappa }_{eff}.(1-p).\rho .C_p}=86\mathrm{J}{\mathrm{s}}^{1/2}{\mathrm{k}\mathrm{g}}^{-1}{\mathrm{m}}^{-2}$, which is consistent with but slightly higher than that in previous works: 20–40 J s^1/2 kg⁻¹ m⁻² from Kührt et al. (1992) and 25–85 J s^1/2 kg⁻¹ m⁻² from Lunine et al. (1982).

As expected, an atmosphereless object where the Sun’s illumination drives the temperature, the warmest point is close to the sub-solar point −20°N : 326°E. Because of the orientation of Phobos with respect to the Sun, the boreal pole regroups the coldest area at the moon’s surface as the simulation date is set in the late northern winter (Ls = 334°).

Retrieved temperatures are consistent with previous works (Giuranna et al. 2011; Lunine et al. 1982), with a range of [130–300] K depending on the season and thus on the heliocentric distance. There is a strong dependence on the surface textural (D_g, p) and thermal properties (specific heat capacity C_p, thermal bulk conductivity κ_m).

Fig. 5 Temperature map at Phobos on 23 October 2007 at 1h17 UTC. The sub-solar point is located on the sub-Mars hemisphere, where the warmest temperatures are recorded (288 K with the given thermal and surface properties). The rectangle in black is the region of interest presented in Fig. 3. The black pentagon represent the position of the sub-solar point at the time of the simulation.

5 Discussion

5.1 The Mars effect

As illustrated in Fig. 4, Mars-shine and thermal emission contributes to a maximum of 77 W m⁻² at sub-Mars hemisphere when Phobos is aligned between the Sun and Mars. This constant illumination of the sub-Mars hemisphere (to 35 W m⁻²) due to Mars thermal emission cannot be neglected in addition to the 42 W m⁻² of the Mars-shine. Over the year, it raises the temperature by almost 7 K for the given surface texture and thermal properties, on the sub-Mars hemisphere.

Although 7 K may appear to be a relatively minor quantity, it nevertheless has a significant impact on the spectral amplitude and position of the thermal tail. By simulating the thermal emission spectra of black-bodies at different temperatures (285, 290, 295, and 300 K; Fig. 6b)), we found that the position of the start of the thermal tail varies by 0.1 μm every 5 K. Similarly, in the craters, the additional 35 W m⁻² illuminating the craters’ inner walls increases the temperature by about 7.5 K for the given surface texture and thermal properties.

Moreover, both of these additional flux sources present the importance of other sources in the correction of the surface thermal emission in order to subtract its contribution from the future MIRS/MMX near-infrared spectra.

Fig. 6 Thermal emission of a rough surface (computed from the Hapke mean slope parameter $\overline{\theta }=24$ for Phobos) showing the strong link with temperature of the thermal tail position.

5.2 Application to data: Comparison with black-body fitting

In this section, we develop the main goal of the paper, which is the removal of the thermal emission component in the reflectance spectra. We compare the results of our work to another method (black-body fitting) developed for the future MIRS data in David et al. (2024).

5.2.1 Removing the thermal emission in synthetic data

Following the protocol to produce synthetic data in Section 3.1, we constructed a set of data (see Fig. 7) simulating several observations to test the robustness of our thermal emission removal method. Thus, we simulated the I/F as observed on a complete day – 20 different local time around October 23, 2007, 1h17 UTC – in the RoI at eight different geometries of observation, varying the phase angle ϕ.

These synthetic data were produced for a grain size D_g = 500 μm and a porosity of p = 50%. Thermal properties κ_m and C_p were set respectively to 1.6 W m⁻¹ K⁻¹ and 760 J kg⁻¹ K⁻¹ to keep a chondritic material composition. These synthetic data are displayed in Fig. 7.

To test our method, we ran our thermal model at the date of the observation (October 23, 2007) for each observation time and each geometry, in order to perform a fit of the synthetic data leaving the grain size D_g and the porosity p as free parameters.

The corrected spectra presented in Figs. 7c and d show a very good removal of the thermal emission contribution for each local time and geometry of observation, the average residual slope being lower than 5 × 10⁻⁴ μm⁻¹.

As presented in Fig. 8, we performed another test by adding absorption bands at 1.4, 2.7, and 3.4 μm, which represents potential absorption respectively from water ice, hydrated minerals (as proposed in Fraeman et al. 2014), and organic materials in order to test the sensitivity of the method to the presence of a band and to compare the variation of the band. The relative depth of the bands vary within the range [10–50]. The sensitivity of the method to the absorption band in the 2–3 μm region is critical in the MMX/MIRS context as its main objective would be to confirm or refute the presence of hydrated material, and thus to constrain the origin of Phobos. Similarly, a good correction after 3 μm will clearly help the identification of bands due to aliphatic and organic materials.

Figure 9 displays the residual slope against the relative band depth after the correction from the thermal emission. It is clear that the location of the band has a strong influence on the quality of the correction, inducing a slope up to 7 × 10⁻³ μm⁻¹ for a deep band (50%) right in the thermal emission range (at 3.4 μm). However, the residual slope is smaller than 5% of Phobos spectral slope (2.99 × 10⁻¹ μm⁻¹ (Fraeman et al. 2012)) in the worst case, a deep band in the thermal emission range.

One can note that the relative band depth remains unchanged in the corrected spectra, and allows a direct analysis of the surface.

Comparison with the David et al. (2024) method. We compared our correction method to the iterative black-body fitting correction described in David et al. (2024). To do so, we used the same data without any absorption band presented in Section 3.1, and ran simultaneously the iterative method and the thermophysical model. In Fig. 10, we display the residual spectra from the thermophysical method (transparent plain lines) and the iterative black-body fitting method (dash-dotted lines). We found a discrepancy between the two methods: the iterative method induces an increasing slope on the complete spectral range. This slope more likely comes from the fit performed by the iterative fitting method to remove the reflective component of the spectra. This reflective component is not simply a straight line, as it depends on the surface reflective properties at each wavelength, specifically the SSA which is increasing on Phobos (Fraeman et al. 2014).

Moreover, for the warmest hours (12AM to 14PM), the spectra present a negative slope after 3 μm which is likely due to this first fit performed to remove the reflective component of the spectrum. An error on this reflective component would indeed tilt the black body, making it more or less steep, thus degrading the fit and inducing an absorption-like feature or a negative slope in the badly fitted areas.

An easy fix to this method would be to compute the reflectance of the surface (similarly to our method) in order to get a better estimation of the reflective component of the spectrum and avoid the apparition of features of this kind.

For the sensitivity test to the presence and depth of absorption feature in the spectra, we performed the same analysis as previously (Fig. 9) and we found a much larger residual slope (5× larger), which is due to fitting a straight line as the reflective component in the spectra, as already pointed out.

The iterative method of David et al. (2024) has a positive slope in the corrected I/F spectrum below 2.6 μm because the reflected component is modelled by fitting a linear function over the limited wavelength range of 2.6–2.9 μm. This linear continuum is then subtracted from the entire synthetic data spectrum to obtain the corrected I/F spectrum shown in Fig. 10a. As a result, the positive slope observed in the corrected I/F spectrum reflects the inherent non-linearity of Phobos’s SSA spectrum (Fraeman et al. 2014), rather than an error in the thermal emission estimation. The flat slopes of the corrected I/F above 2.6 μm (Fig. 10a) show that both methods can remove the thermal emission with adequate accuracy.

Based on these results, thermophysical modelling seems to be a more appropriate solution to remove the thermal emission from the surface as it does not induce any slope, takes into account several physical surface parameters, and is less sensitive to the presence of absorption features in the spectrum, specifically when they are located outside of the thermal emission spectral range. The main issue of this thermal model lies in the knowledge about the surface properties needed to get a good estimation of the surface temperature, and its computation time (several days for a global coverage). However, the first issue will progressively be fixed with the accumulation of data from both the IDEFIX rover (Michel et al. 2022) and the other instruments on board the MMX spacecraft.

Thus, we propose that both models can be used for the future MMX/MIRS data depending on the precision expected. For on-the-fly corrections, we recommend using the iterative fitting method (David et al. 2024) as it is much faster; instead, for precise estimations of band depths, and precise fitting and removal of the thermal emission, our thermophysical approach seems more appropriate.

Fig. 7 Synthetic data created with our model following the protocol in Section 3.1 for (a) 20 different local times and (b) eight different geometries of observation varying the phase angle ϕ, with the fitted spectrum for the correction. Panels (c) and (d) show the residual spectra respectively for synthetic data after subtracting the fitted thermal emission and the assumed reflected component without the artificial absorption band.

Fig. 8 Synthetic data created with our model following the protocol in Section 3.1 for (a) 20 different local times and (b) eight different geometries of observation, varying the phase angle ϕ, with the fitted spectrum for the correction. Panels (c) and (d) show respectively the residual spectra after removing the fitted thermal emission and the assumed reflected component with an artificial absorption band.

Fig. 9 Slope of the corrected spectra. The slope increases with the relative band depth. The location of the band also strongly changes the quality of the correction and, as expected, a band in the range 3–3.5 μm region has a strong impact on the correction.

Fig. 10 (Panel a) Residual spectra (with offsets) after correction with our method (plain transparent), and the iterative black-body fitting from David et al. (2024) (dash-dotted lines). The synthetic data used for the comparison are the same as that presented in Fig. 7 earlier. (Panel b) Evolution of the slope residual with the presence of an absorption band in the spectra.

Fig. 11 (Panel a) CRISM spectrum of Phobos and thermophysical simulation of the surface thermal emission fitted with grain size and porosity considered as free parameters. (Panel b) Corrected spectrum with our thermal model (fitted). The method show a very noisy result after 2.9 μm.

5.2.2 Removing the thermal emission in CRISM data

We tested both methods on real data with the spectral cube FRT00002992_03_IF162S_TRR7 observed on Phobos by CRISM (Murchie et al. 2007) on board the Mars Reconnaissance Orbiter (MRO, Zurek & Smrekar 2007) on October 23, 2007, around 1h18. Although the data are known to bear a calibration artefact around 3.18 μm (Murchie et al. 2009) in the zone 1–zone 2 filter boundary, we decided to try to remove the thermal emission in this data.

Figure 11 presents the CRISM and the simulation from our thermophysical simulation at the time of the observation (panel a). Overall, the simulation manages to reproduce quite well both the continuum [1–2.3 μm] and the thermal emission part [2.5–3.6] μm. However, we note a slight slope (2.5 × 10⁻² μm⁻¹) induced in the range [1–2.3 μ]m by our simulation after the correction in panel b. Since this spectral range bears the reflective component only, we suppose this slope is due to uncertainties in the Hapke parameters ($\overline{\theta }$, SSA, B₀, h_s) or due to averaging the spectra in the RoI.

The correction after 2.5 μm seems reliable on the global slope displaying the same slope before 2.5 μm; however, due to the very noisy data in this region, we cannot assess the presence of absorption features. This high noise is a consequence of a calibration issue (raised in Murchie et al. 2009) due to the limit of the filters on the IR sensor blocking the grating’s high orders of diffraction. Similarly, the spike at 2.75 μm is an artefact due to the filter’s limit (Murchie et al. 2009).

6 Conclusion

We developed a model, in support of the analysis of MIRS data by the MMX mission, dedicated to estimating the surface temperature of the Martian moons in order to remove the contribution of the surface thermal emission in the future MIRS observation. This model computes the flux absorbed by the surface of the Martian moons using the available shape models, and derives the sub-surface temperature profile over a complete seasonal cycle from knowledge of the surface textural and thermal properties.

Using this model, we showed that the Mars reflected and emitted light are increasing the temperature of about 7 K and have to be considered as it would significantly change the surface thermal emission in the NIR range.

Using this thermal model, we could simulate the flux emitted by the surface of Phobos and tried to remove this emission component from the spectral data measured by CRISM. The method has been proven to be less sensitive to the presence of absorption features in the spectrum than its iterative black-body fitting counterpart, and induce fewer artefacts. However, this approach presents the inconvenience of being strongly dependent on the surface properties, and takes much more computation time than the black-body fitting method. Thus, we propose this thermophysical model as a second approach of removing the surface thermal emission in the future MIRS observation in order to perform more detailed analyses of the surface composition.

In addition to the MMX mission, this model may be seen as a proxy to probe the surface and near sub-surface texture since the grain size and the porosity will strongly influence the thermal conductivity, and thus the heat transfer inside the surface. In that case, and to improve the fit of the black body, a larger spectral range would be needed.

Acknowledgements

This work was supported by the DIM ORIGIN – Ile de France (2023) and CNES.

References

All Figures

Fig. 1 Simple schematic of the model pipeline. The surface incoming flux is first computed in order to derive the surface temperature based on the ephemeris, the Mars albedo mosaic from Christensen et al. (2001), the shape model of the targeted moon, and its albedo. The self-heating is computed iteratively after a first run of the thermal model until the convergence is reached, reinjecting the previous estimation to compute the new surface temperature. The thermal model takes the grain size Dg, the surface porosity p and density ρ, and the non-porous thermal conductivity km and specific heat Cp as inputs. Eventually the surface thermal emission is computed for a rough surface knowing the single scattering albedo (SSA) ϖ; the surface Hapke’s average slope angle $\overline{\theta }$; the opposition effect amplitude B0 and angular width parameter hs; the phase function asymmetry parameter b; and the geometry of observation regrouping i, e, and ϕ for the incidence, observation, and phase angle, respectively. This thermal emission is to be used as a reference to remove the emission component in the spectra.

In the text

	Fig. 2 Simple drawing gathering all flux considered in the simulation. Solar flux (direct and reflected by Mars) is indicated in yellow. Thermal emission (Mars and self-heating) is represented in red.
In the text

Fig. 3 Simulation of the fluxes absorbed by the surface of Phobos on October 23, 2007. Stickney is the large crater on the right (310°E). This simulation shows the total absorbed flux on Phobos (direct and reflected solar flux, Mars thermal emission and self-illumination). The rectangle in black is the region of interest (RoI) where we plotted the absorbed flux over the complete revolution in Fig. 4. The black pentagon represents the position of the sub-solar point at the time of the simulation.

In the text

	Fig. 4 In the RoI (black square in Fig. 3), the absorbed flux reaches 600 W.m−2, while on Phobos the maximum is 713 W m−2. A minimum of 35 W m−2 is constantly illuminating the Mars sub-point due to the Mars thermal emission.
In the text

	Fig. 5 Temperature map at Phobos on 23 October 2007 at 1h17 UTC. The sub-solar point is located on the sub-Mars hemisphere, where the warmest temperatures are recorded (288 K with the given thermal and surface properties). The rectangle in black is the region of interest presented in Fig. 3. The black pentagon represent the position of the sub-solar point at the time of the simulation.
In the text

	Fig. 6 Thermal emission of a rough surface (computed from the Hapke mean slope parameter $\overline{\theta }=24$ for Phobos) showing the strong link with temperature of the thermal tail position.
In the text

Fig. 7 Synthetic data created with our model following the protocol in Section 3.1 for (a) 20 different local times and (b) eight different geometries of observation varying the phase angle ϕ, with the fitted spectrum for the correction. Panels (c) and (d) show the residual spectra respectively for synthetic data after subtracting the fitted thermal emission and the assumed reflected component without the artificial absorption band.

In the text

Fig. 8 Synthetic data created with our model following the protocol in Section 3.1 for (a) 20 different local times and (b) eight different geometries of observation, varying the phase angle ϕ, with the fitted spectrum for the correction. Panels (c) and (d) show respectively the residual spectra after removing the fitted thermal emission and the assumed reflected component with an artificial absorption band.

In the text

	Fig. 9 Slope of the corrected spectra. The slope increases with the relative band depth. The location of the band also strongly changes the quality of the correction and, as expected, a band in the range 3–3.5 μm region has a strong impact on the correction.
In the text

	Fig. 10 (Panel a) Residual spectra (with offsets) after correction with our method (plain transparent), and the iterative black-body fitting from David et al. (2024) (dash-dotted lines). The synthetic data used for the comparison are the same as that presented in Fig. 7 earlier. (Panel b) Evolution of the slope residual with the presence of an absorption band in the spectra.
In the text

	Fig. 11 (Panel a) CRISM spectrum of Phobos and thermophysical simulation of the surface thermal emission fitted with grain size and porosity considered as free parameters. (Panel b) Corrected spectrum with our thermal model (fitted). The method show a very noisy result after 2.9 μm.
In the text