Multiple moist climate equilibrium states on arid rocky M-dwarf planets: A last-saturation tracer analysis

Terrestrial-type exoplanets orbiting nearby red dwarf stars (M dwarfs) are the first potentially habitable exoplanets suitable for atmospheric characterization in the near future. Understanding the stability of water in cold-trap regions on such planets is critical because it directly impacts transmission spectroscopy observations, the global energy budget, and long-term surface water evolution. Here we diagnose the humidity distribution in idealized general circulation model (GCM) simulations of terrestrial-type exoplanets. We use the `tracer of last saturation' technique to study the saturation statistics of air parcels. We find that on synchronously rotating planets, the water vapor abundance in the nightside upper troposphere depends weakly on planetary rotation, while more water vapor builds up in the nightside lower troposphere on fast-rotating planets. We then discuss how last-saturation statistics can elucidate the multiple moist climate equilibrium states on synchronously and asynchronously rotating arid planets. We show that the multiple moist climate states arise from the cold-trapping competition between the substellar upper atmosphere and cold surface regions. We find that fast synchronously rotating planets tend to trap surface water on the nightside as a result of their weak atmospheric and strong surface cold traps compared to the slow rotating case. These results elucidate the nature of the water cycle on arid rocky exoplanets and will aid interpretation of atmospheric observations in the future.


INTRODUCTION
As hundreds of Earth-sized planets are discovered beyond our solar system 1 , the new frontier of exoplanet science will soon be characterization of atmospheres on these terrestrial planets. For the next few years, the most observable of these planets will be those orbiting in the habitable zone around nearby M-dwarfs. Potentially habitable exoplanets around M-dwarfs may be very different from Earth because they are close to their host stars and expected to be tidally resonant or locked due to strong gravitational tidal interactions. In the most extreme case, close-in planets in circular orbits will be synchronously rotating, with a permanent dayside and nightside (Kasting et al. 1993;Barnes 2017).
M-dwarf stars have a prolonged pre-main-sequence phase (Baraffe et al. 2015), which has been predicted to cause extensive early water loss on any rocky exoplanets orbiting in their habitable zones (Ramirez & Kaltenegger 2014;Luger & Barnes 2015;Tian & Ida 2015). As a result, many such planets may be water-poor. Understanding the potential atmospheric water content of such planets is vital for several reasons. First, for transit observations, water vapor profiles near the terminator directly impact the transmission spectroscopy of water vapor. Second, water vapor is critical to the energy budget of the climate system: for example, sub-saturated atmospheric regions can behave like "radiator fins" and stabilize the climate (Yang & Abbot 2014;Ding & Pierrehumbert atmospheric water vapor concentration above the surface water reservoir determines whether surface condensation or evaporation should occur, which controls the long-term evolution of the surface water inventory. For example, Ding & Wordsworth (2020) showed that on arid tidally locked planets covered by ice caps on the nightside surface, the difference between nightside atmospheric water vapor concentration and the saturation value at the surface determines whether the surface water will experience sublimation or condensation and hence whether the remaining water inventory will migrate towards the substellar region, forming an 'oasis', or still be stably trapped on the nightside surface. While Ding & Wordsworth (2020) only focused on slowly rotating planets, here we conduct GCM simulations to explore multiple moist climate equilibrium states in a broader parameter space.
Explaining the emergence of multiple moist climate systems requires a full understanding of water vapor distribution in the atmosphere. In a moist climate, the primary source and sink for atmospheric water vapor are evaporation of surface water and precipitation in the atmosphere, respectively. Beyond this, the atmospheric humidity distribution is strongly influenced by atmospheric dynamics. An idealized model of non-local large-scale control of tropospheric humidity, called the last-saturation model or advection-condensation model, has been proposed for understanding the humidity distribution in Earth's atmosphere (Pierrehumbert 1998;Pierrehumbert et al. 2007;Pierrehumbert & Ding 2016). This model assumes the mixing ratio of water vapor in an air parcel is conserved along any segment of a trajectory unless it becomes saturated and loses water by precipitation. In general, there are two ways of describing advective-condensation processes: the Lagrangian approach and the Eulerian approach. For the Lagrangian approach, Pierrehumbert (1998) used the backward Lagrangian trajectory technique to reconstruct relative humidity in the Earth's subtropics and showed close agreement with Meteosat satellite observations. Yang et al. (2019) also used the same method to interpret the differences in nightside humidity distributions on tidally locked planets among some GCM simulations. The back trajectory technique can reconstruct the humidity distribution by the three-dimensional (3D) wind and temperature fields with a high spatial resolution, but it only provides the humidity distribution as a single snapshot. For climate systems with strong transient motions, especially on tidally locked planets (Merlis & Schneider 2010;Pierrehumbert & Hammond 2019), this technique is not optimal for presenting last-saturation statistics in the equilibrium climate system averaged over a long time period.
Here we instead take an Eulerian approach to describe the last-saturation statistics in GCM simulations with tracers of last saturation for the first time. These tracers of last-saturation are advected by the air flow in the same way as water vapor in the GCM, and their concentration intrinsically contain information regarding water vapor's last-saturation statistics over a long time period. Galewsky et al. (2005) first developed this method to study the dryness of the Earth's subtropical troposphere in both idealized GCM simulations and a reanalysis dataset. They found that the subsaturation is primarily controlled by the isentropic transport by midlatitude eddies and that diabatic descending transport from the tropical upper troposphere plays a secondary role. However, this method has not been applied to tidally locked planets so far.
We briefly introduce this method and our GCM setup in Section 2. We show the multiple moist climate equilibrium states on slowly-rotating and rapidly-rotating tidally locked arid planets in Section 3. We diagnose the humidity distributions on slowly-rotating and rapidly-rotating tidally locked planets in Section 4. We then discuss how the last-saturation statistics elucidates the multiple moist climate equilibrium states on arid planets in Section 5 and present our conclusions in Section 6.

Idealized moist general circulation model
Previously, to study the hydrological cycles on arid planets, we developed a moist GCM that calculates the radiative transfer in the line-by-line approach and uses the Manabe moist convective adjustment scheme in Ding & Wordsworth (2020). Here, to study a wider range of climate states in an efficient way, we choose instead to simplify the radiation and convection calculations and make the moist GCM more idealized.

Two-stream gray radiation scheme
Here we use a gray-gas radiation scheme, which is a standard approach for exploring climate states over a wide range for Earth-like atmospheres (O'Gorman & Schneider 2008) and for atmospheres on terrestrial exoplanets (Merlis & Schneider 2010;Koll & Abbot 2015. Following Merlis & Schneider (2010), the longwave optical depth τ in our model has two contributing factors Here κ wv and q are the absorption cross-section and the mass concentration of water vapor, g is surface gravity, p 0 is the global mean surface pressure and τ dry is the total optical thickness of the non-condensible absorber, which is assumed to be well mixed in the atmosphere. In Merlis & Schneider (2010), κ wv = 0.1 m 2 kg −1 and τ dry = 1.2. Our model uses the same value for κ wv . Note that in this approach, the angle of propagation is implicitly incorporated into Eq. (1).

Convective parameterization
Following Galewsky et al. (2005), our idealized GCM only considers the large-scale condensation of water vapor. In any super-saturated grid box, the excess moisture is removed immediately. Frierson (2007) studied the effect of various convection schemes on Earth's zonally averaged climate: the largescale condensation (LSC) only scheme used here, the Manabe moist convective adjustment (MCA) scheme used in Ding & Wordsworth (2020), and the simplified Betts-Miller (SBM) scheme used in Merlis & Schneider (2010). The global surface temperature distributions are similar in simulations with those convection schemes, while tropical precipitation is less in the SBM simulation than in the LSC and MCA simulations. Use of the simple LSC convection scheme allows us to focus entirely on the transport role of large-scale dynamics, which is the main aim of this study.

GCM setup for arid planet simulations
The other aspects of our model setup are similar, but not identical, to those used in Merlis & Schneider (2010). The planet has the same radius and surface gravity as Earth's and both the eccentricity and obliquity are zero. The surface of the planet is flat and the incoming stellar radiation above the substellar point is 1367 W m −2 . The atmosphere is made of N 2 and condensible water vapor.
We assume the simulated planets are in synchronous rotation, which implies the rotation period equals the orbital period. We calculate the orbital period P self-consistently using Kepler's third law (Wordsworth 2015;Kopparapu et al. 2016;Haqq-Misra et al. 2018) as where L/L is the luminosity of the host star scaled by the luminosity of the Sun, and M/M is the mass of the host star in solar mass unit. Two planetary rotation periods are investigated here: 50 Earth days for a synchronously rotating terrestrial planet orbiting in the habitable zone around a M1V type star, and 10 Earth days around a M5V type star.
To explore the multiple moist climate equilibrium states on arid planets first studied in Ding & Wordsworth (2020), we perform an additional set of simulations with our idealized moist GCM. For each simulation, we first place the initial surface water inventory within 40 • around the substellar point, and run the GCM for 1000 days. If no condensation occurs on the nightside surface, we refer to this climate state as the substellar oasis state. Next, we place the initial surface water inventory on the nightside, and run the GCM with the same parameters for 1000 days. If no condensation occurs at the substellar tropopause, we refer to this climate state as the nightside icecap state. While Ding & Wordsworth (2020) only reported GCM simulations on a slowly rotating planet with two different CO 2 concentrations, here we explore a broader parameter space using the idealized moist GCM. We run these simulations with various surface pressures (p 0 in Eq. 1), optical thicknesses of non-condensible components (τ dry in Eq. 1) and planetary rotation rates to verify the surface water state (parameters listed in Table 1).

Last saturation tracer calculation
For the last saturation tracer calculations, the global model domain is divided into N axisymmetric subdomains along the axis from the substellar point to the antistellar point, D i (i = 1, . . . , N ). In other words, these N subdomains are zonally symmetric in the tidally locked coordinate in which Each subdomain D i is associated with a unique tracer of last saturation, T i . If a grid box reaches water vapor saturation in one subdomain, D i , we set the concentration of tracer associated with D i to unity and all the remaining tracers to zero in that saturated grid box, This is performed whenever a grid box reaches saturation in the GCM. Otherwise, all the tracers are simply advected by the large-scale circulation, in the same way as water vapor.
In general, an unsaturated grid box in the free atmosphere above the planetary boundary layer consists of air parcels last saturated in multiple subdomains. Since each subdomain is associated with its own tracer, the value of an individual tracer T i in that grid box can be interpreted as the probability that the air parcel was last saturated in its corresponding subdomain. Then the time mean and zonal mean water vapor distribution in the tidally locked coordinate can be approximately reconstructed by the time and zonal mean probability of last saturation in each subdomain D i and the saturation concentrations of water vapor q sat,i determined by the average temperature in the respective subdomains, where r is the position of the grid box. The two terms in the product of Eq. 4, the concentration of last saturation tracer and the saturation concentration of water vapor, represent the effect of atmospheric circulation (i.e., advection) and condensation on the water vapor buildup in the grid box.

GCM setup for last-saturation tracer analysis
The main purpose of our analysis with tracers of last saturation is to investigate the key physical processes for water vapor buildup on tidally locked planets, and how the humidity distribution is affected by planetary rotation. For this part of the analysis, we therefore assume a global ocean surface or aqua-planet, which is a common approach in GCM simulations of habitable tidally locked planets ( First, we present our moist climate state results in the arid case as a function of rotation rate, surface pressure and atmospheric optical thickness. Figure 1 shows the states of surface water inventory on arid synchronously rotating planets under surface pressures ranging from 0.5 bar to 5 bar and optical thicknesses of non-condensible components ranging from 0 to 3. The surface water on slow rotators is trapped on the nightside when the atmospheric greenhouse effect is weak and in the substellar region when the greenhouse effect is strong, which is consistent with GCM simulation results with realistic radiative transfer calculation in Ding & Wordsworth (2020). In addition, Figure 1 indicates that a thicker atmosphere tends to trap surface water on the dayside more effectively on slow rotators, because stronger upwelling motion at the the substellar tropopause strengthens the cold-trapping effect there. Comparing Figure 1a and b indicates that planetary rotation has a big impact on the moist climate states. Specifically, under the same conditions, fast rotating planets trap the surface water inventory on the nightside as ice more effectively than slow rotating planets. In the parameter space we explore here, fast rotators have no substellar oasis state.
In addition to the substellar oasis and nightside icecap states discussed in Ding & Wordsworth (2020), we find another climate state when the surface ice is stably trapped on the nightside but water vapor condensation occurs near the substellar tropopause. Here we refer to this state as the transition regime. Table 2 summarizes the key characteristics of the three surface water states.  To motivate the last saturation analysis further, the time mean and zonal mean water vapor concentration calculated from the GCM in the tidally locked coordinate for both the slow rotating and fast rotating aquaplanet simulations are shown in Figure 2. There are two important features in the nightside distribution of water vapor. First, in both simulations, the water vapor concentration on the nightside decreases with height. Second, comparison between the two distributions shows that more water vapor builds up in the nightside lower atmosphere (p > 700 hPa) in the fast rotating simulation, but water vapor concentrations are similar in the upper troposphere in both simulations.
As discussed in Section 1, the water vapor abundance in the nightside lower atmosphere is important for the behavior of the nightside "radiator fin" on synchronously rotating planets and the long-term surface evolution on arid planets with limited surface water inventories. In this section, we diagnose the water vapor distribution using tracers of last saturation to explain the above features.   But for water vapor in the lower atmosphere, other than the substellar tropopause, the air parcel has nearly the same chance of last being saturated in the region right above the core of the overturn-ing circulation (Figure 4b), where both air temperature and the saturation concentration of water vapor are higher than at the tropopause. As a result, on the nightside of slowly and synchronously rotating planets, there is a higher water vapor abundance in the lower atmosphere than in the upper atmosphere. In Eq. 4, the water vapor concentration at the reference point can be interpreted as the summation of the product of the probability of last saturation in each subdomain and the saturation concentration of water vapor in the respective subdomain. Because of this, we can reconstruct the water vapor concentration in the atmosphere only using tracers of last saturation and air temperature in each subdomain. The reconstructed water vapor field in Figure 5 is very similar to the field calculated by the GCM (Figure 2a), indicating our diagnostics with tracers of last saturation captures essential physical processes for nightside water vapor buildup on slowly and synchronously rotating planets.

Fast rotating simulation
In the fast rotating simulation with planetary rotation period of 10 days, the atmospheric circulation is not zonally symmetric in the tidally locked coordinate. Although the meridional circulation is still dominated by the overturning circulation, perturbations by stationary planetary waves are nonnegligible, which was discussed in detail by applying the Helmholtz decomposition on atmospheric circulation in Hammond & Lewis (2021).   (Figure 6b). Meanwhile, two surface cold spots appear westwards of the west terminator associated with anticyclone-type circulations. Figure 6d shows super-rotating equatorial jet and mid-latitude stationary planetary wave patterns with the zonal wavenumber 1 in the lower atmosphere of the fast rotating simulation. The interactions between the mean flow and planetary waves and the shifts of the temperature extrema were discussed in both a linear shallow water model and an idealized dry GCM by Hammond & Pierrehumbert (2018). Hammond & Pierrehumbert (2018) showed during the spin-up stage of their GCM that the planetary wave pattern is Doppler-shifted eastwards to the equilibrium position as the zonal jet forms. The parameters of the default run in Hammond & Pierrehumbert (2018) are identical to those of our fast rotating simulation (i.e., an Earth-sized planet with a 1 bar atmosphere dominated by N 2 , the same incoming stellar radiation as on Earth and a 10 day orbital period) except that their atmosphere is dry with a longwave optical depth of 1. Their discussions on global circulation therefore qualitatively apply to our fast rotating simulation. The cold-trapping statistics of water vapor are greatly influenced by this change of atmospheric circulation compared to the slow rotating simulation. This becomes clear when we analyze the probability distributions for locations of last saturation at reference points in the nightside troposphere ( Figure 7a and b). In Figure 7, the reference points are the same as in the slow rotating simulation (Figure 4a and b). For the upper troposphere on the nightside, water vapor was still last saturated in the tropopause region (Figure 7a), but shifted eastwards from the substellar point associated the eastward shift of surface hot spot in Figure 6b. The tropopause temperatures are very close in the slow and fast rotating simulations, so that the water vapor concentration in the nightside upper troposphere are similar (Figure 2). For the lower troposphere on the nightside, the last saturation statistics becomes more distributed compared to the slow rotating case. Water vapor in the lower troposphere was still last saturated in the dayside upper atmosphere above the core of overturning cell, but the probability is only one third of that in the slow rotating simulation, as can be seen by comparing Figure 4b and Figure 7b. Figure 7b shows that the water vapor in the lower troposphere was twice as likely to be last saturated in the middle to lower troposphere near the terminator, where the air is much warmer than in the upper troposphere in both simulations (Figure 4c and Figure 7c).

Sub-stellar Terminator
Therefore, an excess of water vapor builds up in the nightside lower atmosphere on fast rotating planets. towards the nightside directly. It resembles the 'atmospheric river' that brings sustained and heavy precipitation to the west coasts of North America in the Earth's climate system (Zhu & Newell 1994, 1998 Noda et al. (2017). The circulation with the insolation pattern relatively in phase with the deep convection region but out of phase with mid-latitude stationary waves of zonal wavenumber 1 is analogous to the 'Rhines' regime defined in Haqq-Misra et al. (2018) and the 'Type-II' circulation regime defined in Noda et al. (2017). The 'atmospheric river' phenomenon is a natural outcome of the misalignment of the deep convection region and mid-latitude stationary planetary waves on synchronously rotating terrestrial planets in such a transition circulation regime (Noda et al. 2017;Haqq-Misra et al. 2018). Having used last-saturation statistics to elucidate the behavior of atmospheric water vapor as a function of rotation on synchronously rotating planets, we now use these results to understand the multiple moist climate states presented in Section 3. While we used the aqua-planet assumption in simulations with tracers of last saturation in Section 4, the basic analysis also applies to more arid situations. In the case where only substellar water is present, the water vapor profile in the nightside troposphere can be roughly separated into two regions by the altitude of the core of the overturning cell, h cell . Above the altitude of h cell , water vapor was mainly last saturated at the substellar tropopause (Figure 4a and Figure 7a). Let the tropopause height be h trop and the saturation water vapor concentration there be q trop . Below the altitude of h cell , water vapor was last saturated in warmer regions in the atmosphere (Figure 4b and Figure 7b). Assume the water vapor concentration there is roughly q natm . On synchronously rotating planets with substellar surface water, h trop > h cell , q trop < q natm . In addition to the cold traps in the atmosphere, the nightside surface also behaves as an important cold trap where atmosphere water vapor may condense on the surface directly forming an ice layer. Let the saturation water vapor concentration at the nightside surface cold trap be q nsf c . It is the comparison of the cold-trapping strength between the substellar atmosphere and the nightside surface (measured by q trop , q natm , q nsf c ) that gives rise to the three moist climate regimes shown in Figure 1.
• Substellar oasis state. When q nsf c > q natm , water vapor never condenses on the nightside surface anymore. Precipitation forms within the upwelling branch of the overturning cell, leading to the formation of a substellar oasis surrounded by dry land (Figure 8a). This climate regime resembles Earth's tropical rain belt (Intertropical Convergence Zone, ITCZ) surrounded by the subtropical deserts and also Titan's polar methane lakes surrounded by tropical dunes (Mitchell & Lora 2016). Similar to on Earth and Titan, deep convective clouds form within the upwelling branch of the overturning cell. Some cloud particles will be carried by the upper-level divergent flow to the terminator, which may have consequences for transit observations (Komacek et al. 2020;Pidhorodetska et al. 2020;Suissa et al. 2020).
• Nightside icecap state. When q nsf c < q trop , the substellar tropopause loses its cold trap role.
Any substellar water area will slowly migrate to the nightside and eventually be cold trapped as ice caps by the freezing nightside surface. In the equilibrium state, there will be no condensation in the atmosphere and the dayside hemisphere is completely dry (Figure 8c). In Figure 1, this state shows up when the nightside surface is cold because of the lack of insolation and weak greenhouse effect.
• Transition regime. The above two climate regimes are very similar to the the bistable moist climate states discussed for close-in arid exoplanets in Leconte et al. (2013), except that Leconte et al. (2013) focused on planets receiving incoming radiation much higher than the runaway greenhouse threshold. Ding & Wordsworth (2020) discussed the above two cases for synchronously rotating planets in the habitable zone, and how they can be regulated by geological processes. However, there is in addition a third regime where q trop < q nsf c < q natm and the cold traps at both the nightside surface and substellar tropopause are effective (Figure 8b). Since q nsf c < q natm , the surface water inventory is cold trapped by the nightside surface, same as in the nightside icecap state. On the other hand, because q nsf c > q trop , tiny amount of condensation still occurs near the substellar tropopause, producing a thin layer of cirrus clouds.
Whether the tiny amount of precipitation that these clouds produce can reach the substellar surface depends on the competition between the falling speed and evaporation of raindrops, both of which are sensitive to the raindrop size (Loftus & Wordsworth 2021). This question cannot be answered in large-scale GCMs, and should be studied in numerical models that both resolves convection scale motions and employs appropriate microphysical schemes. So this transition climate regime is more complicated than the regimes with a single cold trap. It bears similarities both with present-day Mars with polar water ice caps (Montmessin et al. 2017), and a cold early Mars scenario with both polar and highland ice (Wordsworth et al. 2013).
To sum up, it is the competition of cold-trapping between the substellar upper atmosphere and nightside surface that gives rise to the emergence of the multiple moist climate states on synchronously rotating arid planets seen in Figure 1. Now we can use this concept to explain why surface water tends to be trapped on the nightside on fast rotating planets. First, q natm is slightly larger on fast rotating planets with the presence of substellar water (Figure 2), which was elucidated by our tracer of last saturation analysis in Section 4. Second, the nightside surface cold spots on fast rotating planets are much colder than the uniform nightside surface temperature on slow rotating ones (Figure 6a and b), indicating a much smaller q nsf c . Both larger q natm and smaller q nsf c on fast rotating planets make the surface water more likely to be trapped on the nightside as ice caps.
In a real climate, many factors that are not included in our idealized GCM simulations could affect the cold-trapping competition between the substellar upper atmosphere and nightside surface by influencing the tropopause temperature and nightside surface temperature. These factors include real-gas radiative transfer (Ding & Wordsworth 2020), strong shortwave absorbers in the upper atmosphere such as ozone, methane and aerosols, radiative effects of clouds, turbulent mixing in the highly stratified thermal inversion layer above the nightside surface (Joshi et al. 2020). It will be interesting to explore these effects with more complex GCM simulations in the future.

Asynchronously rotating arid planets
Arid planets with non 1:1 spin-orbit resonance will experience diurnal cycles, and will also have seasonal cycles if their obliquities are non-zero. Here we briefly discuss their likely multiple moist climate states. Although this topic is worth exploring in detail using 3D GCMs with more complex land models in the future, the basic cold-trapping principles for determining the equilibrium moist climate state described in the last section will still apply.
When the surface cold trap (i.e., the coldest surface area in the long-term averaged climate) dominates the hydrological cycle, permanent ice caps will inevitably be drawn there. On arid planets with low obliquity, this will be the polar region; on planets with obliquity higher than 45 • , it moves to the equator (Pierrehumbert 2010, Chapter 7.3). The former case resembles present-day Mars with polar water ice caps. Because of the lack of a massive moon and the proximity of Jupiter, the obliquity of Mars evolves chaotically and can drift to values as large as 47 • (Laskar et al. 2002). When Mars's obliquity is high, the surface ice migrates to the tropics -a process which has resulted in the formation of glaciers on Mars in the recent past (Head et al. 2005;Forget et al. 2006).
When the substellar tropopause cold trap instead dominates the hydrological cycle, surface water should appear under the upwelling branch of the thermally direct overturning circulation, which is usually the warmest place on the planet. On arid planets with obliquity less than 10 • , it is the tropics; on planets with obliquity higher than 45 • , the warmest and coldest place over the course of an astronomical year are the summer and winter pole, respectively (Pierrehumbert 2010, Chapter 7.3). The latter case resembles present-day Titan, with methane replacing water as the working condensible component. For instance, the Imaging Science Subsystem (ISS) and the Visual and Infrared Mapping Spectrometer (VIMS) on Cassini observed rapid migration of convective methane clouds from the high latitudes of one hemisphere to the other during the solstice season. Cassini ISS also observed a large dune field near the equator, implying dry surface condition at low latitudes (Mitchell & Lora 2016). Recently, MacKenzie et al. (2019) observed the disappearance of polar 'phantom lakes' or shallow ponds in the winter season and interpreted this as a result of methane evaporation and infiltration into a porous regolith. This climate regime on arid planets with strong seasonal cycles strongly resembles the substellar oasis state on synchronously rotating planets, with the 'substellar oasis' oscillating between the two poles in this case.
The transition regime subject to diurnal and seasonal cycles when both the atmospheric and surface cold traps are effective becomes even more complex. But it is relevant to many interesting climate evolution topics and should be investigated in more detail in future studies. These topics include the ice distribution and episodic warming on early Mars (Wordsworth et al. 2013, and the delayed onset of runaway greenhouse on terrestrial planets only covered by polar surface water (Abe et al. 2011;Kodama et al. 2018).

CONCLUSIONS
We studied cold trapping in idealized GCM simulations of synchronously rotating exoplanets using tracers of last saturation and discussed the implications for long-term surface water evolution on arid planets. Our main findings are as follows: 1. For the nightside upper troposphere above the center of the overturning cell, water vapor is mainly cold trapped by the substellar tropopause, regardless of the planetary rotation rate.
2. For the nightside lower troposphere, water vapor is mainly cold trapped by the region right above the center of the overturning cell on slow rotating planets. On fast rotating planets, the water vapor abundance in the lower atmosphere is higher than on slow rotating planets given the same condition, because the last saturation locations include the warm lower atmosphere near the terminator due to planetary wave perturbations.
3. On synchronously rotating arid planets with limited surface water inventories, there are multiple moist climate states. We explored the parameter space over a wide range of surface pressures, optical thicknesses of the non-condensible components, and planetary rotation rates. Our results show that fast rotating planets tend to trap the surface water inventory on the nightside surface as ice caps.
4. Multiple moist climate equilibrium states emerge from the competition between cold-trapping strength in the substellar upper atmosphere and on the coldest regions of the surface. This broad concept can be applied to both synchronously rotating and asynchronously rotating planets.

ACKNOWLEDGMENTS
We thank the two referees for thoughtful comments that improved the manuscript. We thank Daniel Koll for providing the Python scripts that convert GCM output from standard/Earth-like coordinates into a tidally locked coordinate system (https://github.com/ddbkoll/tidally-locked-coordinates). To validate that our analysis with tracers of last saturation on aqua-planets can be equally applied to arid planets with substellar water, we conduct another two simulations. The GCM parameters are the same as in both aqua-planet runs, except that the surface water only distributes within 40 • around the substellar point. Figure 9 shows the time and zonal mean states of water vapor distribution. Compared to the two aqua-planet runs, the lower atmosphere above the dayside dry land around the substellar water becomes much drier because of the lack of water vapor fluxes from the surface. However, the water vapor distribution elsewhere is very close to those in the two aquaplanet runs shown in Figure 2. This indicates that the key physical processes controlling the nightside water vapor buildup in the arid planet simulations with substellar water is very similar to those in the aqua-planet simulations. To better understand the physical mechanism of water vapor transport to the nightside, we continue the aqua-planet simulations for another 30 days and use the 6-hourly snapshots to create an animation of water vapor distribution at 670 hPa ( Figure 10). In the slow rotating simulation, Figure 10a shows the deep convection region is controlled by convergent flow towards the substellar point. So water vapor in the remaining subsidence region has to be advected in the upper atmosphere, which is consistent with the last saturation analysis in Figure 4b. However, in the fast rotating simulation, water vapor in the lower atmosphere is carried by the horizontal wind flowing away from the deep convection region both north and south of the deep convective region, resulting in higher water vapor concentrations on the nightside. Other than stationary flows, this animation also present interesting transient wave activities in both simulations, especially on the nightside, which may be responsible for the weak horizontal moisture gradient on the nightside (Figure 2a and b). There are also granular structures in the specific humidity convection cells that are probably associated with grid-scale condensation (our idealized GCM has no sub-grid convection parameterization). Similar granular structures were also observed in the idealized GCM simulations that explored the effects of idealized convection schemes on Earth's tropical circulations (Frierson 2007