# Composition, Roughness, and Topography from Radar Backscatter at Selk Crater, the Dragonfly Landing Site

Léa E. Bonnefoy, Antoine Lucas, Alexander G. Hayes, Sébastien Rodriguez, Valerio Poggiali, Daniel E. Lalich, Ralph D. Lorenz, Alice Le Gall

► To cite this version:

Léa E. Bonnefoy, Antoine Lucas, Alexander G. Hayes, Sébastien Rodriguez, Valerio Poggiali, et al.  
Composition, Roughness, and Topography from Radar Backscatter at Selk Crater, the Dragonfly  
Landing Site. The Planetary Science Journal, 2022, 3, 201 (19pp). 10.3847/psj/ac8428.  
hal-03765696

HAL Id: hal-03765696

https://hal.science/hal-03765696v1

Submitted on 31 Aug 2022

**HAL** is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
OPEN ACCESS

<figure><img src="image_2.png" /></figure>

CrossMark

# Composition, Roughness, and Topography from Radar Backscatter at Selk Crater, the Dragonfly Landing Site

Léa E. Bonnefoy<fnref n="1" />,<fnref n="2" /> ID, Antoine Lucas<fnref n="1" /> ID, Alexander G. Hayes<fnref n="2" /> ID, Sébastien Rodriguez<fnref n="1" /> ID, Valerio Poggiali<fnref n="2" /> ID, Daniel E. Lalich<fnref n="2" /> ID, Ralph D. Lorenz<fnref n="3" /> ID, and Alice Le Gall<fnref n="4" />,<fnref n="5" /> ID

Received 2022 April 21; revised 2022 July 19; accepted 2022 July 20; published 2022 August 29

## Abstract

The Selk crater region is the future landing site of NASA's Dragonfly mission to Titan. The region was imaged by the Cassini RADAR at incidence angles from 5° to 72° and at various polarization angles. Using this data set, we mapped six terrain units and assembled a backscatter curve for each, providing normalized backscatter cross section ($\sigma^0$) as a function of incidence angle. By fitting these backscatter curves with a sum of a quasi-specular and diffuse terms and evaluating three alternative formulations of the first and two for the second, we extracted the best-fit surface effective dielectric constant, rms slope, and scattering albedo. Although the parameters' absolute values are model dependent, relative values between terrains indicate real variations in surface properties. The results are consistent with the impact exposing and fracturing a low-loss tangent material such as the water-ice bedrock, which is likely also present in the hummocky terrains and to a lesser degree in the plains and interdune regions. The dunes and dark terrains are composed of smooth, uniform material with low dielectric constant (1.5–2.3 median values for all models) compatible with organic sand. A diffuse single-scattering model enabled independent derivation of the dielectric constant from high-incidence observations, leading to low values (<2) over all terrains, indicating a depolarizing (sub)surface. Finally, radarclinometry revealed lateral variations in rim height, which remains below 300 m along the SARTopo profile but reaches up to 600 m at other locations, hinting at a rim less eroded than previously thought.

*Unified Astronomy Thesaurus concepts: Radar astronomy (1329); Titan (2186); Planetary surfaces (2113)*

*Supporting material: figure set*

## 1. Introduction

The Cassini–Huygens mission, which orbited Saturn from 2004 to 2017, revolutionized our understanding of Titan, revealing a remarkably Earth-like surface featuring mountains, dunes, valley networks, craters, and a complex methane/ethane cycle involving lakes, seas, rivers, clouds, and rain (e.g., Elachi et al. 2005; Tomasko et al. 2005; Lorenz et al. 2006; Hayes et al. 2018; Turtle et al. 2018; Lopes et al. 2019). The Dragonfly mission, selected through NASA's New Frontiers program, will deepen our understanding of Titan's chemistry and geology by sending a rotorcraft to its equatorial dune fields in the mid-2030s (expected launch in 2027), approximately one Titan year after the landing of ESA's Huygens probe (Turtle et al. 2019; Lorenz et al. 2021). The landing site, in the Shangri-la dune field near Selk crater (6.5°N, 161.5°E), was chosen primarily to examine the prebiotic chemistry of the site, where interactions between organics and liquid water likely occurred within the impact melt, and for its Earth-facing equatorial position allowing direct-to-Earth (DTE) communication (Lorenz et al. 2021). Furthermore, the data coverage from Cassini over this location is substantial and suggests the

presence of smooth interdunes and/or plains that will constitute a safe landing site for Dragonfly.

The Selk crater region was first examined by Soderblom et al. (2010) using high-resolution Cassini VIMS observations and available Cassini SAR images; Malaska et al. (2016b) later mapped the same region, with a more complete data set that included additional Cassini SAR data. The crater is about 80 km wide, with a rim-to-floor height of 470 m as indicated by the SARTopo data set (Werynski et al. 2019; Hedgepeth et al. 2020). This shallow topography for such a large crater is similar to that found on other Titan craters and indicates degradation through both fluvial and aeolian erosion, as also evidenced by the presence of dunes on the crater floor (Hedgepeth et al. 2020) and possibly river channels incising some portions of the crater rim (Soderblom et al. 2010). The crater and the radar-bright terrain encircling it (likely ejecta blanket; e.g., Malaska et al. 2016b) are surrounded by extensive dune fields, including rare cross-hatched dunes located to the north and south of the crater (Malaska et al. 2016a).

The highest-resolution images on Titan (240 m to >1 km pixel⁻<fnref n="1" />) have been acquired by the Cassini RADAR, which operated at 13.78 GHz (2.16 cm wavelength) in Synthetic Aperture Radar (SAR) mode. The radar brightness of a surface, quantified by the normalized radar cross section (NRCS) $\sigma^0$, varies with incidence angle in different ways for different terrains, depending on the local topography, the composition (dielectric constant), and the structure (surface roughness, grain size, and subsurface homogeneity) of the medium. By combining observations over the same terrain at a variety of incidence angles, it is

<figure><img src="image_3.png" /></figure>

Original content from this work may be used under the terms of the [Creative Commons Attribution 4.0 licence](https://creativecommons.org/licenses/by/4.0/). Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.

<footnote n="1">¹ Institut de physique du globe de Paris (IPGP), CNRS, F-75005, Université Paris Cité, Paris, France</footnote>

<footnote n="2">² Cornell Center for Astrophysics and Planetary Science (CCAPS), Cornell University, Ithaca, NY 14850, USA</footnote>

<footnote n="3">³ Applied Physics Laboratory, Johns Hopkins University, Space Exploration Sector, Laurel, MD 20723, USA</footnote>

<footnote n="4">⁴ Laboratoire Atmosphères, Milieux, Observations Spatiales (LATMOS), UVSQ/CNRS/Paris VI, UMR 8190, F-78280 Guyancourt, France</footnote>

<footnote n="5">⁵ Institut Universitaire de France (IUF), Paris, 75005, France</footnote>
possible to build a backscatter curve, which can then be compared with scattering models to derive surface and volume properties. This approach has been applied using different models on Titan's dune fields (Le Gall et al. 2011; Poggiali et al. 2012; Paillou et al. 2014; Lucas et al. 2019), its empty lakes and transient liquids (Hayes et al. 2011; Hofgartner et al. 2014; Michaelides et al. 2016), and globally from the low-resolution scatterometry data set (Elachi et al. 2006; Sultan-Salem & Tyler 2007; Wye et al. 2007). Herein we apply similar methods of analysis to the Dragonfly landing site, which has been observed at incidence angles from 5° to 72°, yielding new constraints on the surface properties of the different terrain units in the region.

In this study, we thus focused on SAR data over the Selk crater region. In Section 2, we describe how we mapped six distinct terrain units from Cassini SAR observations: the crater rim, its ejecta, hummocky terrains, plains, dune fields, and dark terrains, which are mostly featureless radar-dark regions embedded within the dune fields. From the SAR data within each unit, we assembled the backscatter curve of each terrain. Section 3 describes the quasi-specular and diffuse models used to simulate the observations, the fitting process to the data, and the model sensitivity to each parameter. The resulting parameter constraints are presented in Section 4, and their implications in terms of composition and texture for the Selk crater region are presented in Section 5. The relationship between the radar backscatter and topography is examined in Section 6 for the crater rim.

## 2. The Dragonfly Landing Site Seen by the Cassini RADAR

In order to characterize the surface properties of the Dragonfly landing site, we examine the Cassini SAR data in the region of interest (ROI) defined by Lorenz et al. (2021). This region is centered on Selk crater and extends from 156° E to 166° E and from 1° N to 13° N (Figures 1 and 2). It includes Selk crater and its ejecta blanket, surrounding plains, dune fields, and hummocky terrains.

### 2.1. Radar Data over the ROI

The Cassini RADAR measured same-sense linearly polarized radiation using five distinct antenna feeds. During a typical Titan flyby dedicated to RADAR, the observation mode depended on the distance to Titan according to an observation sequence described in Elachi et al. (2004) and Janssen et al. (2009). The radar transmitter was turned ON at about 40,000 km from Titan, and scatterometry scans were performed over part of Titan; these data were processed as real aperture radar (RAR), yielding a single backscatter value per beam footprint and resulting in resolutions of tens to hundreds of kilometers per footprint. Most terrains in the ROI are not resolved in the RAR data, which was thus not considered in this study. Around 4000–12,000 km altitudes, RADAR focused on altimetry for observations at nadir (0° incidence); however, no altimetry data have been acquired within the ROI. At closest approach, all five beams were used to look at oblique incidence angle in SAR mode. During a typical acquisition of SAR data, the resolution of resultant maps varied from as low as ≈240 m near closest approach to ≈1.5 km near the end of the main SAR segments. The reverse observation sequence then occurred as Cassini receded from Titan. SAR data generally include a main swath and several high-altitude (HiSAR) segments; they are denominated by Titan flyby number, followed by the segment number.

Selk crater has been imaged by the Cassini radar in SAR mode on five occasions: T36, T95, T98, T120, and T121, with T121S06 being a very high altitude (>27,000 km) scan. The ejecta was further imaged at relatively poor resolutions during T39 and T83. Luckily, these observations over the Selk crater region include incidence angles varying from 5° to 72°, making it one of the best regions on Titan to examine backscatter variations with incidence angles. There are two additional SAR swaths within the ROI that do not directly cover Selk crater or its ejecta: T41 and T61. For the units observed within these swaths, namely, dune fields and hummocky terrains, these two extra swaths add crucial low incidence angle information. All nine swaths within the ROI are shown in Figure 1, and their resolutions and observation angles are listed in Table 1.

The observation angles include the incidence angle $\theta$ (angle between the antenna look vector and surface normal), the look direction (the projection of the antenna look vector in the plane tangent to the surface at the observation point, counterclockwise from east), and the polarization angle $p$. The polarization angle is the angle between the look direction and the projection of the electric field direction (x-axis of the Cassini spacecraft) onto the local tangent plane (see Table 1). The Cassini radar measured the same-sense linearly polarized signal; this can correspond to HH (perpendicular, 90°/270°), VV (parallel, 0°/180°/360°), or any value in between. For a typical SAR main swath observation (e.g., T61 and T95) using all five beams, the (side-looking) radar aligns the look direction with the y-axis, leading to a perpendicular (HH) polarization. However, this is not true at the beginning and end of SAR main swaths (e.g., T41 and T120 here) and during HiSAR observations, which use only one beam.

### 2.2. Mapping the ROI

Terrain units were mapped with the goal of extracting their backscatter curve. Consequently, the criteria used are based primarily on the terrain's SAR brightness ($\sigma^0$), its texture (uniform, patterned by dunes or likely hills), and the geomorphological context. This classification remains similar to that used by previous geomorphological maps of Titan's equatorial regions (e.g., Williams et al. 2011; Lopes et al. 2019; Schoenfeld et al. 2021), particularly Malaska et al. (2016b), whose map of the Afekan crater region includes the Selk crater area, and Lorenz et al. (2021), who mapped the Dragonfly landing site based on SAR, ISS, and VIMS data separately. Herein, we used only six units, defined hereafter and mapped in Figure 2. Areas with poor data resolution and high noise were mapped based primarily on radar brightness (hummocky for bright, plains for medium brightness, and dune fields for dark areas), but are identified as uncertain by cross-hatching in the map; these areas were not used for backscatter analysis.

To create the map, we used SAR data corrected for incidence angle, downloaded from the Planetary Data System (PDS Cartography and Imaging Sciences Node 2021), at the highest resolution available for each swath and segment. These data were projected to an equiangular Titan projection using ISIS3 software. All SAR swaths were co-registered to T95 (the highest-resolution SAR data directly over Selk crater) using the software ArcGIS Pro. Within the ROI, a simple translation of the data was sufficient to ensure the best possible correspondence of features in overlapping SAR data. Mapping was then performed at a viewing resolution of 1:500,000 on ArcGIS Pro.

The six mapped units are defined in order of decreasing radar brightness as follows:
<figure><img src="image_4.png" /><figcaption>**Figure 1.** The nine SAR swaths within the ROI. The incidence angle $\theta$ is given, and the look direction is shown with a magenta arrow. The NRCS $\sigma^0$ is linearly scaled from 0 to 1.5 for all nine swaths. Note that the aspect of the Selk crater ejecta changes with incidence angle, appearing brighter and more widespread in T95 than in T98, for example. The position of the ROI on Titan is shown in the context map at the top, over the VIMS/ISS background from Seignovert et al. (2019), overlaid with a global radar mosaic.</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>**Figure 2.** Left: mosaic of the incidence-angle-corrected SAR swaths within the ROI. Right: geomorphological map of the ROI. The criteria to differentiate between the six units are detailed in Section 2.2. For easy comparison with previous work, we use a color scheme similar to that of Malaska et al. (2016b).</figcaption></figure>

<table><caption>**Table 1**

Properties of the Cassini SAR Data in the Region of Interest over Selk Crater</caption><thead><tr><th>Date</th><th>Swath ID</th><th>Incidence Angle θ (deg)</th><th>Look Direction (deg CCW from East)</th><th>Polarization Angle p (deg CCW from Look: 90° is ⊥)</th><th>Azimuth Resolution (km)</th><th>Range Resolution (km)</th></tr></thead><tbody><tr><td>2007 Oct. 2</td><td>T36S4*</td><td>35–50</td><td>109–113</td><td>192–194 (≈||)</td><td>2.6–3.5</td><td>2.1–2.7</td></tr><tr><td>2007 Dec. 20</td><td>T39S3*</td><td>16–27</td><td>40–52</td><td>164–178 (||)</td><td>1.4–1.7</td><td>0.5–0.7</td></tr><tr><td>2008 Feb. 22</td><td>T41S2</td><td>5–10</td><td>160–205</td><td>115–155</td><td>0.7–1.4</td><td>1.3–2.5</td></tr><tr><td>2009 Aug. 25</td><td>T61S1</td><td>18–26</td><td>111–118</td><td>268–277 (⊥)</td><td>0.3–0.4</td><td>0.4–0.5</td></tr><tr><td>2012 May 21</td><td>T83S5*</td><td>6–12</td><td>41–42</td><td>137–138</td><td>2.8–2.9</td><td>7.0–7.7</td></tr><tr><td>2013 Oct. 13</td><td>T95S1</td><td>21–29</td><td>332–2</td><td>72–104 (⊥)</td><td>0.2–0.3</td><td>0.3–0.4</td></tr><tr><td>2014 Feb. 2</td><td>T98S2*</td><td>60–72</td><td>231–242</td><td>123–155</td><td>1.2–3.6</td><td>1.3–1.5</td></tr><tr><td>2016 Jun. 7</td><td>T120S1</td><td>62–72</td><td>150–155</td><td>120–137</td><td>0.5–0.8</td><td>0.2–0.3</td></tr><tr><td>2016 Jul. 25</td><td>T121S6*</td><td>31–48</td><td>277–295</td><td>320–350 (≈||)</td><td>3.6–6.7</td><td>1.7–2.4</td></tr></tbody></table>

**Note.** Swath IDs include the Titan flyby number and the segment number; HiSAR observations are indicated by an asterisk. The look direction is the direction of the projection of the antenna look vector in the plane tangent to the surface at the observation point, counterclockwise (CCW) from east. The polarization angle is the angle between the look direction and the projection of the electric field direction onto the same local tangent plane. Note that parallel polarizations correspond to VV and perpendicular to HH. Observations angles and resolutions are as given in the SBDR files on the Planetary Data System (PDS Cartography and Imaging Sciences Node 2021), for the burst numbers in the ROI.

1. The crater rim, which delineates the circular shape of the Selk impact crater, is the brightest terrain in the ROI, at all incidence angles. It is textured with important dark/bright contrasts, caused by the presence of topography (degraded crater rim) and/or by local variations in surface compositional and structural properties (e.g., organic material within valleys eroded through the rim). It is much brighter than the surrounding crater ejecta at high incidence angles (e.g., in T98), whereas at lower incidence angles (e.g., T95) their brightness is similar and the two units are distinguished mainly by the rim's dissected aspect. The crater rim corresponds to the topographically high region identified in

the SARTopo data set (Stiles et al. 2009; Hedgepeth et al. 2020). We note that, since this is a large and degraded impact crater, this extended, rugged, topographically high region may be proximal ejecta; nonetheless, we call it the “crater rim,” both for consistency with previous work (e.g., Williams et al. 2011; Malaska et al. 2016b) and to avoid confusion with the distal ejecta region.

2. The crater ejecta, which surrounds the crater rim, is bright and uniform, similar to hummocky terrains. It likely consists partly of distal ejecta from the Selk impact, partly of preexisting terrain, and partly of organic material brought through aeolian, fluvial, and pluvial processes.
3. Hummocky terrains are radar-bright and uniform in texture. They are often surrounded by dunes that seem to wrap around them, consistent with topographic highs, as indicated in some cases by SARTopo (e.g., Malaska et al. 2016b; Schoenfeld et al. 2021) and altimetry (Poggiali et al. 2019; Lalich et al. 2022).

4. Plains in this area are of intermediate brightness and mostly featureless. They are often located downwind of likely topographic obstacles such as Selk crater and hummocky terrains (assuming that dune direction indicates dominant sediment flux direction; Lucas et al. 2014b; Malaska et al. 2016a). Plains vary somewhat in brightness over the region, being generally brighter near radar-bright units such as hummocky terrain and darker near dune fields (Lopes et al. 2016, 2019). Part of the crater floor shows a similar backscatter response and was mapped as plains; this is further discussed in Section 5.1.1.

5. Dune fields feature linear radar-dark dunes separated by interdune regions of varying brightness. The dunes are oriented mostly west-east (longitudinally) to the south of Selk and mostly south-north (latitudinally) to the north of Selk (Malaska et al. 2016a). In some areas close to the crater, the dunes form a cross-hatched pattern, with both orientations superimposed.

6. Dark terrains are the most radar-dark terrains within the ROI and feature little to no clear dune patterns. They are likely dune fields or sand sheets: either the sand-free interdune regions are unresolved, or the sand cover is thick enough that the radar cannot detect the underlying bedrock. This unit, which is widespread on Titan and especially in the Belet sand sea, is usually included within the linear dunes unit (Malaska et al. 2016b; Schoenfeld et al. 2021). The fact that there are no apparent radar-bright interdune regions in this unit leads to very low backscatter values, often below the noise floor. As a result, the backscatter curve is distinct from that of the linear dunes unit and had to be mapped separately. This unit was mapped as "probable dunes" by Lorenz et al. (2021); however, they also included poorly resolved dune fields (e.g., in T41), which we map as dune fields with uncertain mapping.

The dunes and interdune regions are included in the “dune fields” unit but were also mapped separately by hand in the grid cells where they were well resolved and clearly separated, similar to Bonnefoy et al. (2016). Examples of areas where dunes were not mapped include dunes near obstacles, cross-hatched dune fields where the interdunes are difficult to discern, and very wide dark streaks likely more similar to sand sheets than to dunes. The dunes and interdune regions could only be mapped in data with high enough resolution: T95, T61, and small portions of T120 and T41. They are hereafter treated as two additional terrains, to be compared with the other six terrains. Note that the “dune fields” unit, which mixes the signal from dunes and interdunes together, can include poorly resolved data (T36, T39, T98, T121) and therefore presents a more complete backscatter curve in terms of incidence angle coverage.

## 2.3. Building Backscatter Curves

Even within a given terrain unit, roughness, composition, and topography are not expected to be completely uniform. This is especially obvious for the dune fields, which include

smooth, organic, sandy, sloped dunes separated by flatter, larger-grained interdune regions with potentially exposed bedrock (e.g., Barnes et al. 2008; Rodriguez et al. 2014; Bonnefoy et al. 2016; Lucas et al. 2019). The crater rim and ejecta blanket are also nonuniform, with possibly organic-rich valleys cutting through the hilly and blocky icy ejecta (Neish et al. 2018; Lorenz et al. 2021). When putting together the backscatter curve of the terrain units described above, it is therefore essential to ensure that the terrain is as homogeneous as possible in terms of backscatter and to remain aware of the unresolved diversity of the terrains during interpretation.

### 2.3.1. Separating Each Terrain Unit

We divided the ROI into a 0.25° by 0.25° grid, as in Lorenz et al. (2021). The NRCS $σ^0$ (uncorrected for incidence angle) and incidence angle $θ$ values were averaged for each unit containing >50 pixels within each grid cell. For each grid cell individually, the geomorphological context, quality of the data, and consistency with other data points were examined in order to keep only points within regions similar in terms of backscattering properties. The resulting curves, shown in Figure 3, show distinct backscattering behaviors for each terrain except the ejecta and hummocky regions, which appear similar. The uncertainties for each grid cell are the standard deviation of the data within, divided by $\sqrt{N}$, where N is the number of independent resolution cells within the terrain unit within the grid cell. The number of resolution cells is calculated not from the images, which have been resampled to provide uniform pixel scale, but rather from the reported SAR range and azimuth resolutions of the radar footprint closest to each grid cell. The SAR range and azimuth resolutions (Table 1) are provided in the SBDR files, available on PDS (PDS Cartography and Imaging Sciences Node 2021). For the lowest-resolution data (T83_S05 and T121_S06), there can be as few as four resolution cells within the terrain unit and the grid cell; this is reflected in the very large uncertainties of these data points.

We can see in Figure 3 that T121S6 exhibits higher $σ^0$ values than T36S4 at similar incidence angles and over the same regions. The T121S6 segment was acquired in scatterometry mode, while Cassini was very far from Titan (>27,000 km), yielding poor SAR resolutions. The radar scanned the ROI quickly during this segment, creating few overlapping looks and leading to high speckle noise. The high noise and low resolution explain the larger scatter seen in the T121S6 data; however, they should not cause a bias as observed. Instead, the most likely explanation is a calibration issue, which appears clearly at poor resolutions but may also be present in other data. The problem may, for example, be due to relative calibration between different beams (T121S6 and T98S2 are acquired with beam #3, whereas T36S4 and T39S3 are acquired with beam #4), to the contribution of sidelobes, or to application of the SAR processing pipeline to distant scatterometry data. Correcting the calibration issue of T121S6 would involve an in-depth revision of the calibration process and is beyond the scope of the present study; consequently, we do not include the T121S6 swath segment in the backscatter analysis after this point. The same phenomenon also appears to be present for T98S2, another distant scatterometry data set processed to obtain SAR images, although the low radar brightness at these incidence angles makes the difference between T120S1 and T98S2 not obvious. Because T98S2 is generally consistent with T120S1, we chose to keep these data; all fits were also applied without these data, resulting in the same parameter constraints within $1σ$ uncertainties.
<figure><img src="image_6.png" /><figcaption>**Figure 3.** Variations of the NRCS $σ^0$ with incidence angle $θ$ within each terrain unit mapped within the ROI. As described in Section 2.3.1, the $σ^0$ is averaged within each grid cell, and the error bars are the standard deviation divided by $\sqrt{N}$, where N is the number of independent resolution cells within the terrain unit within the grid cell. Each color corresponds to a different SAR swath, as described in the legends. Note that the crater rim backscatter curve presented in this figure is not corrected for local slopes, although it is in Figure 5. The dune and interdune region curves correspond to part of the area of the dune fields unit, with dunes and interdune regions separated only in the swaths in which they are resolved.</figcaption></figure>

In the backscatter curve of both the dunes and interdune regions, there is an anomaly near 20° incidence. Although T95 and T61 observed the same region at the same incidence angles, the resulting backscatter is different, with T61 being generally brighter. This is also true, though less visible owing to the presence of other data, for the dune fields unit (which includes dunes and interdunes mixed together). This anomaly is most likely caused by differences in backscatter with look direction over dunes, which are very geometrical regular features. The look direction is roughly perpendicular to dunes in T61, which is therefore affected by local slopes, leading to apparently thinner dunes in this swath. Due to this uncertainty in the local incidence angle, the T61 data are hereafter excluded from the dune and interdune region backscatter curves (but kept in all other units). A global analysis of overlapping SAR swaths is ongoing to confirm the azimuthal dependence of backscatter from dunes and extract their morphology.

## 2.4. Local Slopes over the Crater Rim

The SARTopo data set (Stiles et al. 2009) indicates a crater rim height of about 280 m relative to the surrounding terrains and 470 m relative to the crater floor, leading to average slopes of a few degrees (Werynski et al. 2019; Hedgepeth et al. 2020).

These slopes must be taken into account when analyzing radar observations, as they modify the local incidence angle. We built a synthetic topography model of the crater rim, which is inspired by both the SARTopo data and the radarclinometry presented in Section 6. The radarclinometry results were not directly used owing to their partial coverage, noise, and high uncertainties when the rim is seen parallel to the line of sight. Instead, we model the rim as a circularly symmetric 500 m tall smooth hill of Gaussian shape. The local slope along the radar look direction of each observation is then added to the incidence angle, and the data over the rim are binned by local incidence angle. This simple model only adjusts the local incidence angle, without accounting for differences in the pixel area caused by surface slopes, but is a good first-order correction. At incidence angles higher than 30° (that is, in swaths T36, T98, T120, and T121), we find that the backscatter values corresponding to slopes toward the radar are higher than at slopes away from the radar, confirming that we are correcting for real surface slope effects. Around 20° incidence (swath T95), we see no clear backscatter variation with slope: such constant backscatter at low incidence is consistent with a purely diffuse rather than specular surface. The backscatter curve thus derived for the crater rim is used hereafter for all fits and is shown in Figure 5. We also attempted
fitting the original data versus the nominal incidence angle, resulting in much higher uncertainties on all parameters; though poorly constrained, these fits remain consistent with the ones presented herein.

### 3. Methods: Modeling and Fitting Backscatter Curves

The dominant mechanism contributing to radar backscatter varies with incidence angle. Within the Fresnel zone, at very low incidence angles, specular reflection can lead to strong radar returns on smooth surfaces (e.g., Wye et al. 2009; Mastrogiuseppe et al. 2014; Zebker et al. 2014; Poggiali et al. 2020). Away from the Fresnel zone, quasi-specular scattering occurs on smooth surface facets oriented toward the observer. This effect is weaker at larger incidence angles, as there are fewer facets orthogonal to the incoming radiation; it becomes negligible above about 30°. At large incidence angles, diffuse scattering dominates and is controlled by surface roughness and/or subsurface heterogeneity. Subsurface, or volume, scattering due to multiple reflections on buried voids or structures is especially significant within icy and/or porous surfaces such as icy satellites, which are very transparent to radar wavelengths (Ostro et al. 2006, 2010; Le Gall et al. 2019). As a general rule, smooth and homogeneous surfaces should appear very bright near nadir and dark at high incidence angles, whereas rough, heterogeneous material should be of roughly uniform brightness for all incidence angles below about 60° (e.g., Wye 2011).

To simulate the different mechanisms responsible for backscatter at various incidence angles, most models combine a quasi-specular component $\sigma_{qs}^0$ and a diffuse component $\sigma_V^0$. When very low incidence angle data are included, a second quasi-specular term or a coherent term can be added (Rodriguez 2003; Sultan-Salem & Tyler 2007; Le Gall et al. 2011; Wye 2011); however, because the lowest incidence angles in the ROI are of 5°, a coherent model is not necessary in our case. To make sure that our results and interpretation are not model dependent, we tested several different models for both the quasi-specular and diffuse components, as described below.

#### 3.1. Quasi-specular Component

Similar to many previous Titan studies, we assumed a Kirchhoff surface scattering model, which simulates the surface as a series of facets varying vertically by an rms height $\xi$ and horizontally over a correlation length $\zeta$ (e.g., Wye et al. 2007; Zebker et al. 2008; Hayes et al. 2011; Le Gall et al. 2011; Hofgartner et al. 2014; Michaelides et al. 2016). Both the height probability function and autocorrelation function are commonly described by either Gaussian or exponential functions, resulting in the three following quasi-specular scattering laws. The Hagfors law (Hagfors 1964) assumes a Gaussian height distribution and an exponential autocorrelation function, whereas the exponential law assumes an exponential form for both and the Gaussian law assumes Gaussian forms for both:

$$ \sigma_H^0 = \frac{C\Gamma(\varepsilon_{qs}')}{2}(\cos^4 \theta + C \sin^2 \theta)^{-3/2} \quad (1) $$

$$ \sigma_E^0 = \frac{3C\Gamma(\varepsilon_{qs}')}{\cos^4 \theta} e^{-\sqrt{6C} \tan \theta} \quad (2) $$

$$ \sigma_G^0 = \frac{C\Gamma(\varepsilon_{qs}')}{\cos^4 \theta} e^{-C \tan^2 \theta}, \quad (3) $$

where $\theta$ is the incidence angle. These three quasi-specular scattering models have the advantage of being simple (allowing for fast fitting procedures) and of depending on only two parameters: the parameter $C$ and the Fresnel reflection coefficient $\Gamma$ at normal incidence, both of which can be related to physical properties of the surface: the rms slope $s$ and, for a low-loss medium, the real part of the effective relative permittivity (or dielectric constant) $\varepsilon_{qs}'$ as follows:

$$ C = \frac{1}{s^2} = \frac{1}{\tan^2 \phi_{rms}} \quad (4) $$

$$ \Gamma(\varepsilon_{qs}') = \left| \frac{1 - \sqrt{\varepsilon_{qs}'}}{1 + \sqrt{\varepsilon_{qs}'}} \right|^2 \quad (5) $$

The rms slope $s$ expresses the variation of the rms height $\zeta$ over the correlation length $\xi$. It can be expressed as the ratio of these two variables, $s = \zeta/\xi$ (Lucas et al. 2019), although other formulations exist, which vary with the autocorrelation function used (e.g., Hagfors 1966; Shepard & Campbell 1999; Shepard et al. 2001). Since we do not have sufficient data to determine separately the rms height and correlation length, we only consider the parameter $s$, which is equal to the tangent of the rms slope of the surface facets $\phi_{rms}$. For easier interpretation of our results, we hereafter discuss the values of the rms slope angle $\phi_{rms}$ in degrees.

For a nonmagnetic medium (such as, most likely, Titan's surface), the real part $\varepsilon'$ of the dielectric constant dictates the speed of light through the material ($v = c/\sqrt{\varepsilon}$), whereas the imaginary part is associated with the amount of radiation absorbed during transmission through the material. The relative dielectric constant is primarily a function of the surface composition: at 13.78 GHz, it is equal to $\varepsilon = 3.13 + i 1.3 \times 10^{-3}$ for water ice and $\varepsilon = 1.7 + i 0.8 \times 10^{-3}$ to $2.4 + i 22 \times 10^{-3}$ for various hydrocarbons and tholins, depending primarily on the degree of compaction (Rodriguez 2003; Paillou et al. 2008). Higher values (on the order of 3–4) may be found for the real part of the relative dielectric constant of water–ammonia ice mixtures. Silicate and carbonate materials have even higher dielectric constants but are not expected in large quantities on Titan's surface. In general, the materials expected on Titan (hydrocarbons, water ice, etc.) are significantly more transparent (i.e., have smaller imaginary components of their dielectric constants) than terrestrial materials (silicates, liquid water, etc.). Hereafter, we assume a low-loss medium and call “dielectric constant” the real part of the effective (i.e., including porosity) relative dielectric constant, whose minimal physical value is 1 (vacuum). The dielectric constant derived using the quasi-specular component is subscripted qs.

The Hagfors, exponential, and Gaussian scattering laws have been applied to several planetary bodies, including the Moon (e.g., Evans & Pettengill 1963; Hagfors 1964), Venus (e.g., Tyler et al. 1992), Mars (e.g., Simpson et al. 1978; Harmon et al. 1982), and Titan (e.g., Sultan-Salem & Tyler 2007; Wye et al. 2007; Hayes et al. 2011; Le Gall et al. 2011; Michaelides et al. 2016). We note that, with $C = 1/2 m^2$, the Gaussian model is the same as the geometric optics model (Hagfors 1966; Fung et al. 1992), which has been applied to Titan by Rodriguez (2003), Paillou et al. (2006), Paillou et al. (2014), and Lucas et al. (2019). These authors also examined two other, more complex and realistic models, using the physical optics method (POM)
and the integral equation method (IEM); however, they found that all Titan terrains, including the dunes, have a high enough roughness relative to the radar wavelength that they fall into the domain of validity of the GOM for the incidence angles in question. Furthermore, both the POM and the IEM require fitting separately the rms height and the correlation length, thus adding one more parameter to an already poorly constrained problem. In consequence, the POM and IEM were not applied to the Selk crater region in this work. Fractal models (e.g., Shepard & Campbell 1999; Sultan-Salem & Tyler 2006) and two-layer models (Paillou et al. 2006; Grings et al. 2021) have also been devised to describe more realistically the surface roughness, but they also add one or more parameters that cannot be independently constrained using the data available at low incidence angles.

## 3.2. Diffuse Component

For large incidence angles, it is necessary to account for diffuse scattering, which is due both to surface roughness and, in the case of low-loss materials, including water ice and hydrocarbons, to subsurface scattering from heterogeneities including voids, cracks, and embedded material (e.g., buried rocks).

The diffuse scattering model most commonly used on Titan and Saturn's other icy satellites is a cosine power law, which takes the following form (Ostro et al. 2006; Wye et al. 2007; Zebker et al. 2008; Hayes et al. 2011; Janssen et al. 2011; Le Gall et al. 2011; Michaelides et al. 2016; Le Gall et al. 2019):

$$ \sigma_D^0 = A \cos^n \theta, \qquad (6) $$

where A is the diffuse scattering albedo and n is the cosine exponent. This model is semiempirical, and its parameters are not directly related to physical properties of the surface and near subsurface. As expressed by Ostro et al. (2006), n is expected to lie between 1 (a uniformly bright scatterer) and 2 (a Lambertian surface). Values higher than 2, as found globally on Titan by Le Gall et al. (2019), would indicate an inappropriateness of the model to describe the observations. We consequently fix the upper bound of n values at 2.5, an already high value that allows for large uncertainties and, if found, points to the need for a different diffuse scattering model. This model, which has the advantages of being widely used and fitting the data well (due to the flexibility allowed by the exponent), is not physical and therefore cannot be directly linked to surface properties.

We therefore also use a simple yet physical model of the diffuse scattering component, proposed by Swift et al. (1985) and applied to sea ice (Swift 1999; Remund & Long 2003), the Greenland ice sheet (Swift et al. 1985; Ashcraft & Long 2006), and Saharan ergs (Stephen & Long 2005). For single scattering from discrete scatterers, we can thus model the diffuse component as (Swift et al. 1985; Swift 1999)

$$ \sigma_D^0 = \eta T^2 (\varepsilon_d', \theta, p) \cos \theta, \qquad (7) $$

where $T(\varepsilon_d', p)$ is the transmissivity ($T=1-\Gamma$, with $\Gamma$ being the reflectivity), which depends on the dielectric constant $\varepsilon_d'$, the incidence angle $\theta$, and the polarization angle $p$ following Fresnel's equations. Fresnel's equations provide the reflectivity for parallel and perpendicular polarizations separately; these are combined for intermediate polarization angles as follows

(Heiles & Drake 1963):

$$ \begin{aligned} T(\varepsilon_d', \theta, p) &= 1 - \Gamma(\varepsilon_d', \theta, p) \\ &= 1 - (\Gamma_{\parallel}(\varepsilon_d', \theta)\cos^2 p + \Gamma_{\perp}(\varepsilon_d', \theta)\sin^2 p). \end{aligned} \qquad (8) $$

Note that observations in parallel polarization are generally expected to be brighter, especially at high incidence angles, due to the Brewster angle (Ulaby & Long 2015). The dielectric constant derived using the diffuse component is subscripted d in order to differentiate it from the one derived from the quasi-specular component; physically, these are the same variable.

Similar to A in the $A\cos^n$ model, $\eta$ is the scattering albedo, as defined by Swift (1999):

$$ \eta = \frac{n\sigma_b}{2\alpha}, \qquad (9) $$

where n is the number density of subsurface scatterers, $\sigma_b$ is the backscatter cross section per particle, and $\alpha$ is the bulk volume attenuation coefficient. These parameters cannot be separated from the available data using this model and are all grouped together in the parameter $\eta$. The scattering albedo $\eta$ derived from the Swift model must be higher than the value of A derived from the $A\cos^n$ model in order to arrive at the same $\sigma^0$ values in spite of the transmissivity (which is <1). The two main limitations of the Swift model are the difficulty in interpreting $\eta$ in terms of physical parameters and the simplifications intrinsic to the model (e.g., no multiple scattering, no preferential orientation of scatterers).

## 3.3. Fitting Method

The uncertainties in the $\sigma^0$ values derived for each grid cell in Section 2 and shown in Figure 3 are caused by speckle noise and unresolved terrain variability. Because speckle noise has a roughly exponential distribution, the distribution of the values of $\sigma^0$ within each grid cell is rarely Gaussian (only for high numbers of looks; Lucas et al. 2014a). Therefore, the standard deviations do not accurately represent the uncertainty in the data, and a weighted fit would not correctly account for the real variability in the observed data.

Instead, we chose to apply a Monte Carlo randomized sampling of the data. For 500 iterations (near the maximum number of resolution cells within each grid cell, for all swaths), we randomly sample a single data point from the data within each grid cell. The resulting (noisy) scattering curve is fit using Matlab's fit function, which uses the Trust Region algorithm to perform a nonlinear least-squares fit. We assume uniform a priori distributions for each parameter, using physically reasonable bounds as follows: $1 < \varepsilon' < 5$, $1 < \phi_{rms} < 50^{\circ}$, $0 < A < 5$, $0 < \eta < 5$, and $0 < n < 2.5$. The parameters and best-fitting curve are recorded for each iteration, resulting in a distribution of best-fitting values for each parameter, from which we derive the median and $1\sigma$ confidence intervals.

The quasi-specular and diffuse components are fit separately for large and small incidence angles. Although each component is dominant at different incidence angles, scattering at all incidence angles is a combination of both: the fit is therefore applied in an iterative manner. We first fit the diffuse model $\sigma_d^0$ to high-incidence (>35°) data. We then fit the quasi-specular model $\sigma_{qs}^0$ to the low-incidence (<50°) data minus the newly modeled diffuse component ($\sigma_{observed}^0 - \sigma_d^0$). The diffuse fit is then repeated on high-incidence (>35°) data minus the quasi-specular component that was just estimated ($\sigma_{observed}^0 - \sigma_{qs}^0$).
<figure><img src="image_7.png" /><figcaption>**Figure 4.** Sensitivity of the fit to each of four parameters, for each combination of quasi-specular + diffuse model. For this figure, the Swift diffuse scattering model assumes perpendicular polarization angle. For parallel polarization the contribution of $\epsilon'_d$ is negligible at all incidence angles (because the transmissivity near the Brewster angles is almost the same for all dielectric constants considered).</figcaption></figure>

This process is repeated until convergence is reached (i.e., the parameters and quality of fit vary little), which we find occurs when the change in best-fitting dielectric constant is $\Delta\epsilon'_{qs} < 0.05$. This iterative process is repeated for all 500 simulated data sets. The incidence angle cutoff values ($>35°$ for high incidence; $<50°$ for low incidence) were chosen such that extreme data are excluded (very low and very high incidence angles), while still keeping as much data as possible for a reliable fit. Note that $35°$ is the cutoff value used by Wye et al. (2007), who used a similar method with only one iteration. For the rim, ejecta, and hummocky terrains, fitting both the quasi-specular and diffuse components at the same time over all the data provided a better fit, most likely because diffuse scattering is significant even at low incidence angles in these regions.

When the Swift model is used for the diffuse component, there are two independent ways to constrain the dielectric constant, as it intervenes in both the quasi-specular and diffuse components. In both cases, the dielectric constant affects the amplitude of the reflectivity through the Fresnel equations. Since the same terrain is observed, these two values of the dielectric constant should be the same, but since they are derived from different models and from different data (at low incidence angles for the quasi-specular model and high incidence angles for the diffuse model), we chose to keep them as two fully independent parameters, named $\epsilon'_{qs}$ and $\epsilon'_{d}$. The dielectric constant derived from the quasi-specular component is linked primarily to the absorbing or reflecting character of the medium (high $\epsilon'_{qs}$ for a highly reflective surface).

Meanwhile, because we have some variety of polarization angles, the value derived from the diffuse model depends primarily on the degree of polarization of the medium (low $\epsilon'_{d}$ for a depolarizing surface). We note that the manner in which the diffuse dielectric constant is derived from the polarization dependence of Fresnel's equations is analogous to the method used by Janssen et al. (2009, 2016) for passive polarized radiometry observations of Titan.

### 3.4. Sensitivity of the Models

Because the incidence angle coverage of the observed backscatter curves is nonuniform, we tested the sensitivity of the models to each parameter for all incidence angles, following Lucas et al. (2019). This analysis was performed using the method developed by Sobol (2001), as implemented in Matlab by Cannavó (2012). The parameters are allowed to vary within the ranges found by the fits, including $1\sigma$ uncertainties (given in Section 4). The sensitivity index is the fraction of the total output variance attributed to each parameter (Sobol 2001; Lucas et al. 2019). The resulting sensitivity plots are shown for each model in Figure 4.

When using the Swift diffuse model, the polarization is an additional parameter. At $180°$ (parallel, or VV polarization), the transmissivity is almost the same at and above the Brewster angles, for all considered dielectric constants. Consequently, if parallel polarized data alone are considered, this method cannot be used to constrain the dielectric constant $\epsilon'_{d}$. If data in perpendicular or intermediate polarizations are included (as is
the case here; see Table 1), it becomes possible to constrain the dielectric constant from the Swift model, as illustrated in Figure 4. Although the sensitivity analysis does not illustrate it, data acquired at similar incidence but different polarization angles would give even better constraints on the dielectric constants.

From Figure 4, we can see that the dominant parameters controlling the shape of the fit are the rms slope for very low incidence angles (below ≈10°) due to the specular character of flat surfaces and the scattering albedo for all other incidence angles. This result is the same as found in the sensitivity analysis of Lucas et al. (2019). In practice, most regions near Selk present little to no low incidence angle data, and the rms slope is often poorly constrained. The scattering albedo is well constrained, especially by data at ≈30°–50° incidence angles. A known scattering albedo A or η allows the other two parameters ($\epsilon'_{qs}$ and $\epsilon'_{d}$ or n, depending on the diffuse model used) to be reasonably well constrained as well, depending on the available data at <30° and >60°. Finally, we also ran the fitting method on reduced data sets by cutting off all data below 15° incidence. Consistently with the sensitivity analysis, we found that all parameters derived were within error bars of those derived with the complete data set, but the rms slopes had such large uncertainties that they were not constrained. Thus, the incidence angle distribution of available data is well reflected in the uncertainties on the results.

# 4. Parameter Constraints from Fitting the Backscatter Curves

## 4.1. Parameter Constraints

We applied the models and fitting methods detailed in Section 3 to all six terrain units described in Section 2. The resulting fits are shown in Figure 5. The values extracted for all four parameters over the six terrains and using all six models are provided in Table 2 and shown (for three of the four parameters) in Figure 6. Table 2 also gives the reduced $\chi^2$ values $\chi_r^2$, which provides an indication of the quality of the fit. The fact that the $\chi_r^2$ values are similar for each terrain regardless of the quasi-specular model used underlines the difficulty of choosing a single model over the others: it is not clear which one fits the data best. The fits are often slightly better for the Acos$^n$ diffuse model than for the Swift model; however, because the Acos$^n$ model is empirical and the Swift model is physical, we prefer to keep both.

First, it is important to note that in some cases certain parameters are very poorly constrained owing to insufficient high- or low-incidence data. This is the case of the rms slope for the crater rim and the dark terrains, both of which do not feature any data below 15° incidence. Meanwhile, the dark terrains, dunes, and interdunes have little to no data beyond 45° incidence, and both the cosine exponent of the Acos$^n$ diffuse model and the dielectric constant $\epsilon'_{d}$ of the Swift model are poorly constrained, as illustrated both by very large uncertainties and by the sensitivity analysis in Section 3.4. Other parameters in several models are poorly constrained owing to insufficient and nonuniform data; this is reflected in the uncertainties.

As visible in Figure 6, we find that even when they are well constrained, all parameters are model dependent. The rms slope is always lowest for the Gaussian quasi-specular model and highest for the exponential model. This effect has already been

observed and analyzed by, e.g., Wye (2011), and also affects the dielectric constants, although it is less apparent owing to large uncertainties. The scattering albedo is expected to be higher for the Swift model than for the Acos$^n$ model, in order to compensate for the transmissivity term (which is <1); this is indeed observed in the extracted A and η values. In spite of the dependency of the parameters on the models, the relative values of most parameters are consistent across all models and allow us to probe into the structure and composition of Titan's surface.

## 4.2. Effective Dielectric Constants from the Quasi-specular and Diffuse Models

Dielectric constants derived from both the quasi-specular and diffuse (Swift) models are effective dielectric constants. In the presence of porosity, they represent a mix between the physical dielectric constant of the material and that of Titan's atmosphere (≈1) in the void space of pores. For many granular materials, this mixture is roughly linear (Hashin & Shtrikman 1962; Mätzler 1996; Sihvola 2000). As a result, the physical dielectric constant should be similar to the effective dielectric constant divided by one minus the average porosity. For a porosity of 30%, consistent with most granular materials on Earth, the physical dielectric constant of the material averaged throughout the footprint should be about 1.4× larger than the derived effective dielectric constants.

The median values of the effective dielectric constant derived from the Swift diffuse scattering model, $\epsilon'_{d}$, are consistently very low, from 1 to 1.97 (always <2). If confirmed, that is, if the Swift model of a uniform, low-roughness, single-scattering surface is representative of Titan's surface, such low dielectric constants would point to a very organic and/or porous surface everywhere; powdered tholins, for instance, have a dielectric constant of 1.17 (Paillou et al. 2008). Meanwhile, the quasi-specular model yields dielectric constants of around 2 for the dune fields and dark terrains but of 3–4 in the other terrains. This discrepancy between the values derived from the quasi-specular and diffuse models implies that at least one of these models does not accurately represent the surface.

In the quasi-specular models, the dielectric constant intervenes in the reflectivity term $\Gamma(\epsilon'_{qs})$ (Equation (5)), which acts as a backscattering amplitude independent of incidence angle (all incidence angle dependence of these models is captured in the rms tilts of the facets). However, not only do these models assume certain height distribution and autocorrelation functions that may not be representative of the surface (e.g., Shepard et al. 2001; Labarre et al. 2017), but they also only consider surface scattering. Indeed, if the medium is transparent enough, subsurface scattering on layers, voids, or cracks can be significant even at low incidence angles, an effect that is not accounted for in either the quasi-specular or diffuse models used herein. Given that the signal is expected to penetrate up to several meters into Titan's subsurface (for both water ice and solid hydrocarbons; Paillou et al. 2008), it is possible that the quasi-specular models considered do not accurately represent Titan's surface behavior.

By attempting to fit Fresnel's equations to the shape of the high incidence angle data at various polarization angles, the diffuse $\epsilon'_{d}$ values primarily express the degree of polarization of the surface. Thus, the low dielectric constants derived in this manner expose the depolarizing properties of Titan's surface. The method used here is analogous to the one applied to Cassini passive microwave
<figure><img src="image_8.png" /><figcaption>**Figure 5.** Model fits using the Hagfors + Swift model. The dots are the mean $\sigma^0$ values within each grid cell, and the black crosses are the predicted values for each data point, at the observation's incidence and polarization angles. The shaded red and blue regions represent the 95% confidence interval of the fitted curve for VV (parallel) and HH (perpendicular) polarizations, respectively. The best-fitting parameters (median +/- 1σ) are given for each fit; these are the same values as in Table 2. The bottom right panel presents the data and 95% confidence intervals for the modeled backscatter curves for the six mapped regions (excluding pure dunes and interdunes) in a single plot for easier comparison; note that the crater ejecta and hummocky terrains are almost perfectly superimposed.</figcaption></figure>

(The complete figure set (6 images), providing all model combinations, is available in the online Journal.)

radiometry observations of Titan (Janssen et al. 2009, 2016) and Rhea (Bonnefoy et al. 2020). Janssen et al. (2016) derived a global map of Titan's effective dielectric constant by comparing the passive thermal microwave emission at two orthogonal polarization angles (measured by rotating the spacecraft 90° in between two scans of Titan) to the brightness temperatures predicted by the White & Cogdell (1973) model. They found

dielectric constants of the same order as derived herein, with a global average of about 1.6, going up to 2.3, which they interpreted as a lower bound. For an average porosity of 30%, these values correspond to dielectric constants of 2.3–3.3 for Titan’s bulk physical material. It should be noted, however, that depolarization of a radar signal can occur for a number of reasons, including (1) multiple subsurface scattering, as found
<table><caption>**Table 2**
Parameter Constraints for All Six Combinations of Quasi-specular + Diffuse Model, for All Six Terrain Units</caption><thead><tr><th>Unit</th><th>Model</th><th>χ<sup>2</sup></th><th>ε'<sub>qs</sub></th><th>φ<sub>rms</sub></th><th>η or A</th><th>ε'<sub>d</sub> or n</th></tr></thead><tbody><tr><td rowspan="6">Crater rim</td><td>Hagfors + Swift</td><td>2.6</td><td>3.46<sup>+0.98</sup><sub>-0.85</sub></td><td>21.98<sup>+7.34</sup><sub>-8.21</sub></td><td>1.69<sup>+0.21</sup><sub>-0.32</sub></td><td>1.01<sup>+0.17</sup><sub>-0</sub></td></tr><tr><td>Exponential + Swift</td><td>2.6</td><td>4.13<sup>+0.57</sup><sub>-1.06</sub></td><td>29.34<sup>+5.5</sup><sub>-7.2</sub></td><td>1.89<sup>+0.11</sup><sub>-0.12</sub></td><td>1.01<sup>+0</sup><sub>-0</sub></td></tr><tr><td>Gaussian + Swift</td><td>2.4</td><td>3.98<sup>+0.7</sup><sub>-0.91</sub></td><td>19.92<sup>+2.51</sup><sub>-4.27</sub></td><td>1.89<sup>+0.12</sup><sub>-0.11</sub></td><td>1.01<sup>+0</sup><sub>-0</sub></td></tr><tr><td>Hagfors + Acos<sup>n</sup></td><td>2.6</td><td>3.61<sup>+0.99</sup><sub>-1.01</sub></td><td>24.27<sup>+7.64</sup><sub>-10.34</sub></td><td>0.84<sup>+0.14</sup><sub>-0.16</sub></td><td>1.05<sup>+0.22</sup><sub>-0.16</sub></td></tr><tr><td>Exponential + Acos<sup>n</sup></td><td>2.5</td><td>3.81<sup>+0.81</sup><sub>-1.46</sub></td><td>23.11<sup>+7.4</sup><sub>-7.43</sub></td><td>1<sup>+0.12</sup><sub>-0.11</sub></td><td>1.02<sup>+0.15</sup><sub>-0.14</sub></td></tr><tr><td>Gaussian + Acos<sup>n</sup></td><td>2.3</td><td>3.85<sup>+0.77</sup><sub>-0.92</sub></td><td>17.14<sup>+3.64</sup><sub>-4.84</sub></td><td>0.95<sup>+0.1</sup><sub>-0.09</sub></td><td>0.96<sup>+0.13</sup><sub>-0.13</sub></td></tr><tr><td rowspan="6">Crater ejecta</td><td>Hagfors + Swift</td><td>5.7</td><td>3.19<sup>+1.08</sup><sub>-0.95</sub></td><td>20.42<sup>+3.65</sup><sub>-4.73</sub></td><td>0.89<sup>+0.33</sup><sub>-0.31</sub></td><td>1.75<sup>+0.37</sup><sub>-0.45</sub></td></tr><tr><td>Exponential + Swift</td><td>5.8</td><td>3.31<sup>+1.15</sup><sub>-0.92</sub></td><td>24.47<sup>+5.02</sup><sub>-6.45</sub></td><td>1.22<sup>+0.29</sup><sub>-0.28</sub></td><td>1.57<sup>+0.78</sup><sub>-0.56</sub></td></tr><tr><td>Gaussian + Swift</td><td>5.4</td><td>3.79<sup>+0.77</sup><sub>-0.93</sub></td><td>15.77<sup>+2.65</sup><sub>-2.49</sub></td><td>1.08<sup>+0.25</sup><sub>-0.16</sub></td><td>1.29<sup>+0.54</sup><sub>-0.28</sub></td></tr><tr><td>Hagfors + Acos<sup>n</sup></td><td>5.1</td><td>3.36<sup>+1.04</sup><sub>-1.06</sub></td><td>21.89<sup>+4.02</sup><sub>-5.02</sub></td><td>0.45<sup>+0.17</sup><sub>-0.15</sub></td><td>1.43<sup>+0.3</sup><sub>-0.24</sub></td></tr><tr><td>Exponential + Acos<sup>n</sup></td><td>5.5</td><td>2.81<sup>+0.78</sup><sub>-0.78</sub></td><td>22.05<sup>+5.06</sup><sub>-7.14</sub></td><td>0.66<sup>+0.15</sup><sub>-0.12</sub></td><td>1.33<sup>+0.16</sup><sub>-0.19</sub></td></tr><tr><td>Gaussian + Acos<sup>n</sup></td><td>5.2</td><td>3.56<sup>+0.93</sup><sub>-0.84</sub></td><td>15.36<sup>+3.05</sup><sub>-2.4</sub></td><td>0.57<sup>+0.11</sup><sub>-0.09</sub></td><td>1.19<sup>+0.18</sup><sub>-0.19</sub></td></tr><tr><td rowspan="6">Hummocky terrains</td><td>Hagfors + Swift</td><td>4.6</td><td>2.92<sup>+1.11</sup><sub>-0.79</sub></td><td>17.75<sup>+3.85</sup><sub>-5.7</sub></td><td>1.17<sup>+0.35</sup><sub>-0.46</sub></td><td>2.08<sup>+0.96</sup><sub>-0.97</sub></td></tr><tr><td>Exponential + Swift</td><td>4.7</td><td>3.37<sup>+1.03</sup><sub>-0.92</sub></td><td>22.51<sup>+6.38</sup><sub>-5.58</sub></td><td>1.5<sup>+0.2</sup><sub>-0.21</sub></td><td>2.34<sup>+0.97</sup><sub>-0.56</sub></td></tr><tr><td>Gaussian + Swift</td><td>4.4</td><td>3.4<sup>+0.89</sup><sub>-0.73</sub></td><td>13.14<sup>+2.9</sup><sub>-2.35</sub></td><td>1.43<sup>+0.22</sup><sub>-0.32</sub></td><td>2.14<sup>+1.22</sup><sub>-0.93</sub></td></tr><tr><td>Hagfors + Acos<sup>n</sup></td><td>3.3</td><td>2.82<sup>+1.33</sup><sub>-0.92</sub></td><td>17.95<sup>+4.7</sup><sub>-6.3</sub></td><td>0.57<sup>+0.23</sup><sub>-0.23</sub></td><td>1.45<sup>+0.3</sup><sub>-0.32</sub></td></tr><tr><td>Exponential + Acos<sup>n</sup></td><td>3.5</td><td>2.68<sup>+0.99</sup><sub>-0.76</sub></td><td>17.94<sup>+6.01</sup><sub>-7.29</sub></td><td>0.78<sup>+0.15</sup><sub>-0.12</sub></td><td>1.51<sup>+0.24</sup><sub>-0.24</sub></td></tr><tr><td>Gaussian + Acos<sup>n</sup></td><td>3.4</td><td>3.1<sup>+1.13</sup><sub>-1.02</sub></td><td>12.78<sup>+2.82</sup><sub>-2.94</sub></td><td>0.71<sup>+0.16</sup><sub>-0.15</sub></td><td>1.38<sup>+0.28</sup><sub>-0.27</sub></td></tr><tr><td rowspan="6">Plains</td><td>Hagfors + Swift</td><td>8.2</td><td>3<sup>+0.67</sup><sub>-0.51</sub></td><td>20.43<sup>+2.88</sup><sub>-3.19</sub></td><td>0.27<sup>+0.16</sup><sub>-0.18</sub></td><td>1.78<sup>+0.34</sup><sub>-0.32</sub></td></tr><tr><td>Exponential + Swift</td><td>6.7</td><td>3.91<sup>+0.74</sup><sub>-0.85</sub></td><td>28.46<sup>+4.11</sup><sub>-4.76</sub></td><td>0.43<sup>+0.19</sup><sub>-0.09</sub></td><td>1.01<sup>+0.92</sup><sub>-0</sub></td></tr><tr><td>Gaussian + Swift</td><td>7.7</td><td>3.07<sup>+0.44</sup><sub>-0.36</sub></td><td>14.46<sup>+2.03</sup><sub>-1.53</sub></td><td>0.56<sup>+0.09</sup><sub>-0.09</sub></td><td>1.58<sup>+0.33</sup><sub>-0.44</sub></td></tr><tr><td>Hagfors + Acos<sup>n</sup></td><td>6.9</td><td>2.64<sup>+0.58</sup><sub>-0.55</sub></td><td>19.71<sup>+3</sup><sub>-3</sub></td><td>0.22<sup>+0.1</sup><sub>-0.11</sub></td><td>1.63<sup>+0.34</sup><sub>-0.29</sub></td></tr><tr><td>Exponential + Acos<sup>n</sup></td><td>6.6</td><td>2.98<sup>+0.72</sup><sub>-0.72</sub></td><td>23.41<sup>+4.74</sup><sub>-5.36</sub></td><td>0.3<sup>+0.08</sup><sub>-0.1</sub></td><td>1.35<sup>+0.27</sup><sub>-0.38</sub></td></tr><tr><td>Gaussian + Acos<sup>n</sup></td><td>6.9</td><td>2.92<sup>+0.55</sup><sub>-0.49</sub></td><td>14.05<sup>+2.23</sup><sub>-1.86</sub></td><td>0.3<sup>+0.06</sup><sub>-0.06</sub></td><td>1.32<sup>+0.23</sup><sub>-0.24</sub></td></tr><tr><td rowspan="6">Interdune regions</td><td>Hagfors + Swift</td><td>9.7</td><td>2.45<sup>+0.2</sup><sub>-0.31</sub></td><td>14.73<sup>+2.43</sup><sub>-2.53</sub></td><td>0.1<sup>+0.11</sup><sub>-0.07</sub></td><td>2.21<sup>+0.44</sup><sub>-0.36</sub></td></tr><tr><td>Exponential + Swift</td><td>9.8</td><td>3.38<sup>+0.44</sup><sub>-0.42</sub></td><td>23.79<sup>+3.58</sup><sub>-3.84</sub></td><td>0.2<sup>+0.07</sup><sub>-0.06</sub></td><td>2.2<sup>+0.46</sup><sub>-0.38</sub></td></tr><tr><td>Gaussian + Swift</td><td>10.7</td><td>3.12<sup>+0.36</sup><sub>-0.36</sub></td><td>12.68<sup>+1.56</sup><sub>-1.08</sub></td><td>0.21<sup>+0.07</sup><sub>-0.06</sub></td><td>2.17<sup>+0.48</sup><sub>-0.38</sub></td></tr><tr><td>Hagfors + Acos<sup>n</sup></td><td>10.1</td><td>2.53<sup>+0.16</sup><sub>-0.27</sub></td><td>15.48<sup>+2.3</sup><sub>-2.63</sub></td><td>0.02<sup>+0.07</sup><sub>-0.01</sub></td><td>0<sup>+1.81</sup><sub>-0</sub></td></tr><tr><td>Exponential + Acos<sup>n</sup></td><td>10.0</td><td>3.81<sup>+0.49</sup><sub>-0.7</sub></td><td>25.14<sup>+4.83</sup><sub>-4.81</sub></td><td>0.04<sup>+0.07</sup><sub>-0.02</sub></td><td>0.01<sup>+1.64</sup><sub>-0.01</sub></td></tr><tr><td>Gaussian + Acos<sup>n</sup></td><td>11.4</td><td>3.37<sup>+0.32</sup><sub>-0.5</sub></td><td>13.29<sup>+1.53</sup><sub>-1.58</sub></td><td>0.05<sup>+0.08</sup><sub>-0.01</sub></td><td>0.02<sup>+1.73</sup><sub>-0.02</sub></td></tr><tr><td rowspan="6">Dune fields</td><td>Hagfors + Swift</td><td>5.9</td><td>1.75<sup>+0.12</sup><sub>-0.12</sub></td><td>10.9<sup>+1.14</sup><sub>-2.39</sub></td><td>0.14<sup>+0.04</sup><sub>-0.05</sub></td><td>1.98<sup>+0.52</sup><sub>-0.41</sub></td></tr><tr><td>Exponential + Swift</td><td>6.0</td><td>2.19<sup>+0.17</sup><sub>-0.14</sub></td><td>17.07<sup>+2.54</sup><sub>-3.01</sub></td><td>0.19<sup>+0.03</sup><sub>-0.04</sub></td><td>1.92<sup>+0.75</sup><sub>-0.48</sub></td></tr><tr><td>Gaussian + Swift</td><td>6.8</td><td>2.08<sup>+0.14</sup><sub>-0.16</sub></td><td>9.67<sup>+1.08</sup><sub>-1.68</sub></td><td>0.19<sup>+0.03</sup><sub>-0.03</sub></td><td>1.98<sup>+0.51</sup><sub>-0.34</sub></td></tr><tr><td>Hagfors + Acos<sup>n</sup></td><td>5.4</td><td>1.7<sup>+0.17</sup><sub>-0.13</sub></td><td>10.4<sup>+2.63</sup><sub>-3.46</sub></td><td>0.08<sup>+0.03</sup><sub>-0.04</sub></td><td>1.69<sup>+0.35</sup><sub>-0.45</sub></td></tr><tr><td>Exponential + Acos<sup>n</sup></td><td>5.4</td><td>2.15<sup>+0.23</sup><sub>-0.19</sub></td><td>16.45<sup>+3.5</sup><sub>-3.29</sub></td><td>0.1<sup>+0.03</sup><sub>-0.03</sub></td><td>1.53<sup>+0.36</sup><sub>-0.43</sub></td></tr><tr><td>Gaussian + Acos<sup>n</sup></td><td>6.2</td><td>2<sup>+0.2</sup><sub>-0.17</sub></td><td>9.3<sup>+1.24</sup><sub>-2.07</sub></td><td>0.11<sup>+0.02</sup><sub>-0.02</sub></td><td>1.65<sup>+0.29</sup><sub>-0.33</sub></td></tr><tr><td rowspan="6">Dunes</td><td>Hagfors + Swift</td><td>6.9</td><td>1.75<sup>+0.13</sup><sub>-0.13</sub></td><td>10.9<sup>+2.43</sup><sub>-3.13</sub></td><td>0.06<sup>+0.04</sup><sub>-0.03</sub></td><td>2.25<sup>+1.47</sup><sub>-0.34</sub></td></tr><tr><td>Exponential + Swift</td><td>7.2</td><td>2.24<sup>+0.2</sup><sub>-0.19</sub></td><td>17.56<sup>+3.14</sup><sub>-2.87</sub></td><td>0.09<sup>+0.03</sup><sub>-0.02</sub></td><td>2.22<sup>+0.92</sup><sub>-0.3</sub></td></tr><tr><td>Gaussian + Swift</td><td>8.6</td><td>2.23<sup>+0.2</sup><sub>-0.21</sub></td><td>10.8<sup>+0.96</sup><sub>-1.14</sub></td><td>0.09<sup>+0.04</sup><sub>-0.03</sub></td><td>2.25<sup>+0.74</sup><sub>-0.32</sub></td></tr><tr><td>Hagfors + Acos<sup>n</sup></td><td>7.8</td><td>1.77<sup>+0.1</sup><sub>-0.13</sub></td><td>11.24<sup>+2.34</sup><sub>-2.68</sub></td><td>0.02<sup>+0.03</sup><sub>-0.01</sub></td><td>0.03<sup>+0.07</sup><sub>-0.03</sub></td></tr><tr><td>Exponential + Acos<sup>n</sup></td><td>7.8</td><td>2.31<sup>+0.2</sup><sub>-0.25</sub></td><td>18.41<sup>+2.96</sup><sub>-3.6</sub></td><td>0.02<sup>+0.04</sup><sub>-0.01</sub></td><td>0.05<sup>+1.96</sup><sub>-0.05</sub></td></tr><tr><td>Gaussian + Acos<sup>n</sup></td><td>10.0</td><td>2.28<sup>+0.22</sup><sub>-0.25</sub></td><td>11.05<sup>+1.32</sup><sub>-1.26</sub></td><td>0.02<sup>+0.04</sup><sub>-0.01</sub></td><td>0.06<sup>+1.95</sup><sub>-0.06</sub></td></tr><tr><td rowspan="6">Dark terrains</td><td>Hagfors + Swift</td><td>4.0</td><td>1.55<sup>+0.09</sup><sub>-0.12</sub></td><td>13.99<sup>+7.18</sup><sub>-2.78</sub></td><td>0.03<sup>+0.03</sup><sub>-0.02</sub></td><td>1.77<sup>+0.98</sup><sub>-0.57</sub></td></tr><tr><td>Exponential + Swift</td><td>4.1</td><td>2.12<sup>+0.91</sup><sub>-0.31</sub></td><td>16.17<sup>+6.01</sup><sub>-5.01</sub></td><td>0.05<sup>+0.02</sup><sub>-0.02</sub></td><td>1.78<sup>+1.22</sup><sub>-0.62</sub></td></tr><tr><td>Gaussian + Swift</td><td>4.3</td><td>1.78<sup>+0.3</sup><sub>-0.17</sub></td><td>12.67<sup>+3.98</sup><sub>-3.06</sub></td><td>0.05<sup>+0.03</sup><sub>-0.02</sub></td><td>1.76<sup>+1.38</sup><sub>-0.59</sub></td></tr><tr><td>Hagfors + Acos<sup>n</sup></td><td>3.9</td><td>1.56<sup>+0.08</sup><sub>-0.09</sub></td><td>14.44<sup>+7.35</sup><sub>-2.91</sub></td><td>0.01<sup>+0.01</sup><sub>-0.01</sub></td><td>0<sup>+1.22</sup><sub>-0</sub></td></tr><tr><td>Exponential + Acos<sup>n</sup></td><td>3.8</td><td>2.12<sup>+0.81</sup><sub>-0.27</sub></td><td>17.08<sup>+6.39</sup><sub>-5.37</sub></td><td>0.02<sup>+0.01</sup><sub>-0.01</sub></td><td>0<sup>+1.14</sup><sub>-0</sub></td></tr><tr><td>Gaussian + Acos<sup>n</sup></td><td>4.2</td><td>1.79<sup>+0.35</sup><sub>-0.15</sub></td><td>13.35<sup>+4.51</sup><sub>-3.37</sub></td><td>0.02<sup>+0.01</sup><sub>-0.01</sub></td><td>0<sup>+0.92</sup><sub>-0</sub></td></tr></tbody></table>

Note. For each parameter, the median and $1\sigma$ uncertainties are provided, following the method described in Section 3.3. These values are also represented graphically in Figure 6. The last two columns correspond to the parameters A and n for the Acos$^n$ model and to the parameters $\eta$ and $\epsilon'_d$ for the Swift model.

on Titan by Janssen et al. (2011, 2016); (2) large-scale and wavelength-scale roughness (White & Cogdell 1973); and (3) a wavelength-scale density gradient, as suggested by Janssen et al. (2016). These interpretations, as well as a very organic and/or porous surface, are all reasonable scenarios for Titan and are probably each partially valid in different terrains. Due

to the little data available (little variety in polarization angle for each incidence angle), the uncertainties on $\epsilon'_d$ are too large to be able to differentiate between terrains. Future work should include the application of this method to the global low-resolution scatterometry data set and comparison with the dielectric constants derived from passive data.
<figure><img src="image_9.png" /><figcaption>**Figure 6.** Values of the dielectric constant, rms slope, and scattering albedo derived for all eight terrains and for all six combinations of quasi-specular and diffuse scattering models. The values shown in this plot are given in Table 2.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Figure 7.** Variations of the NRCS $\sigma^0$ with incidence angle $\theta$ within dune fields and plains, both inside and outside of Selk crater. As described in Section 2.3.1, the $\sigma^0$ is averaged within each grid cell. Note that the backscattering behavior of both the dune fields and the plains units is consistent inside and outside the crater, pointing to likely similar composition and structure. As in Figure 3, the high $\sigma^0$ values around 40° incidence correspond to T121S6 and are likely caused by a calibration anomaly.</figcaption></figure>

# 5. Geological Interpretations

Given the results and model limits outlined above, we interpret the relative parameter evaluations in different terrains in order to derive new constraints for the surface properties and terrain evolution, both at the Selk impact site and in the surrounding plains and dune fields.

## 5.1. The Selk Impact Crater

### 5.1.1. The Presence of Crater Infill

The dune fields and plains located within Selk crater were initially mapped separately from those outside the crater, resulting in the backscatter curves shown in Figure 7. From this figure, we see that the backscattering behavior of the terrains

inside the crater is perfectly consistent with that outside the crater and helps complete the backscatter curves by adding data near 45° and 70°. At the resolutions of the SAR (hundreds of meters to several kilometers), these units are equivalent inside and outside the crater, indicating similar compositions and structure. The crater floor has likely been filled with windblown organic sand, organic material infalling from the atmosphere, and alluvial material eroded from the rim, which may be similar to material eroded from mountains and hummocky terrains near plains.

### 5.1.2. The Selk Crater Rim

The Selk crater rim is the brightest terrain in the ROI, resulting in the highest scattering albedo. With a $\sigma^0$ value near
−1 dB at 45° incidence, the crater rim is brighter than most of Xanadu and the inselbergs among the dune fields (Wye et al. 2007; Janssen et al. 2011; Lucas et al. 2019). It may, however, be similar to small regions within Xanadu, or other very bright terrains on Titan whose backscatter curves have not yet been examined. Within the crater rim region, a quasi-specular component is not actually necessary to provide a good fit to the data; this highlights the importance of diffuse scattering in the crater rim but may also simply be due to the unavailability of low-incidence data.

The high amplitude of diffuse scattering can be explained with one or several of the following surface properties: (1) very high centimeter-scale roughness; (2) multiple subsurface scattering on organized structures, such as cracks, voids, or embedded rocks; (3) a very low loss medium such as porous water ice, which would aid subsurface scattering; and (4) high local slopes, which would cause quasi-specular reflection at many incidence angles and globally increase the backscattering efficiency. Although these hypotheses leave many possibilities regarding the properties and formation of the surface, they indicate the presence of roughness at centimeter to kilometer scales, and/or an icy, cracked material. Both of these are generally consistent with a crater rim or proximal ejecta, with an unknown degree of erosion both modifying the roughness and bringing non-icy organic material. The presence of water ice appears likely, as it has both a higher real component of the dielectric constant than tholins and a low-loss tangent, allowing for longer path lengths in the subsurface and significant volume scattering. Thus, the results in the crater rim are consistent with the impact excavating water ice from the subsurface, while also creating a rough, heterogeneous, and fractured surface, a geologically plausible scenario.

### 5.1.3. The Selk Crater Ejecta: Similar to Hummocky Terrains, Brighter Than Plains

As apparent in Figures 3 and 6, the hummocky terrains and the Selk crater ejecta have similar backscatter curves, and all parameters are generally within $1\sigma$ uncertainties for both regions. These two units indeed look similar in the SAR images and were mapped separately only because of geological context: the crater ejecta surrounds the crater rim, whereas the hummocky terrains are isolated within dune fields or plains. It is possible that they are actually the same, for instance, if the impact took place within a large hummocky terrain, and the unit mapped as “ejecta” is actually preexisting unmodified terrain. An alternate explanation would be that the crater exposed the same substrate present throughout the region (and possibly exposed in the hummocky terrains), which would have been progressively eroded to look identical (to the radar) to the older hummocky terrain.

Within the plains, only the scattering amplitude ($A$ or $\eta$) is lower than in the ejecta and hummocky terrains; both the dielectric constant $\epsilon'_{qs}$ and the rms slope $\phi_{rms}$ are within $1\sigma$ uncertainties (Figure 6, Table 2). This is surprising, as the plains are expected to be smoother (lower $\phi_{rms}$) and less icy (lower $\epsilon'_{qs}$) than the hummocky terrains and crater ejecta. Although it is possible that the data are simply insufficient to separate these parameters, all three of these regions have reasonably good incidence angle coverage, including data at 7°–8°. Values of $\epsilon'_{qs}$ and $\phi_{rms}$ in plains, ejecta, and hummocky terrains point to similar surface compositions and roughness, whereas the higher diffuse scattering of the ejecta and hummocky terrain implies more

subsurface heterogeneity. These results would be consistent with a uniform superficial layer covering all three terrains. Another possibility is that, although the bulk composition is equivalent in all three units, the surface and subsurface roughness within the plains is at a smaller scale ($s = \tan \phi_{rms} \propto \xi/\zeta$), leading to higher diffuse scattering. Decimeter-scale rounded rocks, for example, as seen on the Huygens landing site and likely present within dry channels (Le Gall et al. 2010), would be very bright to the Cassini RADAR and could be present on the plains.

The effective dielectric constants $\epsilon'_{qs}$ derived in these three regions, though poorly constrained, are very high, with medians around 3–4. Unless substantial amounts of ammonia are present in the water ice, such high values are not realistic for Titan's surface, which is expected to be composed mainly of water ice ($\epsilon' = 3.13$) and hydrocarbons ($\epsilon' < 2.4$) with some degree of porosity (see Section 4.2). Nonetheless, the fact that the dielectric constant measured in these three regions is higher than in the dune fields and dark terrains seems to indicate a relative enrichment in water ice and/or a lower-porosity surface. This interpretation is also consistent with the higher scattering albedos (relative to the dune fields and dark terrains), as subsurface scattering occurs more easily in water ice than in hydrocarbons.

### 5.2. Dune Fields near the Selk Impact Crater

The dune fields of Titan cover about 17% of Titan's surface (Rodriguez et al. 2014; MacKenzie et al. 2021) and change little over hundreds of kilometers. As a consequence, their backscatter curve has already been studied extensively, on a global scale, in different regions separately, and separating dunes from interdune regions (Wye et al. 2007; Le Gall et al. 2011, 2014; Paillou et al. 2014; Lucas et al. 2019). Here we examined the backscatter curve of the dune fields as a single terrain, but also separately as dunes and interdunes, similar to Lucas et al. (2019).

From the parameters derived (Table 2 and Figure 6), the dunes, dune fields, and dark terrains have similar properties, whereas the interdune regions appear different. The fact that bulk dune fields are similar to individual dunes indicates that the signal from sandy dunes dominates, in spite of the presence of radar-bright interdune regions. This can be explained by the apparent abundance of sand around Selk crater (especially to the south), as evidenced by very thin interdunes and the proximity of the dark terrain unit. We note that, though the data north of Selk crater seem to align with the rest of the dune field data, it is practically only observed during one flyby (T95): our dune field's backscatter curve really represents the dune fields located to the south and directly around the crater.

The backscatter from the dark terrains is so low, at almost the noise level, that it is similar to the lakes and seas. Nonetheless, we deem liquid hydrocarbons or damp solid hydrocarbons (similar to the Huygens landing site; Atkinson et al. 2010) to be an unlikely interpretation for these dark terrains because (a) dune fields always surround these regions, pointing to an arid climate; (b) methane/ethane lakes have a smooth and uniform backscatter from the detection of the lake floor, whereas here dune-like structures are sometimes distinguished within the dark unit; and (c) although precipitation in the equatorial regions may lead to rivers, small lakes and pools, and a damp surface (Turtle et al. 2011), the presence of such large and persistent liquid bodies at these latitudes is unlikely (e.g., Lora et al. 2015). Instead, our preferred interpretation for the dark terrains is a thick sand sea, with few if
any distinguishable bedforms at the SAR resolution and with sand-filled interdune regions. A homogeneous covering of sand-sized organic particles would lead to both the low observed dielectric constant (median values of 1.55–2.12, although it is sometimes poorly constrained) and the very low scattering albedo. Possibly due to a topographic low or a convergence of winds, this spot acts as a sand sink. Either the dunes do not follow their usual (for Titan) linear morphology, or the interdunes are filled with enough sand for them to be indistinguishable to the radar. Given that even hydrocarbons are still a relatively low loss medium with several meters of penetration expected (Paillou et al. 2008), this indicates a covering of over a meter of sand.

### 5.2.1. Dielectric Constant $ε'_{qs}$

The Fresnel equations, used in all three quasi-specular models considered herein and in the Swift diffuse model, assume a smooth, homogeneous, and single-scattering surface. This appears to be valid over the sandy regions: dunes, dark terrains, and bulk dune fields if the dunes dominate. As a consequence, the absolute values of the dielectric constants may be examined and interpreted in these units.

For these three sandy units, the effective dielectric constant we derive is low for all models (median values of $ε'_{qs} = 1.55 - 2.31$), supporting the hypothesis of an organic sand composition. These values are similar to $ε' = 1.95$, derived using the Hagfors model over dune fields up to T16 (Wye et al. 2007), and $ε' \sim 2$, derived with the Gaussian model (or GOM) over pure dunes (Lucas et al. 2019), but slightly lower than $ε' = 2.4 \pm 0.4$ using simultaneously the Gaussian and exponential models over global dune fields (Le Gall et al. 2011). Dielectric constants of $ε' = 1.5 \pm 0.2$ were derived using polarimetry from passive radiometry over dune fields (Le Gall et al. 2011; Janssen et al. 2016), and $ε' = 1.7 \pm 0.2$ was calculated from emissivity values assuming a Kirchhoff surface for dark (sand-filled) dune fields by Le Gall et al. (2011) and for pure dunes by Bonnefoy et al. (2016). These values are all consistent with one another and point to an organic composition of moderate porosity (20%–50%, consistent with Earth dunes), with little to no water ice. Meanwhile, in the interdune regions the dielectric constants are higher, for some models as high as in the plains. This is consistent with the wide range of values found in the interdune regions by Bonnefoy et al. (2016) and Lucas et al. (2019) and with their interpretation of a distinct interdune composition featuring some icy materials and surface roughness (e.g., icy rocks). As in other regions of Titan, very high dielectric constants (>3) indicate a likely failure of the model to accurately represent the surface.

The Permittivity, Waves and Altimetry-Mutual Impedance Probe/Huygens Atmospheric Structure Instrument (PWA-HASI) experiment aboard the Huygens probe measured a dielectric constant at 45 Hz of $2.5 \pm 0.3$, compatible with a majority of solid organics (50%–60%) with at least 10% water ice and 15% porosity, and maybe liquid methane or ethane within pores (Hamelin et al. 2016; Lethuillier et al. 2018). A similarly low dielectric constant in sandy regions (dunes, dune fields, and dark regions) is thus consistent with a primarily organic composition, in spite of the differing morphology of dunes and the Huygens landing site. However, the dielectric constant of tholins decreases with increasing frequency (Brouet et al. 2016; Lethuillier et al. 2018), preventing any direct comparison between values at 45 Hz and 13.78 GHz.

### 5.2.2. rms Slope $\phi_{rms}$

The rms slope is very model dependent but is always at its lowest for the dune fields and the dunes: the surface is smooth at scales near the radar wavelength. Low roughness is consistent with a smooth sandy surface, while still allowing for significant topography (dunes or superimposed bedforms), which would be much larger than the wavelength. The rms slope of the interdune regions is slightly higher than in the dunes but remains lower than in the brighter terrains, consistent with the presence of both rocks and sand, as found by Lucas et al. (2019).

### 5.2.3. Scattering Albedo $A$ or $\eta$

The scattering albedo is almost zero for the dunes and dark regions: there is little to no diffuse component in the backscatter curve. This would occur for a very smooth surface, with either very high loss (all the signal is absorbed) or a homogeneous subsurface (the signal travels into the medium but is never scattered back). The dark regions' and the dunes' very low scattering albedo points to uniform sand, whereas the slightly higher values in the dune fields are caused by diffuse scattering in the interdune regions, as previously suggested by Le Gall et al. (2011) and Paillou et al. (2014). Yet the scattering albedo of the interdune regions remains very low (<0.2), sometimes even lower than for the dune fields. There are several possible explanations for the low interdune diffuse scattering, including (1) imperfect mapping, which might include some dune signal in the interdune regions; (2) a small (<1 m) covering of sand over some interdune regions, as has been observed from infrared data in radar-dark interdunes by Rodriguez et al. (2014) and Bonnefoy et al. (2016) but seems inconsistent with the radar-bright interdunes seen here; and (3) a relatively homogeneous subsurface even in the interdune regions, which could then be composed, for example, of larger-grained sand or sedimentized dunes (a sort of hydrocarbon sandstone).

## 6. Rim Topography from Radarclinometry

Alternating bright/dark regions perpendicular to the radar look direction are apparent over the Selk crater rim (Figure 8(a)), a feature characteristic of the presence of topography and commonly associated with mountains and crater rims on Titan (Malaska et al. 2016b; Schoenfeld et al. 2021). Assuming uniform surface properties (composition, roughness, etc.) over a region, it is possible to use the observed backscatter to infer surface slopes and thus topography using radarclinometry. This technique has previously been applied to Titan's dune fields in regions where the interdune regions are radar-dark, revealing dune heights of 45–180 m (Neish et al. 2010), and to several of its mountains, which can reach heights of 3.3 km (Radebaugh et al. 2007; Mitri et al. 2010; Liu et al. 2011). Herein we use radarclinometry to derive the topography of the Selk crater rim.

### 6.1. Method

Of the five SAR images that include the rim (see Figure 1 and Table 1), only T98 is suitable for radarclinometry. Indeed, T95 has high speckle and a calibration problem in beam 3 over the rim, leading to lower $σ^0$ values in part of the data; T120 covers only a small section of the rim; and T36 and T121 both have a very poor resolution and a lot of noise. Meanwhile,
<figure><img src="image_11.png" /><figcaption>Figure 8. Radarclinometry example over the Selk crater rim from T98, showing the NRCS σ⁰ (left), the extracted slopes over the crater rim (middle), and the altitude (right). The slopes are calculated assuming only diffuse scattering from the Swift model, with ε'₉ = 1 in this case. Radar look angle is as indicated in the blue arrow.</figcaption></figure>

during the T98 flyby, RADAR imaged most of the crater for a long period of time, with the number of looks varying from 10 to more than 200, resulting in a very low level of speckle noise. The high number of looks also implies that the incidence angle and possibly the calibration varied slightly over the flyby, causing a slight heterogeneity of backscatter over the image (for example, the seam to the left of the blue arrow in Figure 8, or the very bright area in the southern part of the rim). The radarclinometry process adjusts for this anomaly by fitting the scattering amplitude independently along each profile, as described below.

For a given diffuse backscatter function f(θ), where θ is the incident angle, the backscatter intensity I in the SAR images is derived from geometry to be (Paquerault 1998; Le Hégarat-Mascle et al. 2005)

$$ 
\begin{aligned}
I(\theta, \alpha, \beta) &= C f(\theta - \alpha, \beta) \frac{\sin\theta \cos(\theta - \alpha) \cos\beta}{\sin(\theta - \alpha) \cos\beta} \\
&= C f(\theta - \alpha, \beta) \sin\theta \cot(\theta - \alpha),
\end{aligned}
\quad (10) $$

where C is a calibration constant and α and β are the local slopes in the directions parallel and perpendicular to the look direction, respectively (along and across the profile). The 1/ sin(θ - α) factor is necessary for the normalization of σ⁰ by the correct area, whereas the cos(θ - α) factor is a correction for the incident power over a sloped surface. At the high incidence angles of T98 (65°–70°), scattering is purely diffuse. Both the Swift and the Acosⁿ models described in Section 3.2 were tested for f(θ), with values of 1 < ε'₉ < 3 (Swift model) and 1 < n < 2 (Acosⁿ model). Along each profile in the look direction, separating the two sides of the rim if the profile crosses the crater, we start by solving Equation (10) for the calibration constant C assuming a flat surface (α = β = 0), while n or ε'₉ is kept fixed. This ensures that the start and end of the profile have the same elevation and accounts for calibration differences between different profiles. Note that the scattering amplitude A or η of the diffuse models is incorporated into C. Assuming that the dependence on β, the slope perpendicular to the look direction, is negligible, we can solve Equation (10) for the local incidence angle θ - α and the corresponding slope α within each pixel (shown in Figure 8(b)). The height differential Δh within a given pixel of resolution R can then be derived as follows:

$$ \Delta h = R \frac{\sin \theta \sin \alpha}{\sin(\theta - \alpha)}. \quad (11) $$

The above equation takes into account geometrical distortions due to topography. The horizontal distance (pixel size) along the profile is also distorted by the topography; however, at such high incidence angles this effect is negligible (about 1% of the resolution). The height is integrated along the look direction while taking the moving median of 5 pixels along the azimuth direction, providing the digital terrain model (DTM) shown in Figure 8(c). In order to reduce the striping effect caused by the propagation of speckle noise and pixel-scale mapping uncertainties, we repeat this procedure for azimuth angles varying between ±20° around the look direction and then take the median of the resulting values, as in Paquerault (1998). The same steps are followed for all values of n and ε'₉ with a step of 0.1.

## 6.2. Interpretation

The topography derived from radarclinometry has significant caveats and must always be used with caution. (1) Radarclinometry is most sensitive to slopes toward or away from the radar. On the southeast and northwest sides of the crater (grayed out in Figure 9), the dominant slopes are oriented orthogonally to the look direction: radarclinometry does not produce useful results in these regions. Therefore, although the slopes and heights in these areas were estimated and are shown in Figures 8 and 9, they should not be used and are not interpreted. On the other two sides of the rim (northeast and southwest), most slopes should face toward or away from the radar, and the topography has the expected shape across the rim. (2) The assumption of a single backscatter function along the profile is a simplification: the crater rim is likely to include valleys and alluvial fans, which may accumulate sand or gravel and lead to radar-dark regions easily confused with negative slopes. (3) The presence of noise and the uncertainties in the mapping of the rim's edges can lead to offsets in altitude that are propagated along the profile, leading to the stripes visible in Figure 8. (4) The range resolution of the SAR image is of 1.3–1.5 km pixel⁻¹, meaning that topography extending over shorter distances would be averaged out. Improving the resolution would likely lead to high rms slopes (Shepard et al. 2001) but would not change the rim height, as long as the ~20 km wide rim is well resolved (which is the case for the radarclinometry but not for SARTopo). Even though no slopes
<figure><img src="image_12.png" /><figcaption>**Figure 9.** The median radarclinometry-derived DTM obtained for different values of *n* or *ε*'<sub>d</sub> is averaged within each SARTopo footprint along the only SARTopo track that crosses the Selk crater rim (from north to south). The heights obtained from SARTopo and from the DTM are plotted together on the left for the Acos° model (top) and the Swift model (bottom), showing good agreement, especially for low values of *n* and *ε*'<sub>d</sub>. On the right, we show the best median DTM, corresponding to *n* = 1 (top) and *ε*'<sub>d</sub> = 1 (bottom). The sides of the rim parallel to the look direction are grayed out: because the radar is not oriented perpendicularly to most of the expected slopes there, the altitudes in these regions are less reliable. The SARTopo pass crossing the crater is outlined, and an example SARTopo footprint is drawn in red. Note that although the SARTopo footprints overlap, they are 8–10 km long.</figcaption></figure>

steeper than 16° are indicated in our radarclinometry results, we cannot rule out topography such as cliffs and deep valleys on scales smaller than the radar resolution. At 70° incidence, the shadow from a single mountain over flat topography would extend over a distance of about 3 times the mountain height in a radar image: the lack of obvious shadowing therefore only excludes single-standing, steep, kilometer-tall mountains.

To verify the validity of the radarclinometry DTM, we compared it to the SARTopo data over the rim, which includes three tracks along the eastern edge and one track through the middle of the crater. The radarclinometry DTM over the Selk crater rim has a resolution of about 1.3 km pixel⁻¹, whereas the SARTopo rim-crossing track has footprints of about 2.5 by 9 km. The radarclinometry DTM was averaged within each SARTopo footprint; the result is shown for the rim-crossing track in Figure 9. The match between radarclinometry and SARTopo is best for the lowest values of *n* (≈1) and of the dielectric constant *ε*'<sub>d</sub> (≈1), consistent with the results of backscatter curve fitting presented in Section 4 and with a very diffusely scattering surface. The radarclinometry averaged in SARTopo footprints matches almost perfectly with the SARTopo heights (Figure 9), even though the radarclinometry altitudes along the track were often slightly higher (at most 1.5 times as high). The consistency between these two very different methods of deriving topography increases the credibility of both, and similar comparisons should be made where possible along SARTopo tracks on Titan.

The higher resolution and the wider coverage of the radarclinometry DTM permit the detection of topography variations invisible to SARTopo. With the lowest values of *n*

and *ε*'<sub>d</sub>, both diffuse models yield maximum heights of 600 m with respect to the edge of the rim, ignoring the unreliable data along the edges (Figure 9). These heights, which remain consistent with SARTopo when averaged over the footprint, indicate that the rim can locally be at least over twice as high as SARTopo's rim-to-surrounding terrains heights of 280 m (Hedgepeth et al. 2020). While the altitudes derived remain coherent with those found in other craters on Titan, the higher topography is significant for analyses of the formation, degradation, and age of Selk crater, which may be less eroded than previously thought (Neish et al. 2013; Hedgepeth et al. 2020; Lorenz et al. 2021). This topography would also affect regional winds, a process that may play a role in the variations in dune orientations to the north and south of the crater (Malaska et al. 2016b). While this result does not change global topography trends, it supports the possibility for local peaks in many mountainous regions of Titan, similar to those found by Radebaugh et al. (2007) and Liu et al. (2011) but invisible to SARTopo and radar altimetry.

## 7. Conclusion

The Selk crater region, where the Dragonfly mission is expected to land in the mid-2030s (Lorenz et al. 2021), has been mapped into six units, and the dunes and interdune regions were separated within dune fields. The backscatter curves for each of these terrains have been assembled from SAR data acquired during nine Cassini flybys of Titan and cover incidence angles from 5° to 72° and polarization angles varying between parallel and perpendicular to the incidence
plane. Six different models, composed of the sum of a quasi-specular and a diffuse scattering term, have been fit to these data to derive variations in dielectric constant, rms slope, and scattering albedo from one terrain to another. The interaction between topography and backscatter was also examined over the crater rim. The main results and interpretations are summarized below:

1. Dunes and plains exhibit the same microwave scattering properties both inside and outside the crater, indicating likely aeolian infilling and/or crater rim erosion bringing the same materials into the crater as are available elsewhere.

2. The crater rim is among the radar-brightest terrains on Titan and exhibits strong diffuse scattering, consistent with an icy (low-loss) subsurface with buried scattering structures, although surface roughness likely also plays a role.

3. Around Selk crater, the plains, hummocky terrains, and ejecta blanket all have similar rms slopes and effective dielectric constants, consistent with a similar composition and surface structure (perhaps due to a thin uniform cover). However, the large uncertainties on these parameters still allow for some variation in composition and structure.

4. The dune fields and more specifically the dunes within them have a low effective dielectric constant consistent with previous work and a low rms slope and show little diffuse backscatter. These properties all point to homogeneous organic sand. Meanwhile, the interdune regions have a higher dielectric constant and rms slope, pointing to a rougher and likely icier surface.

5. The dark regions, mapped as very radar-dark terrains with no clearly identifiable dunes, have a low dielectric constant and little to no diffuse scattering. This is consistent with organic sand over depths thicker than ~1 m and likely corresponds to a sediment sink due to converging winds or low topography.

6. The active radar data can be used to derive the effective dielectric constant not only from the shape of the quasi-specular component but also from high-incidence data at different polarization angles, similar to the method used on passive microwave radiometry (Janssen et al. 2009, 2016). Future work should apply this method globally to the low-resolution scatterometry data, which were acquired at a large variety of polarization angles. This method highlights the importance of polarimetry in radar studies and the need to integrate radar polarimetry in future missions, both for a possible Titan orbiter (Rodriguez et al. 2022) and for radar-focused missions planned toward other targets like Venus (VERITAS/NASA and EnVision/ESA).

7. The effective dielectric constants derived from the quasi-specular and diffuse models agree in presumably sandy regions. The smoothness and homogeneity of these terrains conform to the assumptions of the models, which use Fresnel's equations and ignore multiple surface or subsurface scattering. The dielectric constants between 1.5 and 2.3 derived in these "sandy" regions (dune fields, pure dunes, and radar-dark terrains) are therefore considered physical and match well with previous estimates, either globally or in other regions of Titan (Le Gall et al. 2011; Bonnefoy et al. 2016; Janssen et al. 2016; Lucas et al. 2019).

8. In the more radar-bright and likely less sandy regions, the quasi-specular model predicts very high dielectric constants

(3-4), much higher than the values <2 obtained from the Swift diffuse model. This inconsistency implies that at least one of these models is not representative of Titan's surface, especially in regions where volume scattering is significant since both components assume only single scattering. The diffuse model, which is sensitive to the degree of polarization, may give artificially low values when the signal is depolarized and the surface does not follow Fresnel's equations. The quasi-specular model, which does not account for subsurface penetration, may overestimate the dielectric constant to compensate. The permittivity probe DIEL within the DraGMET instrument aboard the Dragonfly mission will measure the permittivity at different locations of the Dragonfly landing site. Although the frequency of this instrument is much lower (10 Hz-1 kHz) than that of Cassini RADAR (13.78 GHz), it will provide a ground truth measurement. By comparing it to dielectric constants measured at different frequencies both in the lab (Brouet et al. 2016; Lethuillier et al. 2016, 2018) and on Titan, we will be able to verify the validity of our models within different terrains and apply them to Titan globally.

9. Radarclinometry over the Selk crater rim provides topography across the rim in good agreement with the SARTopo data set but reveals higher altitudes (up to around 600 m for the best-fitting model) owing to wider coverage and a finer resolution. These higher altitudes are significant both for the geology of the region (crater degradation, dune-forming winds) and for the local meteorology, hence providing useful new information for the preparation of the Dragonfly mission.

Funding support for this research was provided by a postdoctoral research grant awarded by the Centre National d'Études Spatiales (CNES) and by NASA's Astrobiology Institute's proposal "Habitability of Hydrocarbon Worlds: Titan and Beyond" (PI R. M. Lopes). Cassini data are available through NASA's Planetary Data System at the following URL: https://pds-imaging.jpl.nasa.gov/volumes/radar.html. The authors wish to thank the L'Oréal-UNESCO For Women in Science Program for their support, as well as an anonymous reviewer for their insightful comments.

## ORCID iDs

Léa E. Bonnefoy https://orcid.org/0000-0003-1842-7531

Antoine Lucas https://orcid.org/0000-0003-2192-4416

Alexander G. Hayes https://orcid.org/0000-0001-6397-2630

Sébastien Rodriguez https://orcid.org/0000-0003-1219-0641

Valerio Poggiali https://orcid.org/0000-0003-1053-3225

Daniel E. Lalich https://orcid.org/0000-0003-2091-6991

Ralph D. Lorenz https://orcid.org/0000-0001-8528-4644

Alice Le Gall https://orcid.org/0000-0002-9023-4868

## References

Ashcraft, I. S., & Long, D. G. 2006, JGlac, 52, 257

Atkinson, K. R., Zarnecki, J. C., Towner, M. C., et al. 2010, Icar, 210, 843

Barnes, J. W., Brown, R. H., Soderblom, L., et al. 2008, Icar, 195, 400

Bonnefoy, L. E., Hayes, A. G., Hayne, P. O., et al. 2016, Icar, 270, 222

Bonnefoy, L. E., Le Gall, A., Lellouch, E., et al. 2020, Icar, 352, 113947

Brouet, Y., Levasseur-Regourd, A. C., Sabouroux, P., et al. 2016, MNRAS, 462, S89

Cannavó, F. 2012, CG, 44, 52

Elachi, C., Allison, M. D., Borgarelli, L., et al. 2004, SSRv, 115, 71

Elachi, C., Wall, S., Allison, M., et al. 2005, Sci, 308, 970
Elachi, C., Wall, S., Janssen, M., et al. 2006, *Natur*, **441**, 709

Evans, J. V., & Pettengill, G. H. 1963, *JGR*, **68**, 423

Fung, A. K., Li, Z., & Chen, K. S. 1992, *ITGRS*, **30**, 356

Grings, F., Franco, M., Spagnuolo, M. G., Janssen, M. A., & Lorenz, R. 2021, *Icar*, **359**, 114319

Hagfors, T. 1964, *JGR*, **69**, 3779

Hagfors, T. 1966, *JGR*, **71**, 379

Hamelin, M., Lethuillier, A., Le Gall, A., et al. 2016, *Icar*, **270**, 272

Harmon, J. K., Campbell, D. B., & Ostro, S. J. 1982, *Icar*, **52**, 171

Hashin, Z., & Shtrikman, S. 1962, *JAP*, **33**, 3125

Hayes, A. G., Aharonson, O., Lunine, J. I., et al. 2011, *Icar*, **211**, 655

Hayes, A. G., Lorenz, R. D., & Lunine, J. I. 2018, *NatGe*, **11**, 306

Hedgepeth, J. E., Neish, C. D., Turtle, E. P., et al. 2020, *Icar*, **344**, 113664

Heiles, C. E., & Drake, F. D. 1963, *Icar*, **2**, 281

Hofgartner, J. D., Hayes, A. G., Lunine, J. I., et al. 2014, *NatGe*, **7**, 493

Janssen, M. A., Le Gall, A., Lopes, R. M. C., et al. 2016, *Icar*, **270**, 443

Janssen, M. A., Le Gall, A., & Wye, L. C. 2011, *Icar*, **212**, 321

Janssen, M. A., Lorenz, R. D., West, R., et al. 2009, *Icar*, **200**, 222

Labarre, S., Ferrari, C., & Jacquemoud, S. 2017, *Icar*, **290**, 63

Lalich, D. E., Poggiali, V., Hayes, A. G., et al. 2022, *Icar*, **374**, 114775

Le Gall, A., Janssen, M. A., Kirk, R. L., & Lorenz, R. D. 2014, *Icar*, **230**, 198

Le Gall, A., Janssen, M. A., Paillou, P., Lorenz, R. D., & Wall, S. D. 2010, *Icar*, **207**, 948

Le Gall, A., Janssen, M. A., Wye, L. C., et al. 2011, *Icar*, **213**, 608

Le Gall, A., West, R. D., & Bonnefoy, L. E. 2019, *GeoRL*, **46**, 11747

Le Hégarat-Mascle, S., Zribi, M., & Ribous, L. 2005, *IJRS*, **26**, 2877

Lethuillier, A., Le Gall, A., Hamelin, M., et al. 2016, *A&A*, **591**, A32

Lethuillier, A., Le Gall, A., Hamelin, M., et al. 2018, *JGRE*, **123**, 807

Liu, Z. Y. C., Radebaugh, J., Kirk, R. L., et al. 2011, *LPSC*, **42**, 2798

Lopes, R. M. C., Malaska, M. J., Solomonidou, A., et al. 2016, *Icar*, **270**, 162

Lopes, R. M. C., Wall, S. D., Elachi, C., et al. 2019, *SSRv*, **215**, 33

Lora, J. M., Lunine, J. I., & Russell, J. L. 2015, *Icar*, **250**, 516

Lorenz, R. D., MacKenzie, S. M., Neish, C. D., et al. 2021, *PSJ*, **2**, 24

Lorenz, R. D., Wall, S., Radebaugh, J., et al. 2006, *Sci*, **312**, 724

Lucas, A., Aharonson, O., Deledalle, C., et al. 2014a, *JGRE*, **119**, 2149

Lucas, A., Rodriguez, S., Lemonnier, F., et al. 2019, *JGRE*, **124**, 3140

Lucas, A., Rodriguez, S., Narteau, C., et al. 2014b, *GeoRL*, **41**, 6093

MacKenzie, S. M., Birch, S. P. D., Hörst, S., et al. 2021, *PSJ*, **2**, 112

Malaska, M. J., Lopes, R. M. C., Hayes, A. G., et al. 2016a, *Icar*, **270**, 183

Malaska, M. J., Lopes, R. M. C., Williams, D. A., et al. 2016b, *Icar*, **270**, 130

Mastrogiuseppe, M., Poggiali, V., Seu, R., Martufi, R., & Notarnicola, C. 2014, *Icar*, **230**, 191

Mätzler, C. 1996, *ITGRS*, **34**, 573

Michaelides, R. J., Hayes, A. G., Mastrogiuseppe, M., et al. 2016, *Icar*, **270**, 57

Mitri, G., Bland, M. T., Showman, A. P., et al. 2010, *JGRE*, **115**, E10002

Neish, C. D., Kirk, R. L., Lorenz, R. D., et al. 2013, *Icar*, **223**, 82

Neish, C. D., Lorenz, R. D., Kirk, R. L., & Wye, L. C. 2010, *Icar*, **208**, 385

Neish, C. D., Lorenz, R. D., Turtle, E. P., et al. 2018, *AsBio*, **18**, 571

Ostro, S. J., West, R. D., Janssen, M. A., et al. 2006, *Icar*, **183**, 479

Ostro, S. J., West, R. D., Wye, L. C., et al. 2010, *Icar*, **206**, 498

Paillou, P., Bernard, D., Radebaugh, J., et al. 2014, *Icar*, **230**, 208

Paillou, P., Crapeau, M., Elachi, C., Wall, S., & Encrenaz, P. 2006, *JGRE*, **111**, E06S11

Paillou, P., Lunine, J., Ruffié, G., et al. 2008, *GeoRL*, **35**, L18202

Paquerault, S. 1998, PhD thesis, Paris Inst. of Technology, https://pastel.archives-ouvertes.fr/tel-00005664

PDS Cartography and Imaging Sciences Node 2021, Cassini RADAR Data Archive, https://pds-atmospheres.nmsu.edu/data_and_services/atmospheres_data/Cassini/inst-radar.html

Poggiali, V., Hayes, A. G., Mastrogiuseppe, M., et al. 2020, *JGRE*, **125**, e06558

Poggiali, V., Mastrogiuseppe, M., Callegari, M., et al. 2012, *Proc. SPIE*, **8536**, 853608

Poggiali, V., Mastrogiuseppe, M., Hayes, A. G., et al. 2019, *ITGRS*, **57**, 7262

Radebaugh, J., Lorenz, R. D., Kirk, R. L., et al. 2007, *Icar*, **192**, 77

Remund, Q. P., & Long, D. G. 2003, *ITGRS*, **41**, 1821

Rodriguez, S. 2003, PhD thesis, Université Pierre et Marie Curie—Paris VI, https://tel.archives-ouvertes.fr/tel-00005551

Rodriguez, S., Garcia, A., Lucas, A., et al. 2014, *Icar*, **230**, 168

Rodriguez, S., Vinatier, S., Cordier, D., et al. 2022, *ExA*, [Online First]

Schoenfeld, A. M., Lopes, R. M., Malaska, M. J., et al. 2021, *Icar*, **366**, 114516

Seignovert, B., Le Mouélic, S., Brown, R. H., et al. 2019, *EPSC*, *EPSC-DPS2019-1063*

Shepard, M. K., & Campbell, B. A. 1999, *Icar*, **141**, 156

Shepard, M. K., Campbell, B. A., Bulmer, M. H., et al. 2001, *JGR*, **106**, 32777

Sihvola, A. 2000, *Subsurface Sensing Technologies and Applications*, **1**, 393

Simpson, R. A., Tyler, G. L., & Campbell, D. B. 1978, *Icar*, **33**, 102

Sobol, I. M. 2001, *Mathematics and Computers in Simulation*, **55**, 271

Soderblom, J. M., Brown, R. H., Soderblom, L. A., et al. 2010, *Icar*, **208**, 905

Stephen, H., & Long, D. G. 2005, *ITGRS*, **43**, 238

Stiles, B. W., Hensley, S., Gim, Y., et al. 2009, *Icar*, **202**, 584

Sultan-Salem, A. K., & Tyler, G. L. 2006, *JGRE*, **111**, E06S07

Sultan-Salem, A. K., & Tyler, G. L. 2007, *JGRE*, **112**, E05012

Swift, C. T. 1999, *ITGRS*, **37**, 716

Swift, C. T., Hayes, P. S., Herd, J. S., Jones, W. L., & Delnore, V. E. 1985, *JGR*, **90**, 1983

Tomasko, M. G., Archinal, B., Becker, T., et al. 2005, *Natur*, **438**, 765

Turtle, E. P., Perry, J. E., Barbara, J. M., et al. 2018, *GeoRL*, **45**, 5320

Turtle, E. P., Perry, J. E., Hayes, A. G., et al. 2011, *Sci*, **331**, 1414

Turtle, E. P., Trainer, M. G., Barnes, J. W., et al. 2019, *LPSC*, **50**, 2888

Tyler, G. L., Simpson, R. A., Maurer, M. J., & Holmann, E. 1992, *JGR*, **97**, 13115

Ulaby, F., & Long, D. 2015, Microwave Radar and Radiometric Remote Sensing (London: Artech House)

Werynski, A., Neish, C., Le Gall, A., & Janssen, M. 2019, *Icar*, **321**, 508

White, T. L., & Cogdell, J. R. 1973, *EM&P*, **6**, 235

Williams, D. A., Radebaugh, J., Lopes, R. M. C., & Stofan, E. 2011, *Icar*, **212**, 744

Wye, L. C. 2011, PhD thesis, Stanford Univ.

Wye, L. C., Zebker, H. A., & Lorenz, R. D. 2009, *GeoRL*, **36**, L16201

Wye, L. C., Zebker, H. A., Ostro, S. J., et al. 2007, *Icar*, **188**, 367

Zebker, H., Hayes, A., Janssen, M., et al. 2014, *GeoRL*, **41**, 308

Zebker, H. A., Wye, L. C., Janssen, M. A. & the Cassini Radar Team 2008, *Icar*, **194**, 704