Investigating the Martian Surface at Decametric Scale: Population, Distribution, and Dimension of Heterogeneity from Radar Statistics

The studied region covers longitudes 135°E to 165°E and latitudes −5°N to 15°N. We chose this region because it gathers in a contained area most of the 30 dB wide dynamic range of SHARAD reflectivity that can be observed planet-wide (Grima et al. 2012). The region includes a range of dim terrains including Medusae Fossae Formation (MFF) outliers and scattered Noachian rises (NR), the moderately bright Hesperian Elysium shield (ES), as well as the middle-to-late Amazonian volcanic complex of the Elysium plains (EP) that account for one of the brightest regions of Mars at SHARAD wavelengths. This broad continuum of reflectivity that reflects the heterogeneity of terrain textures is essential in making physical sense of the relative coherent and incoherent signal behavior described in Section 4 that frames the methodology we use to approximate surface roughness. The studied region also encompasses the Curiosity and Insight's landing sites.

SHARAD transmits a 10 MHz bandwidth centered at 20 MHz (λ = 15 m) (Seu et al. 2004; Croci et al. 2011), and has been in operation aboard NASA's Mars Reconnaissance Orbiter (MRO) since 2006. The SHARAD signal is reflected by each dielectric gradient on its propagation path until extinction. The resulting echoes are recorded at the antenna, providing the electric field amplitude and phase along the transmitted polarization axis. The transmitted signal is characterized by a nearly isotropic illumination pattern so that the first recorded echo is usually originating from the closest set of surface reflectors–scatterers, which, at first order, is assumed to be near the spacecraft nadir point. SHARAD provides an effective pulse-repetition interval every ∼20–40 m along the track.

We use the raw data products (also called Experimental Data Records, hereafter EDR) available on the Planetary Data System platform. ⁷ The data set represents a global coverage from more than 11 yr of observations at the time of processing. The usage of the raw products ensures that the data are free from ground-based post processing that could alter the radiometric gain through a terrain-dependent component hard to predict a priori. For instance, focusing techniques, commonly integrated in derived radar data products, have a radiometric stability dependent on the along-track surface backscattering function produced by roughness and local slopes (Peters et al. 2007; Scanlan et al. 2021). A SHARAD radargram can be of any length along the track, but has a fixed vertical size corresponding to a receiver window of 135 μs digitized over 3600 samples.

Preprocessing. The EDR data preprocessing is identical to what is described in more detail by Steinbrugge et al. (2022, Section II-B). It includes (i) a range compression using an ideal reference chirp following the methodology by Campbell et al. (2011), and (ii) a correction of the dispersive phase shift and the related signal loss due to the propagation through the ionosphere by using an approach maximizing the signal-to-noise ratio (S/N) to focus the surface echo. Absorption losses that occur during the ionosphere propagation are not corrected. Depending on radio solar flux conditions, such loss is contained below 1.5 dB at SHARAD frequencies (Grima et al. 2012). That is of the second order in comparison with the geologic signal dynamic range dealt with in the rest of this study. One-third of our data set has been acquired on the day side and might suffer such signal absorption.

Surface echo extraction. We extracted the surface echo strength from the preprocessed EDR using the detection criteria proposed by Grima et al. (2012). This criteria looks for the maximum of the stronger surface echo strength weighted by the derivative preceding it in the fast-time domain (i.e., the vertical dimension). However, as we use the range-compressed EDR data that have a low S/N, we avoid unreliable detection by computing the criteria over a restricted 100 sample long vertical window centered on a location where the surface is expected. This predetermined location is obtained through pulse retracking techniques adapted from terrestrial ocean altimetry to identify the rising edge of the surface return as applied by Steinbrugge et al. (2022) on SHARAD.

Gain corrections. MRO quasi-polar orbits are characterized by a periapsis (resp. apoapsis) reaching ∼250 km (resp. ∼320 km) above the Martian geoid over the north (resp. south) polar regions. The total variation in gain due to geometric losses reaches 3 dB (Grima et al. 2012; Campbell et al. 2021). We adjusted each echo strength using a relative geometric propagation factor of ${f}^{2}={(h/{h}_{\mathrm{ref}})}^{2}$ where h represents the observed altitude, and h_ref is a common reference arbitrarily set to 250 km (Grima et al. 2012). The configuration of MRO's solar arrays (SA) and high-gain antenna (HGA) alters the radiation pattern of the SHARAD antenna and is therefore a major source of signal-to-noise variations. We correct for this effect by applying the empirical gain model derived by Campbell et al. (2021). MRO's attitude has varied across the mission to accommodate specific observations. It has been shown that the S/N for SHARAD is strongly dependent on the roll angle of the spacecraft, but an accurate correction model is yet to be obtained (Campbell et al. 2021). To limit any nonquantifiable S/N variations due to MRO's attitude, we excluded data with an absolute roll angle greater than 0°.5 around the axis corresponding to the velocity vector of MRO's orbit around Mars.

The surface echo, or surface return, is the signal reflected by the atmosphere–surface interface. It is the earliest arrival recorded at the antenna affected only by ionospheric attenuation (Mouginot et al. 2008; Campbell et al. 2014). We mitigate this effect as described in Section 2.3.

Scattering volume. The strength of the surface echo is the summation of all the elementary electric fields reflected and scattered by the dielectric contrasts within a volume bounded horizontally by the circular radar footprint and vertically to the near-surface depth (<15 m, depending on the dielectric properties of the ground; Grima et al. 2012; Figure 1(a)). The size of the radar footprint for a perfectly smooth surface is given by the Fresnel zone, an area from which the scattered signals interfere constructively (Haynes 2020). However, since the incoherent scattering is fully integrated into our analyses, we shall consider the larger pulse-limited footprint that accounts for the contribution of peripheral scatterers. The resulting effective footprint in the absence of along-track focusing is about 5 km in diameter for SHARAD (Seu et al. 2007). Therefore, the echo strength holds crucial information about both the surface roughness and the near-surface composition and structure (Figure 1(b)).

Zoom In Zoom Out Reset image size

Components. It can be demonstrated analytically that the total power of the surface return is composed of a coherent component (or reflectance, P_c ) and an incoherent component (or scattering, P_n ), so that the total power received is P_t = P_c + P_n (e.g., Ishimaru 1978; Tsang 2001; Ulaby & Long 2014). By opposition to the coherent component, P_n pertains to a diffused random field that is related to the scatterers randomly distributed within the scattering volume (volumic inclusions) and at its interface (surface roughness; Figure 1(b)). In the absence of volume scatterers, P_c and P_n are equally scaled by the surface reflection coefficient that is determined by the apparent permittivity (ε) of the ground (see Section 4.2 for basic analytic relationships).

Apparent permittivity. ε is a first indicator of the ground composition since geologic materials (e.g., ices and sedimentary, basaltic, metamorphic rocks) lie in various, although often overlapping, ranges of permittivity (e.g., Campbell & Ulrichs 1969; Telford 1990). ε is also dependent on the volume fraction of impurities (including void spaces) in the medium (Sihvola 1999; Martinez & Byrnes 2001; Grima et al. 2009).

Layering. The apparent permittivity sensed by the radar is further modulated by the deterministic and nearly horizontal structure of the ground that generates specular reflections that may interfere with each other. The term deterministic has to be understood as a structure that is not randomly varying across the footprint such that the outgoing wave front of the electromagnetic field is spatially undisturbed and steady. Depending on the geologic context, such deterministic structures could come from the overburden floor, the ice table, or aquifers present in the near surface (e.g., Mouginot et al. 2009; Grima et al. 2016; Scanlan et al. 2022). The overall effect of near-surface layering is extremely sensitive to the ratio between the signal wavelength and the layers' thicknesses (Mouginot et al. 2009; Lalich et al. 2019).

Surface roughness. Surface roughness at a horizontal scale that equals the wavelength is usually the primary source of randomly distributed scatterers. The tilted elements of the surface, usually associated to the random surface slopes, are mainly responsible for diffusing the backscattered field into different, not-specifically specular directions. Distinctly, the vertical distribution of surface heights is the main roughness attribute that is ubiquitous to both the deterministic and random surface responses through the P_c and P_n components, respectively. Its physical effect on ¶_c is very similar to a layer with a thickness that approximates the root mean square of the heights.

Volumic inclusions. Incoherent scattering could also come from near-surface inclusions such as ice lenses, buried blocks, or strongly brecciated media with typical dimensions representing a fraction of the radar wavelength (Rutishauser et al. 2016). As a first approximation, volume scattering is isotropic (Lambertian; Nayar et al. 1991), therefore with a minimal specular contribution. Such volume scattering is ignored in Section 4 related to surface roughness approximation, but might contribute to enhance the incoherent component in very specific settings.

Introduction. One can attempt to describe the randomness of the surface echo amplitude from its envelope distribution with a probability density function (PDF) assuming specific assumptions for the statistical properties of the surface. Grima et al. (2012) demonstrated that the statistics of the SHARAD surface echoes for dominantly coherent signals (i.e., smooth terrains) are well reproduced by a Rice envelope (Rice 1945). Conversely, when the coherent signal shrinks over rougher terrains, the K-distribution (Jakeman & Pusey 1976) better reproduces the signal statistics over the most common Rayleigh distribution used for homogeneous scattering medium (Grima et al. 2012). This observation is expected since a K-noise is a generalization of Rayleigh scattering that allows the scatterers to be clustered, an assumption certainly more robust to model scattering from planetary bare surfaces. The Rice and the K-distribution have been bridged together by the homodyned K-distribution (HK; Jakeman 1980; Jakeman & Tough 1987; Dutt & Greenleaf 1994) that we therefore propose to use in this study as a more suitable model for systematic application across the Martian surface. The flexibility of the HK has been leveraged in many diverse applications ranging from ultrasound imagery, synthetic aperture radar imagery, and ice-penetrating radar reflectometry (e.g., Tison et al. 2004; Grima 2014; Tresansky 2020). Grima (2014) provides a brief overview of the HK relationships with basic PDFs in the context of nadir-pointing radar surface reflections. Comprehensive and extensive reviews for other applications can also be found in Destrempes & Cloutier (2010) and Tresansky (2020).

Definition. The HK is analytically built from the summation of N phasors (reproducing the effects of N interfering scatterers) in a 2D random walk following a negative-binomial distribution, and into which is added one deterministic (or homodyning) component (Dutt & Greenleaf 1994). Physically this allows the N scatterers to be clustered within the radar footprint (i.e., not necessarily independent and identically distributed). The law of large number of scatterers, typical for Rayleigh scattering, does not have to be fulfilled. The HK also appears to be the only model for which the parameters keep their physical meaning in the limiting case of vanishing homogeneity, while remaining valid for a high density of uniformly distributed scatterers (Dutt & Greenleaf 1994; Destrempes & Cloutier 2010). The HK distribution is given by

$$\begin{eqnarray}&&{P}_{\mathrm{HK}}(A| a,{s}^{2},\mu )=A{\int }_{0}^{\infty }{{uJ}}_{0}({ua}){J}_{0}({uA}){\left(1+\displaystyle \frac{{u}^{2}{s}^{2}}{2}\right)}^{-\mu }{du},\end{eqnarray} \tag{ 1 }$$

where A is the signal amplitude, J₀(•) is zeroth-order Bessel function of the first kind, u is a nonspecific variable of integration. The coherent and incoherent powers are given by P_c = a² and P_n = 2s² μ, respectively. Figure 2 illustrates the typical shapes of HK distributions for various sets of of parameters.

Zoom In Zoom Out Reset image size

Physical meaning. a is the constant amplitude added to the random walk and characterizes the coherence of the signal. It controls the location of the HK envelope mode. Physically, it results from the coherent summation of superposed fields reflected by flat interfaces nearly perpendicular to the propagation path: usually the surface boundary in the first place, but also any sharp dielectric gradient in the near surface as well as any stair-like shape of those interfaces within the footprint. 2s² is the scale parameter of the HK. It somewhat represents the scattering strength of an individual scatterer. The total scattering strength is obtained by multiplying 2s² by μ.

μ is a parameter of the negative-binomial distribution used to model the 2D random walk associated with the scattering (Dutt & Greenleaf 1994). It is morphologically related to the kurtosis of the HK envelope and is called the scatterer number density parameter or clustering parameter. This double designation underlines its intricate meaning for the disorganized structure of the target within a resolution cell. Studies suggest that a resolution cell could be translated to a λ³ volume, where λ is the sensing wavelength. However, most of the observations reported in this paragraph are derived from the context of ultrasound biologic tissue characterization where the scattering media is usually a volume, although a 2D surface-like space is often considered for model simplicity. A comprehensive effort to understand the quantitative meaning of μ in the context of radar sensing of bare planetary surfaces still remains to be undertaken, but advancements made in other fields provide a valuable framework. Simulations and phantom experiments suggest that μ increases monotonically with the number of scatterers per resolution cell so that it provides a relative indication of the density of scatterers (e.g., Abraham & Lyons 2002; Cristea et al. 2020). The relationship appears almost linear for independent and identically distributed random samples (Dutt & Greenleaf 1994; Cristea et al. 2020). It is noteworthy that similar results are obtained in the field of seafloor sonar sensing (Abraham & Lyons 2002). Geologic processes do not necessarily produce independent and identically distributed random scatterers. In this case, the HK theory also predicts μ to vary with the clustering of the scatterers, i.e., the degree to which the scatterers tend to be spatially gathered in patches (Dutt & Greenleaf 1994). This is validated by experiments showing that μ increases with the decreasing degree of clustering (Hu et al. 2017). This trend appears to break when the media loses any signature of entropy either by reaching a certain saturation threshold for the scatterer population or when the scatterers become regularly distributed (i.e., without randomness) [ibid.].

The literature often makes use of a ratio that combines a, s, and μ as a relative measure of coherence within the signal. For instance, Destrempes et al. (2013) uses an expression of the coherent-to-diffuse signal ratio $k=a/(s\sqrt{\mu })$ that appears to increase linearly with the periodicity of scatterers, while Tresansky (2020) proposes the "coherent fraction" as the relative power in the total signal P_c /(P_c + P_n ) = a²/(a² + 2s² μ) that carries a more intuitive meaning. One common motive is to produce a scale-independent parameter that is insensitive to signal attenuation as well as Fresnel reflection and transmission losses. Here we use the "coherent content" defined by P_c /Pn = a²/(2s² μ) as proposed by Grima et al. (2012,2014b). In a way similar to other ratios, the coherent content is theoretically only dependent on roughness and near-surface disorder although its relationship with heterogeneity descriptors is not monotonic (Section 4.2). The coherent-to-incoherent balance illustrated by the coherent content is also a first-order indicator to drive the choice of a relevant backscattering model for a further physical properties inversion as it indicates what type of power is dominant in the signal. Morphologically, a negative coherent content or fraction shifts the HK envelope morphology toward a positive skewness.

Standard estimation techniques for the derivation of the HK parameters include a variety of methods of moment (MoM) algorithms based on the combination of moments for the assumed HK-distributed empirical data (e.g., Dutt & Greenleaf 1995; Hruska & Oelze 2009; Destrempes et al. 2013; Haynes 2019), as well as maximum likelihood estimators (e.g., Hu et al. 2017). Here, we use a simple curve-fitting estimation algorithm ⁸ based on a nonlinear least squares fitting of the surface echo amplitude distribution (Markwardt 2008). Initial conditions a₀ and s₀ are obtained for each fit from the mean and standard deviation (STD) of the empirical data, respectively, while μ₀ is initially set to unity. More relevant initial conditions could also be fed from a MoM estimator at the cost of computation time (Tresansky 2020). Each iteration of the curve-fitting process is constrained by the condition that the sum of the derived P_c and P_n values should equal the total signal power (defined by the mean of the set of echoes considered) in order to match the theoretical relationship between the signal components (Section 4.2). The bin width for an amplitude histogram is determined with the rule from Freedman & Diaconis (1981) that estimates the scale of the distribution from its interquartile range. The apparent robustness of curve-fitting estimators against noise makes them more appropriate when the size of the data set is low, typically on the order of 100–1000 s of points as is the case for our data set (Section 3.5; Tresansky 2020). The error estimation of our algorithm is assessed in the following section.

We investigate the error with which our algorithm determines both the coherent and incoherent powers of the surface return as well as the surface clustering parameter. Our core methodology is to draw a synthetic set of 1000 pseudo-random amplitudes following the inverse transform sampling (ITS) technique (Devroye 1986) applied to an HK noise with predefined parameters. We add an arbitrary Gaussian noise of ∼1 dB to simulate the effect of nongeophysical noise (e.g., instrumental) introduced in the measured signal. Haynes (2019) proposed to account for the instrumental noise by adding a random variable (i.e., an additional unknown) to the analytic expression of the PDF used. While we do not do that, our approach enables to assess how our estimator behaves when an unknown and limited non-HK instrumental noise is introduced. We then apply our algorithm to this set of 1000 synthetic amplitudes to derive the predicted HK parameters. This operation is repeated N = 1000 times for a given combination of predefined (true) HK parameters within a range encompassing 99% of the empirical measurements in our study area.

Results are presented in Figure 3 in terms of the relative bias $(E[\hat{y}]-y)/y$ and the normalized STD $\sqrt{\mathrm{Var}[\hat{y}]}/y$ of a predicted value $\hat{y}$ with respect to a true value y. The relative bias and normalized STD can both be interpreted in terms of scaled measurement accuracy and precision, respectively (Walther & Moore 2005). Figure 4 provides an alternative visualization through the true-versus-predicted plots for the coherent content and μ.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Overall, the measured coherent content (P_c /P_n ) emerges as a metric to assess the reliability of the predicted values, as noted by former studies (e.g., Hruska & Oelze 2009; Destrempes et al. 2013). In accordance with the uncertainty principle, when a power component weakens relatively to the other one, the uncertainty on its predicted value worsens. Still, confidence ranges exist: When the coherent content is measured between ±5 dB (resp., <10 dB), the relative bias on the predicted ${\hat{P}}_{c}$ and ${\hat{P}}_{n}$ is <1 dB (resp., <2 dB). The normalized STD is <0.5 for coherent contents down to −5 dB. The prediction of μ is highly sensitive to the coherent content in a manner that might preclude using it for further qualitative interpretations. However, comparative analysis can be attempted for surfaces with similar coherent contents as their true-versus-predicted plots are monotonic (Figure 4, right).

Application. The region of interest is spatially segmented into circles of diameter d = 017 (∼10 km at the equator) and spaced every d/2 in both longitude and latitude. Our RSR algorithm is then run over the set of surface echoes contained within each circle. A single set gathers about 1000 echoes in average. Results for the RSR-derived P_c , P_n , P_c /P_n coherent content, and μ are presented in Figure 5. The goodness of fit between the amplitude distributions and the corresponding HK fits is also assessed through their respective correlation coefficients in the bottom right panel of Figure 5 with most of it being >80% with a median at 98%. For illustration, typical amplitude distributions and their associated HK fit in this region can be found in Putzig et al. (2017, their Figure 5(b)).

Zoom In Zoom Out Reset image size

Observations. We observe that the relative dynamic range of the total signal spans about 30 dB, similar to what was observed by Grima et al. (2012) across a significant portion of the Martian surface, confirming that our selected case-study region is mostly representative of the power dynamic that could be encountered throughout the planet. Noticeably the full extent of this 30 dB dynamic range is mainly driven by P_n whereas the coherent signal P_c is due to geologic variations (i.e., not considering the background noise signal as determined below) and only spans ∼15 dB. The total radiometric range could be artificially reduced by the RSR-induced bias highlighted in our sensitivity analysis (Section 3.4) that is especially important for low coherent contents associated with high μ values. However, the facts that the derived μ values are mostly lower than 2 and that Grima et al. (2012) measured a similar dynamic range from the raw surface return (i.e., without resorting to the RSR methodology) suggest that this range underestimation is a minor effect.

Another important observation comes from the bulk of the measured coherent contents in decibels being positive, implying that the signal return at the antenna is often dominantly coherent. This is especially of interest for driving SHARAD data users toward an appropriate backscattering model for a surface properties inversion in a specific region. In particular, backscattering models for a dominantly coherent signal have a more tractable form and a reduced number of surface property unknowns than for a more balanced coherent–incoherent power location (e.g., Ulaby et al. 1981; Ogilvy 1991). Interestingly, rough terrains such as the MFF and some Noachian units exhibit a strong coherent content despite intuitively being expected to scatter the signal. In some cases, this coherency could be explained by the incoherent power fall-off being greater than for the coherent one (Haynes et al. 2018), so that the distance to the target naturally filters out the relative strength of the incoherent signal. In some other cases, a positive coherent content is simply the manifestation of the surface signal being so weak that our surface picker samples the random noise in a biased manner. The regions where these effects are dominant can be identified and discriminated as described in the following section.

Background noise. In locations where the surface signal is too dimmed to be reliably detected, the derived P_c and P_n components might be misinterpreted as proxies for the surface properties while they are in fact signatures of the background noise. To predict this effect and avoid geophysical interpretation pitfalls, we characterize the coherent signature of the background noise as seen by our surface picker. The purpose is not to characterize the noise itself, but rather to give insights on how the noise is rendered in the P_c − P_n space by the selective signal determination of our picker. For each radargram crossing, our region of interest, we run our surface picker into the first 100 earliest fast-time samples where the surface is usually absent. This is repeated 1000 times along track over equally spaced range lines. The obtained set of 1000 amplitudes is used to determine a single pair of P_c and P_n values per radargram through the RSR algorithm. The results reported in Figure 6 are within a range of relative power defined by P_c < 60 dB and P_n < 50 dB, and restricted into a narrow P_c /P_n ratio of 12 ± 0.5 dB.

Zoom In Zoom Out Reset image size

Thanks to this range of values, one can identify the regions that correspond to background noise. The bottom left panel of Figure 5 represents the distribution of our samples in the 2D P_c − P_n space. The bottom left end of this distribution is split in two groups, both with P_n < 50 dB, but one with P_c > 60 dB while the other is coherently weaker than 60 dB. This latter group of points corresponds to the range of power for the background noise as outlined by the overlapping red contours. Interestingly, those two groups of points can be both geolocated within the MFF, indicating that the S/N of the surface echo there oscillates around the limit of detection for our surface picker. A closer look at the map of the incoherent power (Figure 5, upper right) shows that similar powers within the MFF tend to be arranged along track. It suggests that the S/N oscillations triggering the surface echo detection are driven by the conditions of acquisition (e.g., SA-HGA configuration, ionosphere attenuation) specific to each track.

The behavior of the distribution in the 2D P_c − P_n space (bottom left panel of Figure 5) exhibits a noteworthy comma-like shape that we interpret in the context of signal growth with increasing surface roughness. This interpretation is detailed in the paragraph below through three types of regimes and illustrated by Figure 7. It will drive our logic for roughness approximation in the Section 4.2.

Zoom In Zoom Out Reset image size

When the roughness increases over that of a flat surface (regime (1)) to a very rough surface with quasi-isotropic scattering properties (regime (3)), the coherent power diminishes until it eventually becomes undetectable. The coherent power lost thereby is progressively transferred into the pool of incoherent power, so that the incoherent power integrated over the half-space above the surface grows continuously with roughness. However, the incoherent power that is intercepted by the spacecraft at normal incidence does not exhibit a monotonic behavior. It first grows with roughness around the specular direction. But then, when the tilted surface elements reach a certain threshold depicted by regime 2, the growth of the incoherent power tends to be concentrated in off-normal directions, and counterintuitively, the incoherent power intercepted by the antenna around the normal then diminishes (Ulaby et al. 1981), going from regime (2) to regime (3), although the incoherent power integrated in the upper half-space increases. The coherent content provides another signature of this entire dynamic. It drops rapidly between regimes (1) and (2) since P_c and P_n evolve in an opposite manner. Then, the coherent content nearly stabilizes or even slightly increases, as both P_c and P_n shrink.

The analytical relationships that bound the coherent and incoherent energies to the physical properties of the surface can be given by backscattering models. In general, at a given altitude z from the surface, at normal incidence, without volume scattering and for a stationary and ergodic surface, those components can be written in the forms

$$\begin{eqnarray}&&{P}_{c}=\alpha {r}^{2}\chi {\left({\sigma }_{h}\right)}^{2},\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}&&{P}_{n}=\alpha {r}^{2}L(z)\zeta {\left({\sigma }_{h},{l}_{c}\right)}^{2},\end{eqnarray} \tag{ 3 }$$

where α is a calibration constant that can be adjusted so that P_c = 1 (or 0 dB loss) when the surface is a perfectly flat and conductive reflector (Fresnel coefficient r = 1). r² is the Fresnel reflection coefficient that is a function of the surface permittivity such as follows:

$$\begin{eqnarray}&&{r}^{2}={\left|\displaystyle \frac{1-\sqrt{\varepsilon }}{1+\sqrt{\varepsilon }}\right|}^{2},\hspace{.2in}\mathrm{for}\hspace{.1in}\varepsilon \geqslant 1\hspace{.2in}\mathrm{and}\hspace{.2in}0\lt {r}^{2}\lt 1.\end{eqnarray} \tag{ 4 }$$

L(z) is the relative geometric propagation loss with respect to a specular reflection. For a flat surface, L(z) = 1/(π z)² (Grima et al. 2012). The surface roughness properties are defined by the rms height (σ_h ) and correlation length (l_c ). For a more intuitive interpretation, we will refer to the effective slope s_e = σ_h /l_c (Campbell & Garvin 1993; Shepard et al. 2001). It is equivalent to the rms slope for a surface profile with an exponential autocorrelation (Ogilvy 1991). $\chi {({\sigma }_{h})}^{2}$ and $\zeta {({\sigma }_{h},{l}_{c})}^{2}$ are functions describing losses due to roughness. These are both dependent on the backscattering model considered. From Equations (2) and (3), it follows that the coherent content is independent on the surface permittivity and is only a function of the two roughness parameters when the altitude is known:

$$\begin{eqnarray}&&\displaystyle \frac{{P}_{c}}{{P}_{n}}=\displaystyle \frac{\chi {\left({\sigma }_{h}\right)}^{2}}{L(z)\zeta ({\sigma }_{h},{l}_{c})}.\end{eqnarray} \tag{ 5 }$$

When the surface backscatter tends toward regime (1) in Figure 7, the coherent power is weakly affected by roughness and is mainly dependent on the surface permittivity through r². Our study area encompasses the transition between regime (1) to regime (2) through a distinguishable nearly horizontal distribution of measurements on Figure 7. The surface backscatter tends toward regime (1) for ${P}_{{c}_{0}}\approx 76\,{dB}$. We use this value as a reference for the surface echo strength reflected from a flat surface (σ_h = 0) with a given, but unknown, reflection coefficient ${r}_{0}^{2}$. A surface with the same r² that is affected by roughness will lose coherent power with respect to ${P}_{{c}_{0}}$ by a factor $\chi {({\sigma }_{h})}^{2}$. This allows to associate each P_c measurements to a σ_h value. To this end, we use the analytical relationship $\chi {({\sigma }_{h})}^{2}\,=\,{e}^{-{(2k{\sigma }_{h})}^{2}}$, where k is the wavenumber, as it tends to be an ubiquitous relationship across standard backscattering models as a first-order approximation (e.g., Ulaby et al. 1981; Ogilvy 1991).

Each measurement is now associated to a coherent content P_c /P_n and an rms height σ_h . Then, both can be injected into Equation (5) to derive a correlation length, where ζ(σ_h , l_c ) is evaluated from the simplified integral equation method (IEM) that is valid for a stationary surface with normally distributed heights, with rms slopes <0.3 (∼17°) and with k σ_h < 2 (i.e., σ_h < 4.8 m at SHARAD wavelength) (Fung & Chen 2010).

The approximation for the surface roughness parameters obtained from this methodology is shown through the red grid on Figure 8. We note that the distribution of measurements lie within the application limits of the simplified IEM (Section 4.2).

Zoom In Zoom Out Reset image size

This inversion method allows for the quantitative approximation of the surface roughness parameters at the radar wavelength scale, and without the support of a third-party data set. However, it requires the fundamental assumption that the reflection coefficient ${r}_{0}^{2}$ (or similarly permittivity, ε₀) of theflat-surface reference is the same for all the terrains investigated. When this is not the case, the bias (or accuracy) induced on the approximated roughness can be estimated, provided that reasonable assumptions can be made regarding the permittivity difference between the two terrains. Let ${\rm{\Delta }}{r}^{2}={r}_{0}^{2}-{r}_{i}^{2}$ be the relative difference between the Fresnel coefficient of the flat-surface reference and another considered terrain r² _i . Remember that both P_c and P_n are equally proportional to r² (Equations (2) and (3)). It follows that the true roughness of the considered terrain can be obtained by (1) locating the position of its measured $({P}_{{c}_{i}},{P}_{{n}_{i}})$ pair on Figure 8; (2) sliding diagonally from their by an amount Δr² in both the x- and y- directions (i.e., parallel to the coherent content isolines); (3) the roughness values given at this new location on the plot are the true ones and can be compared with the roughness obtained at the initial $({P}_{{c}_{i}},{P}_{{n}_{i}})$ location to estimate the bias.

For example, let us consider a terrain with measured surface components of $({P}_{{c}_{i}},{P}_{{n}_{i}})=(65,60)$ dB. The approximated roughness from Figure 8 is σ_h = 1.8 m and s_e = 12. Let us now consider ε₀ = 9 and ε_i = 4. These permittivities correspond to typical end values for materials of igneous origins (Ulaby et al. 1981), which we assume to be covering most of the studied regions as a first approximation. Δr² is then equal to +3.5 dB, and the derived roughness parameters corrected from Δr² are σ_h ≈ 1.4 m and s_e ≈ 10. Therefore, the accuracy of the derived roughness grows with Δr², but it remains limited, especially when the surface is already rough. For nonvolcanic materials like nearly clean H₂ O ice, ε_i = 3.1 (Grima et al. 2009), ${{\rm{\Delta }}}_{r}^{2}$ with a reference volcanic terrain of ε₀ = 9 is 5.2 dB. In that case, a new flat-surface reference over an icy terrain would certainly provide better accuracy. When the surface is very smooth, the roughness bias due to Δr² grows faster and might require other approaches to better assess the surface properties, such as the inversion technique combining the small angle approximation with the assumption of large correlation length (Grima et al. 2012, 2014a).

The bias due to the relative ignorance of the permittivity difference, between the reference terrain and the studied terrain, has to be combined with the uncertainties on P_c and P_n if the coherent content is outside of the confident ranges defined in Section 3.4.

The methodology described in Section 4.3 provides a pair of roughness parameters (σ_h , s_e ) for each measured pair of coherent and incoherent surface echo powers (P_c , P_n ). The results are mapped on the top row of Figure 9. The bottom row provides comparable parameters independently obtained from the MOLA, but with different inversion techniques and at larger horizontal scales than SHARAD (Kreslavsky & Head 2000; Neumann 2003). A detailed study of the relationship between the SHARAD-derived roughness with the depositional and erosional geologic processes particular to this region will be addressed in a follow-on study. Our aim here is to discuss the validity of the so-derived roughness parameters by comparison with MOLA products. This work can be seen as an extension to the discussion initiated by Grima et al. (2012, Section 3) on the comparison between SHARAD total reflectivity and MOLA roughness.

Zoom In Zoom Out Reset image size

Surface slopes. The MOLA median differential-slope product at 600 m scale (Figure 9, bottom left) is obtained from the median absolute values of the slopes measured at a 600 m horizontal baseline length from altimetry, but detrended from the regional slope measured at twice this scale (Kreslavsky & Head 2000). The removal of the regional topographic trend is responsible for the low, generally subdegree, median differential slope, making it hard to quantitatively compare with SHARAD roughness. However, the MOLA differential-slope product provides a good indicator outlining with fine sensitivity the spatial variations of intrinsic surface roughness that are really due to features at the baseline scale considered. The SHARAD-derived s_e at a 15 m baseline exhibits similar trends. The large outliers of the MFF (known to be rough at decameters scale; Mandt et al. 2008) and cratered NR appear both highly rough on SHARAD and MOLA slopes. More spatial variations can be distinguished within the MFF on the MOLA data set, while the SHARAD signal strength saturates at about s_e = 5°. This corresponds to the limit beyond which the signal strength reaches the background noise (Figure 5, bottom right), then becoming barely detectable. This behavior demonstrates that the upper limit of SHARAD sensitivity to roughness is controlled by its S/N capability to overcome extinction from strong surface scattering. Outside of those rough areas, many spatial nuances suggesting lava flow margins and/or structures in the smooth EP are mostly recovered in SHARAD-derived slopes similar to MOLA's, outlining small-scale variations in EP between 150°E and 155°E, for instance. The ES appear in both maps as terrains of intermediate roughness getting locally rougher within crater rims and ejecta.

Surface heights. The MOLA rms height (Figure 9, bottom right) is obtained from the relative spreading of the altimetric laser return after reflection onto the surface Neumann (2003). The result is a measure of vertical roughness at a horizontal baseline similar to the laser footprint diameter of ∼75 m, i.e., the smallest baseline at which roughness can be retrieved with global coverage so far. Theoretically, the linear additivity of the various sources of errors predicts that rms heights as low as 1 m can be resolved. Yet, the nuances of lava flows within EP are not recovered, while the transition from EP to ES is barely noticeable. The only terrains where inner variations are distinguishable are within the rougher MFF. By contrast, the SHARAD-derived rms heights at 15 m baseline still exhibit the subtle contrasts observed with the MOLA median differential slopes at 600 m baseline with the various geologic units considered. EP especially shows a large spectra of roughness from intermediate values of σ_h = 1.8–2 m gathered in a north–south band centered at 157° longitude that transitions abruptly eastward to a very smooth and narrow terrain with σ_h < 1 m, also seen on the MOLA median differential slopes. This latter terrain mostly contributes to the flat reference surface used in Section 4 to produce our roughness approximation grid. Its coherent content up to 15 dB is higher than that for the south polar layered deposits measured at 12 dB by Grima et al. (2012), and thus is one of the smoothest terrains observed on Mars.

SHARAD-MOLA complementarity. The rms heights derived from SHARAD and MOLA are quantitatively compared in Figure 10 (left). Both data sets measured at different baselines are not expected to exhibit self-similarity (i.e., lying around the identity line) since bare planetary terrains are mostly observed to be self-affine (Shepard et al. 2001). However, Figure 10 displays specific behaviors worth mentioning. Two group of points can be distinguished that illustrate saturation effects related to both data products and their respective sensitivity range. A horizontal group of points aligned around the y coordinate 2.25 m outlines the detection limit beyond which the SHARAD surface return is lost into the background noise (Section 3.5.0.0). A vertical group of points mainly slightly below x-coordinate 1 m corresponds to measurements within ES and EP. These illustrate the lower-resolution limit of MOLA pulse-width measurements.

Zoom In Zoom Out Reset image size

Those detection thresholds define ranges where both data sets are sensitive to roughness variation. Again, this range for SHARAD is <2.25 m, defining where the surface signal is not pure noise. It is >1 m for MOLA, a range where pulse-width variations are distinguishable from the pulse expected from a perfectly flat surface (Neumann 2003). Mapping the locations where those sensitivity ranges are achieved (Figure 10, right) provides highlights on how complementary the MOLA pulse-width and SHARAD roughness are when it comes to quantify relative rms heights at decametric scales. The roughest terrains (MFF and NR) are mainly within the detection limits of MOLA only (blue on Figure 10, right). By contrast, SHARAD provides the only measurement able to distinguish roughness variations throughout the smoothest terrains (red) like Elysium Planitia (Kreslavsky & Head 2000). Various transitional terrains between those two roughness behaviors can be characterized by both data sets (yellow). Thus, SHARAD complements MOLA at decametric scales since it can sense roughness variations over smooth terrains where the MOLA pulse-width technique is below its detection threshold. This can be informative for a finer quantitative comparison of smooth terrains to, for example, support the selection process of landing site candidates (e.g., Golombek et al. 2017) while small-scale measurements from MOLA pulse-width do not exhibit any contrast between sites. Note that two landing sites are covered by our study. SHARAD roughness of the most nearby measurements are about σ_h = 1.5 m and s_e = 13 for the Curiosity landing site, and σ_h = 1.9 m and s_e = 09 for the Insight landing site. MOLA provides σ_h 1.2 m and σ_h < 1 m for Curiosity and Insight, respectively.

SHARAD-MOLA discrepancy. It is remarkable and counterintuitive that the rms heights are usually higher for SHARAD than for MOLA, while the SHARAD baseline length of 15 m is lower. Indeed, fractal laws predict that height variations of self-affine surfaces decrease with the horizontal scale at which they are measured (Shepard et al. 2001). We propose that this discrepancy is explained by the 2 orders of magnitude difference between the footprint size of SHARAD and that of MOLA. Indeed, although the SHARAD signal interacts with wavelength-sized (15 m) objects, the electromagnetic field received at the antenna is the coherent summation of all the fields scattered within the larger 5 km wide footprint. Hence, any large-scale trend of the surface, like kilometer-wide undulations, would increase the height range of the scatterers that contribute to the total phase of the wave front reaching the antenna, then increasing the approximated rms height. The effect of large-scale topography is much less pronounced for the MOLA pulse-length since the footprint contributing to the laser signal is only 70 m wide. To some extent, this effect is very similar to the MOLA median differential slopes giving much smaller values than those from altimetric roughness measurements not detrended from the regional slopes from Aharonson et al. (2001), Kreslavsky & Head (2000).

This hypothesis suggests that any quantitative multiscale comparison of radar-derived roughness should be wary of the footprint extent at which measurements are obtained in addition to the sounding wavelength. In addition, it also suggests that very small rms heights approximated from SHARAD (e.g., the submeter reference terrains within EP) are not only smooth but also necessarily extremely level across the 5 km footprint. These observations advocate for the use of fractal (self-affine) backscattering models whose surface roughness parameterization accounts, in essence, for variations at different scales (Franceschetti et al. 1999; Shepard & Campbell 1999; Campbell & Shepard 2003.