Potential Surface Ice Distribution on Close-in Terrestrial Exoplanets around M dwarfs

Previous studies suggested that surface ice could be distributed on close-in terrestrial exoplanets around M-dwarfs if heat redistribution on the planets is very inefficient. In general, orbital and atmospheric parameters play an important role in the climate on terrestrial planets, including the cold-trap region where the permanent surface water reservoir can potentially be distributed. Here, we develop a simple coupled land-atmosphere model to explore the potential surface ice distribution on close-in terrestrial planets with various orbital and atmospheric parameters, assuming that the planets are airless or have a thin \ce{N2} atmosphere. We find that the most significant factors in deciding the surface cold trap region are the spin-orbit ratio and obliquity. The incident stellar flux and the surface pressure play a limited role in the thin \ce{N2} simulations for incident flux smaller than Mercury's and surface pressure lower than 10$^4$ Pa. Our result illustrates the possible distribution of surface ice on arid terrestrial planets and can help to understand the climate of these exoplanets.


INTRODUCTION
The successful launch of the James Webb Space Telescope (JWST) opens a new era for exploring the climate of terrestrial planets around M dwarfs (Morley et al. 2017;Batalha et al. 2018;Krissansen-Totton et al. 2018;Lustig-Yaeger et al. 2019;Gialluca et al. 2021;Wordsworth & Kreidberg 2022).Close-in terrestrial planets, such as TRAPPIST-1 b and c, have become key targets for atmospheric detection with JWST due to their close proximity to parent stars, and preliminary observations suggest that both TRAPPIST-1b and c have thin atmospheres or are airless (Greene et al. 2023;Zieba et al. 2023).
Close-in M dwarf planets are subject to strong tidal dissipation, whose rotation states are likely to be in spin-orbit resonance (Goldreich & Peale 1966;Heller et al. 2011;Barnes 2017), making their insolation pattern similar to Mercury rather than Earth.Previous studies suggest that surface water ice is allowed to accumulate in the cold surface region (usually referred to as the surface cold trap) even for planets receiving stellar flux that is several times higher than Earth's (Leconte et al. 2013;Ding & Wordsworth 2020, 2022), when the atmosphere has low infrared (IR) opacity and therefore cannot redistribute heat efficiently (Wordsworth 2015;Koll & Abbot 2016).Long-lived liquid water has also been proposed to possibly be present near the edge or at the bottom of the ice cap (Leconte et al. 2013;Lobo et al. 2023), raising an intriguing question regarding the regional habitability of such types of close-in exoplanets.
In this paper, our aim is to explore how the potential surface ice distribution on close-in terrestrial exoplanets is impacted by orbital and atmospheric parameters with a simple coupled land-atmosphere model.We briefly introduce our land-atmosphere model and the setup of the model in Section 2. We then apply it to airless M dwarf planets in Section 3 and to planets with thin N 2 atmospheres in Section 4. We present our conclusion in Section 5.

Coupled land-atmosphere model
We develop a simple coupled land-atmosphere model that calculates the evolution of atmospheric, surface, and subsurface temperatures over a full orbit with various orbital parameters, assuming that the atmosphere is transparent in IR with a low N 2 content.We assume that the surface of the planet is made of rocks, with thermal conductivity = 2.9 W m −1 K −1 ) and thermal diffusivity = 1.43× 10 −6 m 2 s −1 (Pierrehumbert 2010, Chapter 7.4).We assume that the simulated terrestrial planet has the same mass as Earth with a gravity value = 9.8 m s −2 and a  uniform surface albedo = 0.2.When there is no atmosphere at all, our model is very similar to the thermal model used to study the thermal phase curve of super-Mercury-type planets in Selsis et al. (2013).
We first test our model with the airless scenario and select six input parameters: star type, average stellar flux ( ), eccentricity (e), spin-orbit ratio ( = / ), obliquity ( ), and internal heat flux ( ) , where and are the orbit and rotation periods of the planet, respectively.Specifically, we apply the luminosity and mass values for M1, M5, and M9 type stars and ensure that the incident stellar flux on the planets is the same as Mercury's or Earth's by varying the planet's semi-major axis (a).In simulations with thin N 2 atmospheres, we investigate two additional variables: surface pressure ( 0 ) and average surface wind speed (U).
To better compare and estimate the effects of each parameter on the surface temperature distributions, we choose a reference simulation: tidally locked planet around a M5 type star with Earth-like incident stellar flux, circular orbit with 1:1 spin-orbit ratio, uniform internal heat flux of 90 mW m −2 , zero obliquity and precession angle.For the thin N 2 atmosphere scenario, the reference surface pressure is 10 4 Pa and the surface wind speed is 3 m s −1 .For each numerical experiment, we change one parameter based on the reference simulation and calculate the temporal evolution of the surface temperature distribution.Note that some of the input parameters in our simulations are not independent.For example, for eccentric orbits, the spin-orbit ratio cannot be 1:1 when the planet's spin becomes pseudo-synchronised and aligned with the orbit (Correia & Laskar 2010;Makarov 2012;Penev 2024).
Then we estimate the maximum coverage of permanent surface ice based on surface cold trap distributions in our simulations with the criterion that the surface ice deposit could stay for 10 Gyr against hydrogen escape (Section 2.5).All input parameters used in our simulations are listed in Table 1.

Incident stellar flux
We assume that the simulated planets revolve around their parent stars with Keplerian orbits.The incident stellar flux at any point on the surface of the planets can be calculated by using Kepler's law and spherical geometry (Pierrehumbert 2010, Chapter 7.5).For a planet orbiting around a star with the mass of M and the semi-major and semi-minor axis of a, b, its period and the angular momentum of the revolution are Because the spin-orbit ratio is given as input parameters in our simulation, the period and angular velocity of rotation , can be calculated as = 2 / = / .Suppose that the planet is at the northern hemisphere solstice when = 0, and is the season angle that describes its position relative to the initial one.The evolution of the season angle ( ) can be derived by the conservation of angular momentum as Then, the longitude and latitude of the substellar point on the surface of the simulated planet ( 0 , 0 ) can be calculated by where , is the precession angle and obliquity.The zenith angle at any point on the surface of the planet ( , ) can be calculated by cos ( ) = cos cos 0 ( ) cos [ − 0 ( )] + sin sin 0 ( ).
For cos ( ) < 0, the star is below the horizon and no incident stellar flux can reach the surface.Let the star's luminosity be L.The incident stellar flux at any point on the surface of the planet is (5) In our simulations, the integration time step is Δ = /200.The horizontal grid resolution is 18 • × 18 • .We integrate our model for 50 planetary years to ensure that the simulated climate reaches equilibrium.

Heat conduction in subsurface layer
As the vertical temperature gradients in the subsurface layer are typically much greater than the horizontal one, only vertical heat conduction is taken into account in our simulations.The energy equation in the subsurface layer is Here , is the density and the specific heat capacity of the subsurface layer.The lower boundary condition is that the diffusive heat flux there is balanced by the internal heat flux .For the boundary condition at the surface in the airless scenario, the diffusive heat flux at the surface is balanced by the net heating by insolation and thermal emission The infrared emissivity of the surface is assumed to be unity and ( ) = 4 where is the surface temperature.We set 35 layers in the subsurface layer and the depth of each layer i is discretized as where = Δ / , which is the characteristic depth to which temperature fluctuation can penetrate within one time step.The maximum depth 35 in our simulations with various orbital parameters is within 10 metres.

Reproducing the surface temperature evolution of the Moon and Mercury
To validate our model in the airless scenario, we use our model to simulate the surface temperature evolution on the Moon and Mercury.For the soil heat conduction equation, the lunar regolith thermal parameters ( = 0.01 W m −1 K −1 , = 10 −8 m 2 s −1 ) are used for the Moon and Mercury (Pierrehumbert 2010, Chapter 7.4).The surface albedo is set to be globally uniform, with 0.136 for the Moon and 0.088 for Mercury 1 , and the horizontal grid resolution is 18 • × 18 • .The simulated surface temperature evolution at the equators of the Moon and Mercury (Fig. 1) is very similar to previous work with similar thermal properties of the regolith (Vasavada et al. 1999).The Moon's surface temperature at the equator varies from ∼ 390 K at noon to ∼ 100 K at night.On Mercury it varies from ∼ 570 K at noon to ∼ 100 K at night.The secondary sunrise and sunset at specific longitudes on Mercury, due to the faster orbital angular velocity than its spin rate at perihelion, are also correctly captured in our simulation.The surface temperature variation in our simulations may be slightly different from the observed values, because the thermophysical properties of the regolith are depth-and temperature-dependent (Vasavada et al. 1999;Hayne et al. 2017).

Estimation of surface cold-trap distribution
We use a critical temperature to estimate the permanent surface cold-trap distribution assuming that an ice layer with thickness of 1 km could stay for 10 Gyr against hydrogen escape.For the airless scenario, we simply assume that all the water molecules sublimated from the ice surface can leave the planet immediately by photodissociation into hydrogen molecules and subsequent hydrogen escape.Thus, the loss rate of the surface ice sheet depends only on the sublimation rate where sticking coefficient = 1, = / H 2 O is the density and = ( / H 2 O ) 1/2 is the molecular speed of water vapour in equilibrium with the surface ice, the background atmosphere pressure 0 = 0 in the airless scenario and the water vapour pressure is equivalent to the saturation vapor pressure = ( ) (Ingersoll et al. 1985).We use the empirical formula (Huang 2018) to estimate the saturation vapor pressure in phase equilibrium with the surface ice with temperature in Celsius In this way, the loss rate of the surface ice sheet can be expressed as a function of surface temperature.For a conservative estimate, we assume that the ice sheet can last for more than = 20 Gyr, attaining the temperature threshold in the airless scenario with 0 ≈ 140 K by ( 0 ) ≤ / .
Here is the mixing ratio of H 2 and is the temperature in the upper atmosphere, = 1.9 × 10 21 ( /300 ) 0.75 m −1 s −1 is the binary diffusion parameter for H 2 in the N 2 atmosphere, , H 2 is the molecular mass of N 2 and H 2 .The mixing ratio of hydrogen in the atmosphere can be estimated approximately as the saturation value of surface water vapour, because the potential ice layer cannot be stably trapped on the surface with an effective tropopause cold trap (Ding & Wordsworth 2020).
is estimated by the skin temperature of the planet.We also assume efficient horizontal mixing of water vapour in the atmosphere.Therefore, the loss rate of the surface ice sheet is relevant to the area ratio of the surface ice sheet.Then the lifetime of surface ice sheet can be calculated as where hydrogen escape is assumed to happen on a global scale with = 42 .The area ratio of surface ice / can be determined by the cold trap region with a surface temperature lower than the critical temperature over a full orbit.Given that ≥ 10 Gyr, the critical temperature can be solved by the surface temperature distribution in simulations with various parameters.

Energy conservation with a thin nitrogen atmosphere
For an IR-transparent atmosphere, the atmosphere can only exchange energy with the underlying surface by the turbulent sensible heat flux (Pierrehumbert 2011), and the atmosphere can still transport heat horizontally by air flow.Due to the long non-dynamical (e.g., radiative) timescale in the pure N 2 atmosphere relative to the short dynamical adjustment timescale, the weak-temperature-gradient (WTG) approximation can be applied to the atmosphere.Using the WTG approximation to implicitly solve for horizontal heat transport in IR-transparent atmospheres on synchronously rotating planets was first proposed in Pierrehumbert (2011), and verified by recent global climate modelling by Wang & Yang (2022).With the WTG approximation, the near-surface atmosphere has a horizontally uniform temperature and the energy conservation equation for the atmosphere is 0 = ( − ). ( Here 0 is the surface N 2 pressure, is the specific heat capacity of N 2 , Parameter a is the turbulent coupling coefficient calculated as , where = 0 / N 2 is the surface air density, is the average wind speed, and is a dimensionless drag coefficient, with a typical value of 0.0015 over moderately rough surface (Garratt 1994).
The sensible heat flux also changes the upper boundary condition of the heat conduction in the subsurface layer compared to that of the airless case Compared to the surface energy budget of typical terrestrial climates, the thermal emission flux of the atmosphere is ignored in Eq. 14 because the atmosphere is assumed to be transparent in IR.However, the existence of sublimated water vapour from the potential ice surface may violate this assumption.To evaluate the greenhouse effect contributed by water vapour, a line-by-line radiation code (PyRADS-shortwave) is used 2 , which has been used to study the greenhouse effect of water vapour on Earth and other exoplanets (Koll & Cronin 2018, 2019).Assuming a background surface pressure of 10 4 Pa, a vertical uniform temperature profile of 200 K and water vapour profile with the saturation value at the surface, we find that the downward infrared flux reaching the surface is less than 6 W m −2 , which is one order of magnitude lower than the maximum sensible heat flux in our reference simulation.This comparison confirms that our assumption with zero IR opacity can be treated as a first-order approximation, although the simulated surface temperature by Eq. 14 may be slightly underestimated if the greenhouse effect of sublimated water vapour is taken into account.
We also use PyRADS-shortwave to evaluate the shortwave radiative effect of the 10 4 Pa N 2 atmosphere with low H 2 O concentrations, because H 2 O has plenty of absorption lines in the near-infrared and can potentially decrease the planetary albedo while the Rayleigh scattering by N 2 can slightly increase the planetary albedo.For a late M-dwarf with an effective temperature of 2500 K, 1.2% of incoming stellar flux is absorbed in the atmosphere by H 2 O and the planetary albedo is reduced by 0.3% relative to the surface albedo.For an early M-dwarf with an effective temperature of 3700 K, only 0.65% of incoming stellar flux is absorbed in the atmosphere by H 2 O and the planetary albedo is reduced by 0.05%.The radiative calculations suggest that the stellar energy distribution has a trivial effect on the energy budget of the terrestrial planets with thin atmospheres, and validate our assumption that ignores the shortwave radiative effect of the thin atmospheres in this work.

AIRLESS SCENARIO
For airless planets, the surface can be warmed by either the incident stellar flux from above or the diffusive heat flux from below and cooled by thermal emission to space (Eq.7).The surface temperature distribution in the reference simulation on airless planets is shown in Fig. 2. As the planet is in the synchronous rotation state with zero obliquity and eccentricity, the incident stellar flux stays constant at any point on the surface, resulting in a time-invariant equilibrium state.The hottest place coincides with the substellar point, reaching approximately 360 K as a result of intense incoming shortwave radiation.The coldest place occurs on the nightside, where the temperature drops below 100 K and is solely heated by the planet's internal heat flux.According to our estimate, the critical temperature to allow the existence of surface ice for 10 Gyr is about 140 K, indicating that permanent surface ice could potentially span the entire nightside in the reference simulation.
Fig. 3 shows the maximum coverage of permanent surface ice under different parameters listed in Table 1.In the reference state with a permanent nightside, the surface ice is allowed to spread throughout the night hemisphere and extend marginally to the dayside near the terminator where the insolation is very low.Therefore, the surface ice coverage is slightly higher than half.Among all the parameters we survey, only the influence of spin-orbit ratio, eccentricity, and obliquity is not negligible, while the star type, average stellar flux, and internal heat flux have a trivial effect.Non-1:1 spin-orbit ratio with non-zero eccentricities or non-zero obliquities can expose more surface to incident stellar radiation and warm the exposed surface.Note that the permanent surface ice distribution in our estimate requires that the surface temperature stays below 140 K over the course of a full orbit.
For planets with a non-1:1 spin-orbit ratio, the substellar point can move along the zonal direction compared to the synchronous rotation state, which forces the permanent surface cold trap to retreat to the polar region.This is similar to Mercury's state with the presence of water ice near the north pole despite the fact that this planet receives an averaged stellar flux seven times as high as the Earth (Lawrence et al. 2013).
Non-zero obliquity also plays an important role in the evolution of the substellar point on the planet, leading to the motion of the substellar point along the meridional direction and the seasonal variation in both the northern and southern hemispheres, similar to the seasonal cycles on Earth or Mars.When obliquity increases, more surface area near the poles will be exposed to the incident stellar flux in the summer season even with the 1:1 spin-orbit ratio, causing the permanent surface cold trap to retreat to mid and low latitudes on the nightside.
The incident stellar flux at the substellar point and the internal heat flux can change the surface temperature distribution, as energy inputs from the top and the bottom of the surface (Section 2.3), respectively.However, the surface temperature of the nightside in our reference simulation is only relevant to the internal heat flux, which is not sufficient to warm the nightside to the critical temperature.The surface temperature near the terminator is also altered marginally by these two parameters.Thus, the incident stellar flux and the internal heat flux have a very trivial impact on the distribution of surface cold traps in our simulations.The star type has no effect on the surface temperature distribution in Fig. 3, because the incident stellar flux, spin-orbit ratio, and surface albedo of the planet are fixed as in our reference simulation, although the actual orbital period and rotation period can change.

SCENARIO WITH THIN NITROGEN ATMOSPHERES
For planets with a thin N 2 atmosphere, the surface temperature is subject to sensible heat transfer between the surface and the atmosphere, in addition to incident stellar flux, diffusive heat flux, and thermal emission flux in the airless scenario (Eq.14).The surface temperature distribution in the reference simulation with 10 4 Pa N 2 and an average surface wind speed of 3 m s −1 is shown in Fig. 4. Compared to the airless case in Fig. 2, the highest temperature at the substellar point decreases, while the surface temperature on the nightside increases.This more uniform surface temperature distribution results from sensible heat transfer between the surface and the atmosphere, which cools the surface region hotter than the atmosphere and warms the surface region colder than the atmosphere.This heat redistribution mechanism is different from that dominated by radiative transfer in the thermal IR described in Koll & Abbot (2016).Furthermore, the critical temperature to allow the existence of surface ice for 10 Gyr, indicated by the pink contour in Fig. 4 , notably increases to 220 K in our reference simulation, because nitrogen atmospheres effectively slow hydrogen escape and thus reduce the loss rate of the surface ice sheet.These combined factors drive the expansion of surface ice coverage to ∼60% under Earth-like insolation, further into the dayside hemisphere compared to the airless case.
Fig. 5 shows the maximum coverage of permanent surface ice for planets with thin N 2 atmospheres.Similar to the airless case, the surface ice coverage reaches its maximum in the reference simulation with 1:1 spin-orbit ratio and zero obliquity.The spin-orbit ratio, obliquity, and eccentricity play a significant role in the surface temperature distributions.Unlike airless simulations, incident stellar flux can affect surface cold trap coverage in a minor way.When the planet receives incident flux as high as Mercury, the intense insolation warms the terminator region significantly while still keeping the nightside cold.It causes the surface cold trap region to retreat to the nightside, similar to the airless case.The effects of atmospheric parameters also appear to be relatively minor, although both surface wind speed and surface air pressure can affect sensible heat transfer between the surface and the atmosphere.For the typical values that we explore in the simulations, their impact is smaller than the orbital parameters we discussed above, even when the surface pressure is varied by two orders of magnitude.

CONCLUSION AND IMPLICATION FOR PLANET SURVEY
We investigate the potential surface ice distribution of close-in M dwarf planets using a simple coupled land-atmosphere model and discuss the effects of orbital and atmospheric parameters in airless and thin N 2 atmosphere scenarios.In general, the 1:1 spin-orbit ratio and zero obliquity favour a global-scale distribution of surface cold traps with possible surface ice.For terrestrial planets that receive incident stellar flux less than that on Mercury and have a N 2 atmosphere thinner than 10 4 Pa, stronger incident stellar fluxes and thinner atmospheres tend to decrease the surface coverage of cold traps, but in a limited way.
Based on our simulation results, we filter the confirmed exoplanets around M dwarfs3 and find several candidates that can potentially have a large-scale surface ice distribution on its nightside.We use the criterion with ≤ 1.5 ⊕ and ≤ 5 ⊕ to ensure that it is a terrestrial planet, eccentricity < 0.05 and the incident stellar flux between that of Mercury and Earth.The filtered result shows only several planets such as Proxima Centauri d and TRAPPIST-1 b-d can be hopeful candidates, and other close-in terrestrial planets either have unknown or large eccentricities, or are too far away from the Earth to be characterised.
In this work, we only estimate the potential surface ice coverage by computing the surface cold trap distributions.Close-in planets may be water poor if they formed during the prolonged pre-main sequence phase of their parent stars (Ramirez & Kaltenegger 2014;Luger & Barnes 2015;Tian & Ida 2015).Without a water supply for the nightside surface, the area can still be completely dry.Recent work discussed the possibility of secondary atmosphere buildup on close-in terrestrial planets by volcanic outgassing from the hydrated interior (Kite & Barnett 2020).Other outgassed molecules such as CO 2 , CO and CH 4 can increase the IR opacity of the atmosphere and thus lower the chance of surface ice formation.Mineral dust is another complicating factor.Dust particles can be lifted to the atmosphere by dust activity on the dayside and transported to the nightside by large-scale circulation, where dust can increase the IR opacity and warm the nightside surface similar to greenhouse gases (Boutle et al. 2020).However, recent GCM simulations suggested that the convective behaviour in the atmosphere on tidally locked planets, which is crucial to trigger dust activity, still requires a substantial background IR opacity (Wang & Yang 2022).At last, it should be noted that other mechanisms have been proposed to increase the water content on terrestrial planets around M dwarfs against atmospheric escapes, e.g., migration of planets in the proto-planetary disks (Alibert & Benz 2017) and oxidation of atmospheric hydrogen by rocky materials from incoming planetesimals and from the magma ocean (Kimura & Ikoma 2022).
Another simplification we made in this work is that the incident radiation from the host star reaches the top of the atmosphere on the terrestrial planet as plane-parallel rays, which works well for terrestrial planets in the solar system.However, for close-in terrestrial planets around their host stars, the finite angular size of stars may make more than 50% of the planet receive the stellar radiation.This hyper-illumination effect is important for the climate dynamics of lava planets that receive intense stellar radiation because of their close proximity to the host stars (Nguyen et al. 2020;Kang et al. 2023).Recently, Carter et al. (2024) showed that even for a temperate terrestrial planet such as TRAPPIST-1 d, the illumination can reach 51.2% and is still slightly higher than one-half.The hyper-illumination effect may make the surface ice distribution on close-in terrestrial planets around M dwarfs more sensitive to the orbital obliquity than discussed in our present work.These factors discussed above should be taken into account in future studies to better understand the climate on close-in terrestrial planets.

Figure 1 .
Figure 1.Simulated surface temperature at the equator of Mercury (top) and Moon (bottom) as a function of local time by our simple land-atmosphere model.The Mercury results are for 90 • E longitude, where secondary sunrise and sunset happen due to the faster orbital angular velocity than its spin rate at perihelion.

Figure 2 .
Figure 2. Surface temperature distribution in the reference simulation on airless planets (unit: Kelvin).The yellow contour marks the critical temperature of ∼140 K permitting the existence of permanent surface ice that spans the entire nightside.

Figure 3 .
Figure3.Sensitivity of the maximum coverage of permanent surface ice to the input parameters defined in Section 2.1.The dash horizontal line marks the maximum coverage of permanent surface ice in the reference simulation.The red vertical lines mark the range of ice coverage when the parameter listed in x-axis varies.Note that the spin-orbit ratio and orbital eccentricity are not independent parameters and thus are taken into account together in the horizontal axis.

Figure 4 .Figure 5 .
Figure 4. Surface temperature distribution in the reference simulation on planets with a thin N 2 atmosphere (unit: Kelvin).The pink contour marks the critical temperature of ∼220 K permitting the existence of permanent surface ice.

Table 1 .
Parameters used in our simulations.Stars mark the values for the reference simulation on airless planets, and daggers for planets with thin N 2 atmospheres.