Stable Existence of Tropical Endorheic Lakes on Titan

Cassini detected numerous hydrocarbon seas/lakes in the polar region of Titan and wide areas of sand seas in the tropical region, which initially led to the perception that Titan's tropical region may be too dry for lakes. Yet, a few tropical lakes were possibly seen on Cassini's near‐infrared images, while they were not seen by other imaging instruments. This study shows by a lake balance model in combination with a global climate model and global topography map that a few lakes can perennially exist in Titan's tropical drainage basins of Shangri‐La, Tui Regio, and Hotei Regio. This is possible because the lakes are fed by precipitation in a huge catchment area, while efficient lake evaporation occurs only in a small area inside of deep topographic depressions. However, tropical lakes may occasionally desiccate due to orbitally forced changes in tropical precipitation.


Introduction
Titan's surface is unique in that parts of the globe are covered by liquid hydrocarbon seas/lakes. At present, all seas and most lakes are found in the polar region, especially in the north (Aharonson et al., 2009;Hayes, 2016;Lorenz et al., 2014;Stofan et al., 2007). By contrast, wide areas of sand dunes were discovered in the equatorial region (Lorenz et al., 2006), which initially led to the perception that the equatorial region may be too dry for lakes ). This view was seemingly supported by climate models, which predicted drying of the equatorial surface due to an excess of evaporation over precipitation (Mitchell, 2008;Schneider et al., 2012).
Yet, possible tropical lakes were seen on near-infrared spectral surface images acquired by Cassini visual and infrared mapping spectrometer (VIMS) (Griffith et al., 2012), although this observation was not confirmed by higher-resolution synthetic aperture radar images or altimetry data sets in this area (Lopes et al., 2019). In addition, a few putative paleolakes were identified in other tropical areas (Tui and Hotei Regiones in Xanadu) (MacKenzie et al., 2014;Moore & Howard, 2010). A wet climate in the far past (Moore et al., 2014) or nitrogen rainfall in a hypothetical pure nitrogen atmosphere in the past (Charnay et al., 2014) was proposed as a possible mechanism for the formation of such tropical paleolakes.
Previously, the possible presence of extant tropical lakes was explained by a subsurface methane source given the dry tropical climate (Griffith et al., 2012). However, possible tropical lakes on Titan may not necessarily point to exotic scenarios in the past or present. In the case of Earth's arid areas, terminal lakes can exist in topographic depressions of enclosed inland basins (endorheic basins) even if the annual potential evaporation (evaporation that would occur if the surface were perennially covered by liquids) per unit area by far exceeds the annual precipitation per unit area (Löffler, 2014). Famous examples for terrestrial endorheic lakes in dry climates are the Caspian Sea, Aral Sea, Dead Sea, Chott Melrhir, and Great Salt Lake. These lakes possess elevated salt content, which decelerates the lake desiccation. In the case of Titan's lakes, dissolved nitrogen and ethane may play a similar role in reducing the methane evaporation rate. This calls for a reconsideration of the cause of possible tropical lakes on Titan. Shangri-La and Xanadu are large tropical endorheic basins whose bottom is up to ∼1 km deeper than the drainage divides surrounding the basins (Figure 1). Not only do such endorheic basins prevent liquid precipitation from flowing off towards the deep polar basins, but they also reduce the lake surface area by collecting the liquids at the bottom of the basin and thereby reduce the total evaporation from the lake, which in turn helps maintain the lake volume.
This study investigates by a lake balance model with meteorological input from a global climate model (GCM) and considering the global topography map where and under which conditions tropical lakes may exist.

Methods
A lake balance model (e.g., Mason et al., 1994) is based on the principle that a terminal lake responds to an increase/a decrease in runoff by expanding in area until the lake surface area becomes large/small enough to balance the total runoff from precipitation in the catchment area by evaporation from the lake surface. As a general rule, the lake volume increases when the precipitation integrated over the entire basin area exceeds the evaporation integrated over the entire lake surface area and vice versa. A similar model was previously applied to Titan's northern seas (Lorenz, 2014) as well as Ontario Lacus near the south pole (Dhingra et al., 2018). This approach bypasses an explicit simulation of detailed hydrological processes within the drainage basin such as surface runoff, percolation, or groundwater flow as was previously done by Horvath et al. (2016) and Faulk et al. (2019).
The stability of Titan's tropical lakes in three similar-sized tropical endorheic basins (Shangri-La, Western Xanadu and Eastern Xanadu) was investigated by a lake balance model considering the global topography map obtained by Cassini (Corlies et al., 2017) and meteorological input from author's GCM (Tokano, 2019). The total volume of liquid methane in an endorheic basin is predicted by solving the equation where V is the total volume of liquid methane, c is the runoff coefficient (a value between 0 and 1), P is the precipitation rate per unit area, E is the evaporation rate from the lake surface per unit area, A is the area over which to integrate, A B is the constant drainage basin area, and A L is the variable lake surface area where A L < A B .
The model assumes that a certain fraction of precipitation is lost to the atmosphere by evaporation from the land surface and methane table in the subsurface. The remaining fraction of precipitation, cP, eventually drains to the lakes either as surface runoff or subsurface flow after infiltration. The model does not distinguish between different forms of runoff and does not explicitly simulate the transport process itself. Infiltration of rainfall within the catchment area and subsequent subsurface methane flow to the lake does not explicitly appear in the above equation but is implicitly contained in the term cP. On the other hand, subsurface methane flow across drainage divides that surround each endorheic basin is neglected. This is because subsurface methane tables are predicted to be located at 25 m or less below the surface (Faulk et al., 2019), but the drainage divides are several hundred meters higher than the predicted lake level, so that subsurface methane flow inside the basin is everywhere directed towards the lakes. In the absence of empirical knowledge concerning c on Titan, it has to be arbitrarily prescribed. The baseline simulation assumes a runoff coefficient of 1, which is the maximum possible value and implies no evaporation of runoff before reaching the lake. Further simulations are carried out with smaller runoff coefficients.
The model performs a time-marching simulation of the seasonal cycle with a time step interval of 10 Titan days. The model does not seek for a steady state at any given time but is run until a repeatable seasonal cycle sets in, which may be regarded as a steady state. The time lag of surface and subsurface runoff behind precipitation is neglected for simplicity. Since the subsurface flow tends to be more uniform in time than the surface runoff, this simplification may cause the seasonal variation in lake volume to be larger than it should be. At each time step, the model first predicts the instantaneous lake volume using equation (1). Once the lake volume is predicted, each endorheic basin is filled hydrostatically from bottom up to the lake level below which the integrated basin volume equals the instantaneous volume of the accumulated liquid methane, analogously to the previous modeling of shrinkage of a hypothetical methane ocean in the past (Larsson & McKay, 2013). The relationship between predicted lake volume, lake surface area, and elevation, which is specific to each basin, is shown in Figure S1 of the supporting information. Small-scale basins on elevated terrains are not filled with liquids before all lower basins are filled. This is a simplification made to avoid further subdividing the catchment area based on small elevation differences. However, since the observed local basins are typically two orders of magnitude smaller than the major endorheic basin (Figure 1b), lakes in local basins would also be two orders of magnitude smaller than the terminal lakes and thus represent a tiny fraction of the total lake budget. The model topography has a horizontal grid resolution of 0.25 • as with the published data set (Corlies et al., 2017). The vertical grid resolution is 10 cm.
The model is applied to three adjacent tropical endorheic basins separately (Figure 1). The drainage divides separating the tropical endorheic basins from the polar drainage basins are found by tracking major mountain ridges on the global topography map. The drainage divides between the three tropical endorheic basins (Shangri-La, Western Xanadu, and Eastern Xanadu) are constructed such that they are not crossed by the observed drainage patterns (Burr et al., 2013;Langhans et al., 2012). The basin of Shangri-La covers an area of 5.2 × 10 6 km 2 (6.2% of the globe) roughly between 45 • S and 10 • N and 150 • W and 200 • W. The observed possible lakes in Shangri-La (Griffith et al., 2012) are located near the center of this basin and the landing site of the Huygens probe (10 • S, 192 • W) is located in the western part of this basin. Western Xanadu covers an area of 4.3 × 10 6 km 2 (5.2% of the globe) roughly between 100 • W and 150 • W and 40 • S and 5 • N. Tui Regio, where putative paleolakes were identified (Moore & Howard, 2010), is located in this basin and is rich in fluvial valleys (Burr et al., 2013;Langhans et al., 2012). Eastern Xanadu covers an area of 5.7 × 10 6 km 2 (5.7% of the globe) roughly between 60 • W and 100 • W and 45 • S and 15 • S. It contains Hotei Regio, another region where putative paleolakes were identified (Moore & Howard, 2010).
The precipitation rate as a function of season and latitude is adopted from the output of the nonuniform geography version of the Cologne (Köln) Titan GCM (Tokano, 2019). This GCM was run considering the observationally constrained global map of topography, surface albedo, surface thermal emissivity, and surface thermal inertia. A cyclic seasonal pattern is assumed; that is, the same precipitation pattern is repeated from year to year without assuming interannual variations. The latitude-season section of precipitation as well as of insolation, near-surface methane mole fraction, near-surface wind speed, and near-surface air temperature required for the calculation of the evaporation rate pattern is depicted in Figures S2, S3, and S4. Rainfall is assumed to consist of 80% CH 4 and 20% N 2 .
The lake evaporation rate per unit area is calculated as a function of lake temperature, wind speed, and air humidity . Evaporation is assumed to stop during rainfall because of the high near-surface air humidity. The lake temperature is predicted by a lake energy balance model (Tokano, 2009) considering insolation, thermal emission, evaporative cooling, and sensible heat flux. The total evaporation rate from the lake surface is given by the product of lake surface area and evaporation rate per unit area.

Equilibrium Lake Distribution
Starting from initially empty basins, precipitation increases the lake volume from year to year for a while (Figure 2a). The year-to-year increase of lake volume gradually diminishes, and eventually, a steady state is reached (by ∼3,000 years), after which a repeatable seasonal cycle is established. The total annual precipitation equals the total annual evaporation in this state. The predicted methane accumulation cannot be explained if one simply compares the annual precipitation (P) and evaporation (E) per unit area with each other. P is only ∼10 mm per Saturn year near the equator (Figure 2d), whereas E amounts to 4,200-4,400 mm per Saturn year, so that P ≪ E. Methane accumulation is only possible because liquids from huge catchment areas are collected in topographic depressions and evaporation occurs only in a small portion of the catchment area. To demonstrate that the lake accumulation is indeed supported by the basin topography, another simulation is run in which precipitation is stored on a hypothetical flat plain with the same area. Under this condition, evaporation occurs everywhere regardless of the lake volume. This also means that the effective area of lake evaporation is typically several orders of magnitude larger than in real basins. Consequently, lake evaporation can balance the runoff from the very beginning on and thereby enable seasonal but no perennial methane deposits (bottom curve in Figure 2a).
If provisionally all precipitation is assumed to drain to the lake without evaporation before reaching the lake (runoff coefficient of 1), the model predicts several large lakes in depressions of all three endorheic basins (Figure 2b). The largest lake in Shangri-La (15-19 • S, 163-174 • W) is located close to the site where possible lakes were spotted (Griffith et al., 2012). The model thus can naturally explain a stable lake in the deepest pit of this endorheic basin as a consequence of a balance between drainage of precipitation to the basin and evaporation from the lake surface, so that a local subsurface methane source is not required. However, the predicted lake in Shangri-La is one order of magnitude too large to be consistent with the observation and the model predicts several lakes in Western and Eastern Xanadu, which were not seen by the Cassini teams. Hence, the predicted lake size for c = 1 may be regarded as unrealistic.  Tables S1 and S2 of Turtle et al. (2018). + denote clouds observed by ISS, * the clouds observed by VIMS. No distinction is made between the observed cloud types (convective, discrete, streak) since it is unknown which clouds produced rainfall. For elongated clouds, the central longitude/latitude is given.
It is not likely that the excessive lake size points to an excessive precipitation from the GCM used. As shown by a GCM intercomparison study , Titan GCMs predict a large variability in the precipitation pattern/rate, but the GCM used (Tokano, 2019) tends to predict overall weaker equatorial precipitation than other GCMs. On the other hand, a runoff coefficient of 1 may well be unrealistically large since it means no evaporation of runoff en route to the terminal lake. The model predicts a monotonic shrinkage of the predicted lake size with decreasing runoff coefficient. A runoff coefficient of 0.1 is found to approximately reproduce the observed lake size in Shangri-La (Figure 2c). This means that a considerable fraction of precipitation in tropical areas should evaporate before reaching the endorheic lake.
According to the model output for c = 0.1, Shangri-La accommodates an oval lake of ∼100 km length (∼820 km 3 volume, ∼8,500 km 2 surface area, maximum depth 238 m) in a deep pit near 16 • S, 172 • W (Figure 3a), and three tiny lakes. This lake is located inside a hummocky terrain (Figure 3d). The largest lake in Shangri-La has a volume comparable to that of Lake Titicaca, the largest lake in South America. Lakes in Western Xanadu (220 km 3 ) are concentrated in Tui Regio (17-21 • S, 110-122 • W) and are distributed in eight separate lakes (Figure 3b). Lakes in Eastern Xanadu (245 km 3 ) are concentrated in Hotei Regio (22-26 • S, 84-91 • W) and, as in Tui Regio, exist as a group of lakes (Figure 3c). The predicted total lake surface area in Tui and Hotei Regiones (7,000-8,000 km 2 each) is comparable to that in Shangri-La because of the similarity in the size of the catchment area and precipitation rate. However, the lakes in Tui and Hotei are shallower (at most ∼100 m) as with the topographic depressions. All predicted lakes in Tui and Hotei are located in areas classified as plains (Figures 3e and 3f). Both Tui and Hotei are rich in evaporites (MacKenzie et al., 2014), which are indicative of the former presence of lakes, but fluvial valleys are more abundant in Tui than in Hotei (Burr et al., 2013;Langhans et al., 2012). The total volume of all tropical lakes amounts to ∼1,300 km 3 , which is ∼2% of the estimated volume (∼70,000 km 3 ) of polar seas/lakes (Hayes, 2016). The lake size varies with season by a few percent because of the seasonal cycle of precipitation and evaporation (Figures 4a and  4b). During most of the Cassini mission, this lake is predicted to have been shrinking but quickly expanding shortly before the end of the mission. The predicted shoreline recession during the Cassini mission is ∼6 km, which is too small to be observed by Cassini VIMS but possibly observable by Cassini RADAR, if the lake exists after all.

Long-Term and Short-Term Variations
The predicted tropical endorheic lakes do not completely desiccate under the present conditions unless a drought lasts for more than 1,200 years ( Figure 4c). Given the fact that equinoctial methane storms were observed within only 13 years of the Cassini mission Schaller et al., 2009;Turtle et al., 2011), such a long drought in the entire catchment area does not appear realistic in the present epoch.
However, on such long timescales, orbitally forced climate variations may come into play (Aharonson et al., 2009). While the orbitally forced insolation pattern changes at low latitudes are moderate, the annual precipitation near the equator is predicted to have varied by orders of magnitude within one Croll-Milankovitch cycle according to the GCM of Tokano (2019) (Figures S5 and S6). This directly affects the tropical lake size (Figure 4d). The present lake size turns out to be close to the long-term maximum but is still 25-40 times smaller than the inferred paleoseas (MacKenzie et al., 2014;Moore & Howard, 2010). Thus, orbital forcing alone is not likely to have caused the observed large tropical paleoseas. On the other hand, tropical lakes are predicted to have been empty 5 Kyr BP (before present) when the predicted precipitation was minimal.
Short-term variations in tropical lake size may be more limited. If one assumes a single equinoctial tropical storm, which temporarily wets a surface area of 2 × 10 5 km 2 (Turtle et al., 2011) with a cumulative precipitation per unit area of 200 mm (Hueso & Sánchez-Lavega, 2006), the lake size in Shangri-La increases by 10.1029/2019GL086166 merely 5% (400 km 2 ) ( Figure 4c). Total precipitation, be it sporadic convective rain or weak steady rain, is more relevant for the lake budget than a singular storm.

Implication for the Methane Hydrology
As already mentioned, the presence of extant lakes in Shangri-La (Griffith et al., 2012) was not confirmed by the radar (Lopes et al., 2019) and no extant lakes in Xanadu were reported. If there is indeed no extant lake in any of the three tropical endorheic basins, this may generally imply that the methane sink in the model is underestimated relative to the methane source. If the subsurface is porous, it appears possible that much of the precipitation infiltrates but is eventually lost to the atmosphere by evaporation from the methane table.
In terms of equation (1), this means that c is even smaller than assumed in this study. If, however, lakes exist only in Shangri-La but not in Tui and Hotei Regiones, this could point to a substantial difference in the ground porosity between the basins. Xanadu was characterized to contain exposed fractured bedrock (Radebaugh et al., 2011); that is, the ground porosity may be higher than in Shangri-La and thereby cause faster evaporation.
On the other hand, it may be more difficult to explain the possible difference in the lake presence between Shangri-La and Xanadu by mesoscale meteorology: If the mountain ridges constituting the drainage divides surrounding Xanadu were higher than the lifting condensation level of ∼5 km (Tokano et al., 2006), forced lifting along the mountain ridge could cause orographic precipitation on the windward side and a dry föhn on the lee side. These mountain ridges, however, have heights of merely ∼1 km at most relative to the basin interior ( Figure 1). Under such condition, orographic precipitation does not develop on the windward side but instead further from the mountain ridge on the lee side (Barth, 2010); that is, there would be no rain shadow characteristic of a föhn. Despite the inclusion of topography, the GCM used in this study indeed predicts no obvious rain shadow in tropical endorheic basins (Figure 2d). Clouds observed in this region by Cassini (Turtle et al., 2018) (Figure 2d) exhibit a slight predominance near the southern drainage divide of the three tropical endorheic basins (30-45 • S), yet clouds occurred on either side of the drainage divide; that is, the cloud distribution is not indicative of föhn. Finally, widespread morning drizzle over Xanadu was suggested based on ground-based observations (Ádámkovics et al., 2007).
As a concluding remark, the composition of minor species and deposits in possible tropical lakes may be substantially different from those of the polar seas since tropical lakes do not collect polar precipitation and they may desiccate more frequently than the northern polar seas. It remains to be seen under which past conditions tropical lakes became so large to leave geomorphologic evidence of shorelines (Moore & Howard, 2010) and vast areas of evaporites (MacKenzie et al., 2014).