Texture and composition of Titan's equatorial sand seas inferred from Cassini SAR data: Implications for aeolian transport and dune morphodynamics

TL;DR

This study analyzes Cassini SAR data to determine the texture, composition, and morphodynamics of Titan's equatorial dunes, revealing differences between dunes and interdunes and insights into active sediment transport and climatic history.

Contribution

It provides the first quantitative analysis of Titan's dune surface properties across a wide range of incidence angles using advanced modeling of SAR data.

Findings

Dunes are more microwave absorbent than interdunes.

Interdunes have higher dielectric constants and multi-scale roughness.

Evidence suggests active sediment transport and organic sand composition.

Abstract

The texture, composition, and morphology of dunes observed in the equatorial regions of Titan may reflect present and/or past climatic conditions. Determining the physio-chemical properties and the morphodynamics of Titan's dunes is therefore essential to understanding of the climatic and geological history of the largest moon of Saturn. We quantitatively derived average surface properties of dune and interdune areas (texture, composition) from modeling of the microwave backscattered signal and Monte-Carlo inversion of the despeckled Cassini/SAR data over Titan's three largest sand seas: Belet, Shangri-La and Fensal. We present the first analysis of the backscatter functions extracted from despeckled SAR images that cover such a large range in incidence angles, including data from the beginning of the Cassini mission up to its Grand Finale. We show that dunes and interdunes have significantly different physical properties. Dunes are found to be more microwave absorbent than interdunes. Additionally, potential secondary bedforms, such as ripples and avalanches, may have been detected, providing potential evidence for currently active dunes and sediment transport. Our modelling shows that the interdunes have multi-scale roughnesses with higher dielectric constants than the dunes which have a low dielectric constant consistent with organic sand. The radar brightness of the interdunes can be explained by the presence of a shallow layer of significantly larger organic grains, possibly non-mobilized by the winds. {Together, our} findings suggest that Titan's sand seas evolve under the current multi-directional wind regimes with dunes that elongate with their crests aligned in the residual drift direction.

Equations

$$\xi =\sqrt{\frac{{\sum }_{i=1}^n(z_i-\bar{z}{)}^2}{n-1}},$$

$$\rho_{N}(\delta)=\frac{1}{\xi^{2}}\big{<}z(x_{i}+\delta)z(x_{i})\big{>},$$

$$\zeta ={\rho }_N^{-1}(\exp (-1)).$$

$$\xi >\frac{\lambda }{8\ast \cos (\theta )},$$

$${\sigma }_{Shh}^0=\Gamma (\theta =0)\frac{\exp (-{\tan }^2(\theta )/2m^2)}{2m^2{\cos }^4(\theta )},$$

$${\sigma }_{Vhh}^0=\frac{3}{4}a{T}_{ij}^2\cos (\theta )\left(1-e^{-2\tau {\cos }^{-1}({\theta }_t)}\right),$$

$${\sigma }_{IEM}^0=\frac{k\exp (-2k^2{\xi }^2{\cos }^2(\theta ))}{4}\sum _{n=1}^{\infty }|I_{hh}^n{|}^2\frac{W^n(-2k\sin (\theta ),0)}{n!},$$

$$I_{hh}^n=$$

$${\sigma }_{POM}^0=$$

$${\sigma }_{GOM}^0=\Gamma (\theta =0)\frac{\exp (-{\tan }^2(\theta )/2m^2)}{2m^2{\cos }^4(\theta )},$$

$$\sigma^{0}_{Shh}=\left\{\begin{array}[]{ll}\sigma^{0}_{IEM}&\textrm{if }{}k\xi<3\textrm{ and }{}k^{2}\xi\zeta<1.5\sqrt{\epsilon},\\ \sigma^{0}_{POM}&\textrm{if }{}k\zeta>6\textrm{ and }{}k^{-1}<\xi<0.06k\zeta^{2},\\ \sigma^{0}_{GOM}&\textrm{if }{}(2k\xi\cos\theta)^{2}>10\textrm{ and }k\zeta>6.\end{array}\right.$$

Full text

\journalname

JGR-Planets

Université de Paris, Institut de physique du globe de Paris, CNRS, Paris, France LMD,UPMC, CNRS, France LATMOS/IPSL, UVSQ Université Paris-Saclay, CNRS, Guyancourt, France Department of Physics, University of Idaho, Moscow, ID, USA Observatoire de Bordeaux, Université de Bordeaux, CNRS, Bordeaux, France

\correspondingauthor

A. Lucaslucas@ipgp.fr

{keypoints}

Physical properties over Titan sand seas are quantitatively derived from Bayesian inference of the Cassini SAR amplitude data

Interdunes terrains are made of fine and coarse parsed grains over a polluted water-ice bedrock

Analogy with terrestrial sand seas in Niger and China suggests that dunes on Titan reflect current climatic conditions

1 Introduction

After 13 years of exploring and revealing the Saturn’s system, in particular its largest moon Titan, the Cassini spacecraft plunged \replacedon Saturninto Saturn’s upper atmosphere on Friday 15th of September 2017. Among the most prominent of Titan’s landforms discovered by the Cassini-Huygens mission are the large sand seas that cover 15-20$\%$ of the global surface of the moon, located mostly in the equatorial regions (Lorenz et al., 2006; Radebaugh et al., 2008; Radebaugh, 2013; Le Gall et al., 2011; Rodriguez et al., 2014; Aharonson et al., 2014). As inferred from the spectral observations of the VIMS (Visual and Infrared Mapping Spectrometer) instrument (Soderblom et al., 2007; Barnes et al., 2008; Rodriguez et al., 2014) and \deletedsupported by the surface properties derived from the radiometry mode of the RADAR instrument (Le Gall et al., 2011), these vast sand seas constitute a major reservoir of organics at the surface. Consequently, they play a prime role in the geological and climatic history of Titan (Aharonson et al., 2014).

Thanks to the high spatial resolution of the 5 beam Synthetic Aperture Radar (SAR) imager mode of the RADAR instrument– operating at 2.17 cm wavelength (Ku-band) in HH polarization (Elachi et al., 2004)–individual linear dunes have been identified within those sand seas and mapped at the global scale (Lorenz et al., 2006; Radebaugh et al., 2008; Le Gall et al., 2011; Rodriguez et al., 2014). These sand seas present a dark signature in visible (ISS camera) and microwave (RADAR instrument) data, while they are “brown\replacedish” in hyperspectral data (VIMS) (i.e. dark at 1.3, 2 and 5$\mu$m) (Figure 1). Geomorphic analyses show that most of the bedforms are 100-to-1000’s km long linear dunes with a crest-to-crest distance around 3 km (Lorenz et al., 2006, 2010; Le Gall et al., 2011; Savage et al., 2014; Lucas et al., 2014a). Dune heights have been measured as 50-150m high with altimetry and radarclinometry (Barnes et al., 2008; Neish et al., 2010; Mastrogiuseppe et al., 2014). In some SAR swaths, interdune areas are clearly distinguishable from the dunes as exemplified in the Belet sand sea shown in Figure 2. This have been confirmed from hyperspectral observations (Barnes et al., 2008; Bonnefoy et al., 2016). Sources for the organic sand have yet to be identified at the surface and may have even disappeared as these aeolian systems are thought to be old (Rodriguez et al., 2014; Barnes et al., 2015; Brossier et al., 2018). Additionally, radar-bright topographic obstacles called ”inselbergs” are embedded in the sand seas. These features are probably outcrops of the bedrock and thus likely have a texture and composition distinct from that of the dunes. Moreover, their properties may vary across Titan: a ”radar-bright unit” is considered as geomorphic group rather than a single geological unit.

As on Earth (Gao et al., 2015; Lucas et al., 2015) and Mars (Fernandez-Cascales et al., 2018), assessing the surface properties of the sand seas, especially of the dunes and their interdune areas on Titan, is critical because distinguishing between dune origin and nature–that is, their growth mechanisms– will depend upon dune and interdune properties (composition, grain size and thickness of the sediment cover) (Gao et al., 2016; Ping et al., 2017). Moreover, evaluating the properties of sand sea ”inselbergs” is of importance because these geological units could be a proxy for the dune substratum that may outcrop in the interdunes areas. Without a more accurate knowledge of the sediment properties (and availability for saltation) in the interdunes and in the dunes, the interpretations of the observed dune shape and orientation in terms of wind regime and dune activity remain ambiguous. Previous studies based on the infrared VIMS observations (Barnes et al., 2008; Rodriguez et al., 2014; Bonnefoy et al., 2016) and microwave radiometry (Le Gall et al., 2011; Le Gall et al., 2012; Le Gall et al., 2014; Bonnefoy et al., 2016) have revealed a possible difference in composition and/or grain size in the interdunes compared to the dunes. However, a statistical inversion of the data with a physical model has yet to be completed.

In this study, we aim to assess Titan’s surface and subsurface properties within the different equatorial sand seas based on their Cassini RADAR SAR signature. By combining state-of-the-art SAR image despeckling, filtering and classification techniques, modeling of the microwave backscatter, and taking advantage of the large coverage in incidence angles over the sand seas now available, we seek to determine the bulk surface properties of the dunes, interdunes and inselbergs at three major sand seas: Belet, Shangri-La and Fensal (Figure 1). We present the geomorphic classification and microwave backscatter extraction workflow for the SAR data in section 2. We describe our modeling approach and Bayesian inference in section 3. We discuss our results and their implications in section 4 and then present our conclusions.

2 Terrain segmentation and backscattering properties from SAR data: Example from Belet

In this section we describe the data reduction workflow to separately extract the backscatter cross-section ($\sigma^{o}$) from the dunes, interdunes and the inselbergs. We detail below the complete workflow for the Belet sand sea as an example; the same analysis has been applied to the Shangri-La and Fensal sand seas (see the Supporting Information).

2.1 Sand seas as seen by the RADAR SAR imager

In SAR mode, the RADAR experiment provides the highest spatial resolution available of Titan’s surface (down to 300 m), allowing the detection of individual dunes as shown in Figure 2. In addition, some regions have now been observed with a wide range of acquisition geometries in between non-circular Cassini flybys (see Figure 3). In particular, Belet has been imaged with incidence angles ranging between 8∘ to 55∘ and \deletedalmost in the four main cardinal directions (i.e., with a wide range of look directions\deleted) (Figure 3) thanks to the partial overlapping of seven RADAR swaths from Titan’s flybys: T08, T21, T50, T61, T84, T91 and T92. For the sake of clarity, only the first four are shown in Figure 3. Note that Shangri-La and Fensal have been similarly observed (see Appendix A1 and the Supporting Information for the summary of Titan’s observations used in this study). No additional data will be available until the next mission to Titan.

Such a wide range of observation geometries offers a unique opportunity to determine the microwave backscatter behavior (i.e., $\sigma^{o}$ as a function of the incidence angle) and therefore the textural and compositional properties of different terrains over the three sand seas considered. We note that a similar breadth of geometries does not exist for Earth’s deserts. Thus, the scope of our comparison between this work and terrestrial studies in section 4 is limited structural rather than radiometric analogies (shape and orientations of the dunes, grain size sorting between dunes and interdunes, sediment availability in the interdunes, etc.).

2.2 Data reduction and classification

Cassini SAR data suffer from speckle noise, hindering the detection of fine details and quantitative analysis. We therefore use the despeckled SAR data, NLDSAR (Non-Local Denoised SAR) in our work. The efficiency and conservative properties of this non-local denoising is discussed in detail in Lucas et al. (2014a). While significantly improving the SAR image quality and sharpness, which is of great use for pattern identification and detection, we show that the denoising process does not alter the physical values of the backscattering cross-section ($\sigma^{0}$). This is essential when using the data for the inversion with a physical model. \replacedAs exemplified in Figure 4, where the backscatter function over the Belet sand sea (extracted from the T08 swath) is shown, we greatly reduce the data dispersion thanks to the use of the NLDSAR image while preserving its mean values.As shown in Figure 4, we greatly reduce the data dispersion in the backscatter function over the Belet sand sea while preserving the mean by using NLDSAR rather than the raw SAR.

Taking advantage of the NLDSAR data, we then develop a complete workflow in order to automatically extract the backscatter functions of dunes, interdunes and inselbergs. Geomorphic unit segmentation is performed by a split-spectrum analysis over the SAR data inspired by Daily (1983). As discussed in that previous study, low spatial frequency content is due to backscattering signal variations over the surface and subsurface while high frequency content is associated with local topographic features. Low-pass and high-pass filters are hence applied to the SAR data in order to separate low and high spatial frequency information. Additionally, we define a fixed saturation threshold that controls the resulting brightness of the split-spectrum map. This parameter plays a negligible role in the geomorphic unit extraction hereafter. The low frequency content can then be assigned to the Hue axis, the high frequency content to the Value axis, and the fixed saturation threshold to the Saturation level, resulting in a colored image defined in the Hue-Saturation-Value colorimetric domain. Finally, we convert the H-S-V values into the R-G-B domain to create a synthetic-colored image based on the spatial frequency content of the $\sigma^{o}$ amplitude, hereafter referred to as the split-spectrum map. The split-spectrum map created from NLDSAR data greatly facilitates classification, as shown in Figure 5. Finally, the backscattering cross-section values as a function of incidence angle are extracted from the NLDSAR data upon the split-spectrum segmentation. The complete data reduction workflow, including the main steps, is summarized on Figure 6 and has been applied to observations of Belet, Shangri-La and Fensal sand seas (Figure 1, Figure A.1, and Supplementary Information).

2.3 RADAR backscattering cross-section as a function of the incidence angle

We computed the mean within a 3-sigma dispersion of the backscatter cross-section every half-degree of incidence angle. The mean and dispersion are calculated over thousands of measurements for each considered incidence angle, guaranteeing the statistical reliability of the extracted backscatter functions. This is done for each of the 3 geomorphic units (dunes, interdunes and inselbergs) isolated over the Belet sand sea.

Note that the incidence angle of observation is not necessarily the local incidence angle, especially for the dunes and the inselbergs which are not flat surfaces. We expect that the dunes have lee sides with avalanches close to the angle of repose (30*∘) (Neish et al., 2010) based on glints observed in T8 and T61 (Supplementary Information), and that the bright-areas might present local slopes as steep as 75∘* (Tomasko et al., 2005). However, the Cassini RADAR topography products (SARTopo, altimetry, stereo-derived digital elevation maps) cannot precisely resolve local slopes over the sand seas because of limited coverage and/or coarse spatial resolution (Loren et al., 2013; Corlies et al., 2017). Direct information on the local scale topography is therefore missing. We choose to deal with this issue by averaging the backscattering cross-sections into 0.5-degree incidence bins, and retain the resulting values when more than 10,000 pixels per bin is reached. Thus, surface facets with different azimuths and slopes are integrated and topographic effects, if any, are balanced. This is especially true for the dunes, where limited spatial resolution of the SAR images generally impedes the distinction between their two sides even at the scale of a single SAR pixel. Our present study therefore focuses on the ”representative” properties of the dunes, interdunes and inselbergs geomorphic units, averaged over the entire Belet, Shangri-La and Fensal sand seas. Local study at the scale of the dune or interdune objects is beyond the scope of this work but will be the subject of future studies.

The observed backscattering cross-section as a function of the incidence angle over the three sand seas is shown in Figure 7. We emphasize that these are the first backscatter functions ever extracted from denoised SAR images that cover such a large range of incidence angles since the beginning of the Cassini mission. Interestingly, the microwave backscatter functions of the 3 geomorphic units do not overlap, suggesting significant differences in their surface and subsurface properties. The inselbergs show a gently decreasing trend with increasing incidence angles: -0.09 dB/∘ when accounting only for the purely non-specular part of the curve ($\theta>15^{\circ}$) and -0.11dB/∘ for the whole curve. This is an indication of a very large roughness with respect to the Cassini RADAR wavelength of 2.17 cm (see Ulaby et al., 1982) and/or a strong volume scattering in the subsurface as proposed by Paillou et al. (2006). Though lower in backscatter cross-section, the interdune areas also present a similar trend of -0.09 dB/∘ on average for the three dune fields. The dunes, on the other hand, have an even lower value of the backscatter function with a strongly decreasing trend with incidence angle of -0.24 dB/∘ on the non-specular part and -0.18 dB/∘ for the whole curve. This explains their darker appearance compared to the interdunes and inselbergs at all viewing geometries across all Cassini SAR observations. This may reflect a textural difference (i.e., lower roughness), a lower relative permittivity and/or a lower diffuse (i.e., roughness and/or volume scattering) component (Ulaby et al., 1982). The overall flatness of the curves for the three units between 10∘ and 50∘ (e.g., 20∘ and 50∘ for the dunes) suggests a significant diffuse component due to roughness or/and volume scattering. These will be further quantified in the following section by the use of microwave backscattering models.

3 Microwave backscatter modeling and inversion results

In this section we first discuss the fundamentals of microwave backscatter modeling and then explain our choices for the subsequent inversion of the surface and subsurface properties of the considered terrains \replacedfromwith Bayesian inference.

3.1 Microwave backscattering forward models

Microwave backscatter modeling has been widely used for terrestrial and planetary applications. We base our analysis on previous studies that have developed models that depend on surface roughness and relative permittivity (represented as $\epsilon=\epsilon^{\prime}-j\epsilon"$), accounting for coherent and single-scattering parts of the signal (Fung, 1992, 1994; Ulaby et al., 1982; Paillou et al., 2006).

The surface roughness is defined by two empirical statistical descriptors: the root-mean-square height ($\xi$), quantifying the vertical variations of the surface, and the correlation length ($\zeta$), characterizing the typical horizontal size of surface bumps. We usually define surface roughness with a Gaussian height distribution and exponential autocovariance function $\rho_{N}$:

[TABLE]

$z$ being the elevation with $\bar{z}$ being its mean, $n$ the number of measurements along the profile, and

[TABLE]

with \big{<}\big{>} being the averaging operator and $\delta$ being the measured distance from a considered point $x_{i}$. The correlation length is given by

[TABLE]

We then choose the surface backscatter term. For the range of very smooth to very rough terrain, 1-layer Integral Equation Method (IEM), PO (Physical Optics) and GO (Geometrical Optics) surface scattering models are appropriate (see Appendix A2 for further details). Here we use the validity domains defined by equation (11) following Paillou et al. (2006) (Appendix Figure A.3). It is important to recall a few limitations of the considered models.

$\bullet$ Note that the valid domains for the surface term (see $\sigma^{0}_{Spp}$ in Appendix A3) differ slightly throughout the literature (e.g., Paillou et al., 2006; Fung, 1992; Zribi and Dechambre, 2002). Moreover, it has been shown that for terrestrial terrains, models are not always able to reproduce ground truth observations in their theoretical domain of validity due to some misconception and/or poorly defined surface roughness field measurements (Zribi et al., 1997, 2000; Zribi and Dechambre, 2002) (See Appendix A3).

$\bullet$ As mentioned above, the non-coherent surface scattering depends on the two statistical descriptors of the surface roughness ($\xi$ and $\zeta$), which, as defined, are bound to one another (see Appendix A2). Consequently, they cannot be inverted independently.

$\bullet$ The ability of these descriptors to characterize natural surfaces is, in practice, not straightforward, as they depend on the scale at which they are defined (i.e., the length of the roughness profile) as shown by Baghdadi et al. (2000). Moreover, the interdependence of the two statistical descriptors is itself a function of this scale and the nature of the terrain (i.e., rough vs. smooth) (e.g., Baghdadi et al., 2000; Zribi and Dechambre, 2002; Bretar et al., 2013). It has been shown that it is more relevant to consider the ratio of the two descriptors in order to characterize the roughness of a natural surface: the root-mean-square slope ($s=\xi/\zeta$). Note that only the GO model depends directly on the surface slope $s$ and hence provides a better representation of the surface roughness. Although previous studies suggest that the ratio $Z=\xi^{2}/\zeta$ may be a more satisfying descriptor of roughness, it cannot be directly related to the backscattering modeling as described above and hence is not considered in our study.

The Snell-Descartes laws assume flat interfaces and that deviations from planarity do not exceed a wavelength fraction. A surface is therefore considered to be rough if the mean quadratic deviation of surface irregularities satisfies the Rayleigh quality criterion (Ulaby et al., 1982):

[TABLE]

$\lambda$ being the wavelength and $\theta$ the incidence angle. Based on the Cassini RADAR wavelength of 2.17 cm, the validity domain for each model (Appendix A.3), the roughness of natural terrains (Appendix Figure A.5) and previous works from Paillou et al. (2014), all Titan equatorial terrains appear rough and fall into the GO model domain of validity. We confirm this by Monte Carlo simulations accounting for 1-layer models with no model parity: ¿90% of our reduced $\chi^{2}$-tests best fits are obtained with the GO model (Appendix A.2). Consequently, subsequent analysis is only done with the GO model for which the backscattering function reads:

[TABLE]

with $m=\sqrt{2}\xi/\zeta$ and $\Gamma$($\theta$ = 0) being the Fresnel reflectivity at normal incidence angle (i.e., $\theta=0$). Recently, Paillou et al. (2006, 2014) applied these models to the Cassini SAR data and emphasized i) that Titan’s surface is rough, with respect to the RADAR’s wavelength confirming our preliminary analysis, ii) the importance of volume scattering from the subsurface. These findings are also supported by scatterometry and radiometry observations (Zebker et al., 2009; Janssen et al., 2016). To account for the volume scattering contribution, we also add the following term:

[TABLE]

where $T_{ij}$ is the Fresnel transmission coefficient between two media (i.e., the air and ground), and $\theta_{t}$ is the transmitted incidence angle\replaced, and. $\tau$ is the optical depth defined as $\tau$ = 1/(1-$a$), with $a$ being the microwave albedo of the first interface (i.e, the surface in case of 1-layer model) which ranges between 0 and 1. Note that the term $\sigma^{0}_{Vhh}$ accounts for the totality of the diffuse part of the backscattering signal. That is, volume scattering and/or surface multiple-scattering processes can not be differentiated.

The overall modeled backscattering coefficient is then obtained by evaluating $\sigma^{0}=\sigma^{0}_{Shh}+\sigma^{0}_{Vhh}$ (in physical units and then converted into dB) with varying $\xi$, $\zeta$, $\epsilon$ and $a$.

Although the Cassini RADAR microwave signal can penetrate a few 10’s of cm beneath the surface, only 1-layer models are considered in this study (Paillou et al., 2006). The number of unknowns in the case of the +2-layer models becomes high and would require additional constraints on the physical properties of Titan’s surface and subsurface that we do not yet have. In Titan’s case, therefore, it would provide no further insights into the physical properties of Titan’s geomorphic units (Paillou et al., 2006).

3.2 Sensitivity analysis over synthetic tests from Bayesian inference

The non-linearity and non-homogeneity of the modeling make the physical parameter inversion non-trivial. As we aim at being as quantitative as possible, we performed a sensitivity analysis with the variance-based Sobol’s method which quantifies the amount of variance that each parameter contributes to the unconditional variance of the model output (Sobol’, 1993). These amounts, caused either by a single parameter or by the interaction of two and more parameters, are expressed as sensitivity indices (SI), which represent fractions of the unconditional model output variance. The first-order SI are the direct variance-based measure of the contribution to the output variance of each parameters. The case of a very rough surface, including the volume scattering term (i.e., GO model + vol. scat.) is shown on Figure 8a for all physical parameters considered. Our sensitivity analysis quantitatively demonstrates that the effect of surface RMS slope ($\sqrt{2}s$, with $s=\xi/\zeta$) dominates at low incidence angles, while the contribution of the albedo $a$ dominates at higher angles. Additionally, the contribution of the real part of the relative permittivity $\epsilon^{\prime}$ remains low ($<$20%). The GO model is insensitive to its imaginary part. However, the low uncertainties of the first order SI emphasize that even low-contribution parameters can still be assessed as long as the incidence coverage is sufficient.

In order to compute the marginal posterior probability of each physical parameter, we use a Bayesian statistical model and fitting algorithms based on the Markov Chain Monte Carlo inversion (MCMC) method (Salvatier et al., 2016). We assess the capability of the method to correctly retrieve the physical parameters\replaced, synthetic tests are doneusing synthetic tests as illustrated in Figure 8b. In this example, the backscatter function is computed from the GO+Vol. model with $\epsilon^{\prime}=1.55$, $s=\xi/\zeta=0.10$, and $a=0.30$ for a variety of incidence angles randomly drawn with significant gaps between 30∘ and 50∘. Gaussian noise is independently added on the mean values for each incidence angle with a standard deviation of 0.3 dB which corresponds to the $\sqrt{N_{RMS}}$ of the derived noise obtained from NLDSAR technique (Lucas et al., 2014a). Uncertainties are simulated from a uniform distribution with 1-$\sigma$ at 0.6 dB. We therefore obtain realistic synthetic data that mimic the actual data from Cassini SAR (Figure 8 and Appendix A3). As shown in the example in Figure 8b, true parameters are accurately assessed by our inversion scheme. This method also presents the advantage of providing confidence intervals for each of the retrieved parameters.

3.3 Titan’s equatorial regions microwave properties from Bayesian inference

In order to characterize the physical parameters $s$, $\epsilon$’ and $a$ which reflect the surface roughness (i.e., RMS slope), the composition of the terrain ($\epsilon$’ may also depends on porosity), as well as the subsurface contribution respectively, similarly to the synthetic tests, we perform a Bayesian inference over the GO+Vol. model based on advanced Markov Chain Monte Carlo method (Salvatier et al., 2016). This analysis has been conducted with 50$\times$106 runs of the backscatter model for $\epsilon^{\prime}\in$ [1,5], which corresponds to the expected range of values for material relevant to Titan’s surface at 90 K and 13 GHz (Paillou et al., 2008). $s=\xi/\zeta\in$ [0.001,1] is based on terrestrial values (see Appendix A1), and $a\in[0.0,1.0]$. Resulting optimal fits and 95% confidence intervals (CI) are shown for the three considered units of the three sand seas in Figure 9. In all cases, the 95% CI falls into the observed dispersion, even when the incidence coverage is small. The mean observed data are correctly retrieved by our simulations (i.e., black lines on Figure 9). As discussed below, the synthetic tests support our results on the real observations. Indeed, these tests consist of both exploring the sensitivity to observation properties (incidence angle) and to the standard deviation and noise function derived from Lucas et al. (2014a). Additionally, the spatial distributions (somewhat randomly distributed) of the physical parameters have been considered in theses tests. As illustrated in Fig. 8 and Fig. A.8, and Appendix A4, the synthetic tests show that our method does not breakdown when accounting for the actual data properties (coverage, standard deviation, noise, etc.).

The \deletedresulting marginal posterior probabilities computed from our Bayesian analysis for each physical parameter are summarized in Figure 10. We obtain nearly Gaussian posterior probabilities for the dune units for each parameter ($\epsilon^{\prime}$, $\xi/\zeta$, and $a$) of the three sand seas. Overall, the resulting distributions are very similar and emphasize the homogeneity of the dune material over the three studied regions.

The results are more subtle in the case of interdune and bright area units. The shapes of their respective marginal posterior probabilities indicate that for all of the considered physical parameters, the interdunes and bright areas are made of more complex media, possibly with mixed contributions. These distributions also suggest an important variability of surface properties across Titan at the global scale. It is worth noting that the radar-bright areas are extracted from only inselbergs in Belet; without interdunes this population is therefore inherently more homogenous. Our results suggest that the radar-bright areas in Shangri-La and Fensal are more heterogeneous. Note that both bound Xanadu (Figure 1). It is unclear at this stage how this proximity may affect our regions of interest by, for example, influencing sediment transport.

For the inselbergs, the model is not capable of reproducing the observations (see ”no-ampVpp” curves in Figure 9). Indeed, the albedo reaches its saturation value of 1. Hence an amplification of the diffuse part of the backscattering signal is required in order to fit the observations (”ampVpp” on Figure 9). This amplification of the diffuse part is in agreement with the findings for the radar-bright region of Xanadu from previous studies (Janssen et al., 2011, 2016). Accounting for this effect by multiplying the subsurface term by a factor $>$ 2 (i.e., 3) after Janssen et al. (2016), we go beyond the coherent backscattering effect and hence obtain the curve ”ampVpp” in Figure 10. As shown, this correction remains to be explained by some mechanism, but it allows us to partly retrieve a Gaussian distribution of the labedo $a$. Note that these bright areas are in many ways similar to the Xanadu region. They also show similar backscatter function values compared to Enceladus for which a few SAR observations have been made (Le Gall et al., 2017). In the absence of any sediment cover, icy bedrock may present high to extreme surface roughness due to high internal friction, in particular on steady surface (Kietzig et al., 2010). Such extreme surface roughness has been shown to enhance the backscatter beyond the coherent backscatter effect in lunar regolith (Campbell, 2012). Another explanation would come from a strong subsurface terms as proposed by the two layer model of Paillou et al. (2006). As stated above, a two layer model is beyond the scope of this work due to the untenable number of unknowns.

The real part of the relative permittivity ($\epsilon$’) significantly differs between the three considered geomorphic units. For the Shangri-La sand sea, the dune’s distribution shape is Gaussian with a mean value around 2. The interdune areas, however, show a non-Gaussian distribution with peaks at 1.2 and $\sim$1.7 and the bright areas are pseudo-Gaussian with a mean value around $\sim$3 with large uncertainties. By comparing to laboratory measurements (Paillou et al., 2008), the derived relative permittivities reflect a homogeneous, organic-based composition for the dunes; the relative permittivity for bright terrains is more compatible with a water-ice/ammonia mixture. Therefore, the wide bi-modal distribution of the interdunes might be explained by the presence of a mixed composition between these two end-members. Thus, when the inversion tries to fit with this non-homogenous composition\replaced and failed in finding, it fails to find one particular value. Note that this is consistent with the presence in the interdunes of a thin and/or partial sediment cover over an icy-bedrock as discussed previously in (Rodriguez et al., 2014).

Contrary to other geomorphic units, the surface roughness assessed through the ratio ($\xi/\zeta$) shows similarity between the dune units of the three sand seas. All three bimodal distributions have peaks at $\sim$0.05 and $\sim$0.25 which give RMS slopes of 0.07 and 0.35 respectively (where RMS slope is equal to $\sqrt{2}{}\xi/\zeta$). The implications of this finding are discussed in the next section.

Lastly, the inferred contribution of the volume scattering mainly depends on the microwave albedo $a$. Low values of $a$ correspond to low volumetric energy return due to a small penetration depth (resulting from high absorptivity of the surface material) of the medium \deletedalong the depth of penetration. High values of $a$ correspond to higher volumetric return due to lower absorptivity and/or volume heterogeneities caused by the presence of a bedrock and/or subsurface defects. Over the dunes, we find relatively low values of $a$ (0.2-0.3) indicating that they strongly attenuate the radar signal (i.e., they appear dark in amplitude data). The albedos of interdunes (and bright areas), however, contain a strong subsurface component (i.e., values of 0.5-1). This implies a significant bedrock contribution beneath the organic-covered terrrains, in agreement with the results from the other two parameter inversions as well as observations by microwave radiometry (Le Gall et al., 2011; Le Gall et al., 2012) (low resolution) and infrared spectroscopy (Barnes et al., 2008; Bonnefoy et al., 2016) (resolved).

To summarize, our Bayesian inferences show a clear distinction between dunes, interdunes and inselbergs in terms of roughness and composition. Ourresults provide more evidence for the global homogeneity of the dune units and show that interdunes and bright areas are composed of more complex and varying types of terrains, with possibly strong regional variations \explainhow can the findings be ”possibly strong”? It seems like it should be one or the other.. Interdunes systematically present a wider surface roughness distribution (with higher roughness), a widely distributed permittivity and a stronger subsurface contribution than the dunes. We discuss what this implies for the aeolian system morphodynamics in the following section.

4 Discussion

We are approaching the problem with both expertise in microwave modeling and aeolian sediment transport. Hence, we discuss here the connections between results of the MCMC and interpretation in terms of aeolian morphodynamics.

4.1 Interpretations on geomorphic units properties

As shown dunes in Fensal have slightly higher relative permittivity and albedo than in Belet and in Shangri-La. Note that the incidence angle coverage goes at lower values (i.e., almost 5*∘*) for Fensal compared to the other regions. According to our sensitivity analysis, the roughness contribution is particularly strong in this incidence domain. So, one explanation would be that this is the only region where the model can very well constrain the roughness. However, this difference may also reflect actual variations due to latitudinal and altimetric controls of Titan’s dune field morphometry (Le Gall et al., 2012). Fensal is indeed at both higher latitude and altitude than Belet and Shangri-La. Sediment availability is expected to decrease when latitude and altitude increase. The relative permittivity inversions show that i) dunes are compatible with a very homogeneous sand made of porous organics, ii) bright terrains are a mixture of contaminated water-ice with tholins in different ratios, and iii) the backscatter function of the interdunes may reflect spatial variations with contributions from both organics and icy bedrock. Our results thus provide a direct confirmation of what was previously suggested from infrared analysis (Barnes et al., 2008; Soderblom et al., 2007) and passive microwave observations (Janssen et al., 2016; Le Gall et al., 2011) compared to laboratory characterization measurements (Paillou et al., 2008). Our derived values for the roughness also provide interesting insights. The distributions have a systematic best-fit with observations despite differences in data coverage (Figure 7).

While we cannot directly link the actual distribution of the roughness with our results over the dunes, compared to other units, reflect a clear non-homogeneous distribution. Note that the synthetic tests emphasize that both model and method are stable even for spatially distributed values over the three main parameters. Moreover, dunes on Earth and on Mars, are represented by smooth terrains (at micro-wave wavelengths) with secondary patterns, such as ripples and avalanches on the lee side. Therefore, we expect that these secondary patterns that should be present on Titan as well – otherwise these features would not be dunes. The explanation for non-Gaussian distribution after MCMC for the dunes is here presented as the ability of the method to detect the systematic effect (i.e., for the three regions) of such patterns in the roughness component of the microwave backscatter function over Titans’ dunes. Note that such distributions are not consistent with the other, non-dune units across regions of interest. Because the distributions we obtained show a systematic pattern across the three regions of interest (again accounting for differences in the observations described above), we suggest that the patterns may reflect distinct surface roughness contributions. The lowest surface roughness contribution would correspond to a smooth surface where vertical roughness is one order of magnitude smaller than horizontal. The highest would be compatible with second order bedforms of higher aspect ratio, such as ripples and/or avalanches. Uncertainties span over the obtained distributions and so do not allow us to be more quantitative on the actual roughness.

Interdunes and bright-areas exhibit wider surface roughness distributions for Shangri-La and Belet. This may be due to multi-scale roughness variations compatible with either mountainous/dissected terrains or sedimentary plains made of a non-homogeneous and widely distributed granularity, with significant contributions from higher scale roughnesses. This distinct difference between the roughness of the dunes and the interdunes is a critical result. It might correspond to differences in sediment size and/or variations of thicknesses of the organic cover in the interdunes, in agreement with recent analysis in the IR spectra from VIMS instrument (Bonnefoy et al., 2016). Fensal roughness distributions are, on the other hand, much narrower with close mean values. This may be indicative that Fensal’s dunes, interdunes and inselbergs differ in composition (as at Shangri-La and Belet) but are similar in terms of texture (which is not the case at the other two dune fields).

To summarize our interpretations, the dunes are organic-rich smooth terrain, potentially with detected secondary bedforms, while interdunes are water-ice rich rough bedrock partly covered by a shallow layer of organic sediments.

4.2 Implications on morphodynamics

All studies of dune morphologies and interactions with topography agree that Titan’s linear and barchanoid dunes propagate eastward (Lorenz et al., 2006; Radebaugh et al., 2008; Rodriguez et al., 2014; Lucas et al., 2014b; McDonald et al., 2016). Due to our limited knowledge on the current climatic conditions, only recent Global Circulation Models (GCM) have succeeded in generating surface winds with eastward mean direction and sufficient strength to allow sediment transport (Lebonnois et al., 2012; Tokano, 2010; Ewing et al., 2015; Charnay et al., 2015; McDonald et al., 2016). Globally averaged surface winds on Titan are predicted by GCMs to be primarily bi-modal, blowing seasonally to the southwest or the northwest with an average speed rarely exceeding $\sim$1 m/s near the surface (Tokano, 2010; Lebonnois et al., 2012; McDonald et al., 2016). Consequently, while it is commonly accepted that these aeolian systems are mature and still currently active with a  50, 000-year timescale for Titan’s dunes to reorient following global climatic changes (Ewing et al., 2015), there is a long-standing debate on whether the dune orientations reflect past or current wind conditions at Titan’s surface (Ewing et al., 2015; Lucas et al., 2014b; Charnay et al., 2015; McDonald et al., 2016). The current wind conditions are a priori incompatible with observed dune migration to the east.

In order to reconcile geomorphic observations with wind conditions, two opposing dune evolution models have been proposed: i) reorientation of dune crests resulting from variations in sediment availability, wind direction and strength on long-term climate cycles controlled by variations of Saturn’s eccentricity (Croll-Milankovitch cycles) (Le Gall et al., 2012; Ewing et al., 2015; McDonald et al., 2016) and ii) equinoctial gusts and/or storms such as those observed in 2010 (Turtle et al., 2011a) generating winds that drive strong eastward sediment fluxes under the current conditions (Tokano, 2010; Lucas et al., 2014b; Charnay et al., 2015). Note that the second hypothesis succeeds in accurately explaining the dune orientation, direction of migration, as well as the equatorial dune confinement (Lucas et al., 2014b). In either case, aeolian transport is clearly the dominant transport mechanism of sediment in the equatorial zones (Malaska et al., 2016). Isolated dunes evolve by recycling their own sediment at the crest. Hence, their sedimentary structure reflects the wind regime responsible for their formation and evolution. However, this effect is limited to a short period of time. At the scale of a vast sand sea, morphologies integrate climatic cycles punctuated by alternations between arid and humid periods.

Therefore, it is challenging to determine if a sand sea is in equilibrium with contemporary wind regimes or whether it is in a transient state reflecting different wind regimes integrated over time (Ewing et al., 2015; Lucas et al., 2014b; McDonald et al., 2016). When the sediment bed is a thick layer of particles that can be entrained by the wind, dune growth can be described as a bed instability selecting the alignment for which the normal-to-crest component of transport is maximized (Rubin and Hunter, 1987; Ping et al., 2014; Gadal et al., 2018). It gives rise to the emergence of a regular dune pattern where topographic highs (crests) and lows (troughs) keep the same composition and texture. In contrast to this transport-limited condition, recent works have shown that dunes may also elongate on a non-erodible surface in the direction of the resultant sediment flux at the crest (Reffet et al., 2010; Courrech du Pont et al., 2014; Gao et al., 2015; Lucas et al., 2015). These finger-like structures occur under multidirectional wind regimes in zones of low sand availability or where coarse grains form an armor layer as a result of size-segregation effects (e.g., selective transport and/or avalanche as described in Gao et al. (2016)). In all these cases, dunes composed of mobile particles can be easily distinguished from the interdune areas (bedrock, desert pavement). The sediment reservoir therefore critically controls the dune growth mechanism and the subsequent dune and surface properties. In regards to these considerations, we hence propose that our results on the shallow surface of Titan’s sand seas (dune and interdune areas) can provide new insight on the local dune morphodynamics.

On Earth, dunes elongating over tens of kilometers are the prominent geomorphic features in the Ténéré desert in Niger (Lucas et al., 2015) and in the Kumtagh desert in China (Ping et al., 2017), where they form isolated linear ridges aligned with the resultant transport direction (Figure 11). These longitudinal dunes are regularly indented by transverse secondary bedforms, which may locally eject barchans. Ping et al. (2017) have shown that non-symmetric multidirectional wind regimes like those currently measured in these deserts produce the observed dune patterns. Indeed, the coexistence of the two dune growth mechanisms can naturally explain how finger-like structures elongate in a specific direction and generate superimposed bedforms with a different orientation as the bed instability develops on their flanks. The instability is not observed in the interdune areas because the bed surfaces exhibit armor layers composed of coarse-grained material and zibars (Figure 11). The transport rate of uniform fine grains along the crest can then be derived from the local wind data to explain the elongation rate of dunes as well as the migration rate of secondary bedforms (Lucas et al., 2015; Ping et al., 2017).

Although there is no established link between the grain size and the surface roughness at the RADAR wavelength, our infrared and microwave analyses suggest that the surface layers of Titan’s dune and interdune areas are different in composition and texture. On many (if not all) dune fields on Earth, such differences between dunes and interdunes is the result of granular sorting causing armoring effect in aeolian systems. A grain size ratio of 5 to 10 and a surface roughness ratio of 10 between the interdune areas and the dunes are typically observed (Figure 11). Accounting for the above considerations and as we find a significantly differences in the RMS slopes, we propose that the sediment beds in the interdune areas have a mean grain-size larger than on the dune bodies. In Titan’s equatorial belt, dunes essentially composed of fine sediment may thus have developed over a coarser/rougher sediment layer as a result of a size-segregation process.

Such an interpretation of the compositional and textural properties of Titan’s dune fields supports the elongating growth mechanism proposed by Lucas et al. (2014b). Considering this scenario, sediment transport occurs mainly along the linear dunes and not in the interdune areas covered with a pavement of coarse particles. The dynamics of the size-segregation process could also explain the origin of the dune fields in the absence of apparent source of sediment in Titan’s equatorial belt. Indeed, the mobile particles that compose the active dunes may be stored under the armor bed and occasionally set in motion by the strong eastward winds associated with equinoctial storms (Tokano, 2010; Charnay et al., 2015). In addition, Titan’s dune ”defects” depicted in Ewing et al. (2015) as evidence of long-term dune re-orientation from ancient to modern winds are in fact similar to the secondary patterns of raked linear dunes documented in the Kumtagh desert (Figures 11 and Ping et al. (2017)). These smaller scale dune features could then be the natural byproduct of Titan’s dune morphodynamics under the current wind regime. In addition to the other geomorphic processes preventing transport (Rubin and Hesp, 2009; Harper et al., 2017), the granular segregation and/or the armoring effects might contribute to the morphodynamics of the linear dunes on Titan similarly to terrestrial cases (Gao et al., 2016).

5 Conclusions

Microwave backscatter functions have been derived from Cassini SAR data, processed (notably corrected from the speckle noise) and analyzed over the three major sand seas in the equatorial belt of Titan. First, we show that the backscatter function from the interdunes is clearly distinguishable from that of the dunes. This demonstrates a clear difference in composition, attenuation, and surface roughness. Our microwave analysis is fully consistent with hyperspectral observations from the VIMS instrument. As demonstrated by a Markov chain Monte Carlo inversion, interdunes have a higher surface roughness, a distinct relative permittivity and a stronger subsurface contribution than the dunes. These properties may reflect a thin and/or partial sediment bed of a mixture of fine and coarser grains over an icy crust interspersed with organics. Our results may confirm that the dunes are primary composed of organics and very homogeneous at the global scale. The bimodal surface roughness distribution we derived for the dunes is also compatible with the presence of secondary bedforms, such as ripples or avalanches.

The derived physical properties of the three major Titanian sand seas suggest that interdunes may act as an armored surface composed of large grains over which the linear dunes may develop. As coarser grains would need higher shear stresses to be moved, this implies that dunes and interdunes significantly differ in their respective contribution to the total sediment flux: the sand of the dunes largely dominates the flux. Such a scenario would be similar to what is observed in Niger or China where linear dunes elongate along the resultant drift direction under a multi-directional wind regime over a sedimentary bed dominated by pebbles and boulders which are rarely mobilized by the winds. The resulting morphodynamics is therefore compatible with the elongating growth mechanism proposed by Lucas et al. (2014b).

Indeed, unless the current wind regime at Titan’s surface is bi-modal and purely symmetric – which is unlikely due to Titan’s non-circular orbit, topography, and albedo variations – one can expect a different orientation of the dune crests based on the sediment supply of the interdunes: (1) depleted and/or coarser-grained interdunes, or (2) densely covered interdunes with the same granularity of the dunes. Furthermore, the state of the interdunes contributes to linear dunes growth as observed in the Rub’ al Khali and in the northern hemisphere of Mars. Additionally, the different dune shapes observed on Titan (e.g., raked linear dunes, star dunes, asymmetric barchans, see Radebaugh (2013); Lucas et al. (2014b); Ewing et al. (2015)) clearly suggest the presence of a complex wind regime associated with a limited sediment supply at Titan’s equatorial region; these observations are also compatible with our results.

Consequently, sand seas on Titan appear to be currently evolving under a complex wind regime, with dune pattern coarsening and interdune armoring similar to their terrestrial analogues either in Saudi Arabia, Niger or China. Future investigation, further investigation into the orientations of potential secondary bedforms (i.e., ripples or raked-shaped linear dunes) and their sub-sequent secondary winds must account for local azimuthal variation of the backscatter function and local slope effects will be necessary. The methods presented here could also allow us to quantitatively retrieve the surface properties of different geomorphic units (e.g. plains, badlands, mountainous, fluvial, and lacustrine terrains) with the highest spatial resolution possible with Cassini data. Such results would enable further investigation into Titan’s geological and climatic processes and inform future in situ missions to Titan such as Dragonfly, which is target to explore these equatorial regions (Turtle et al., 2017). Based on our results, the Shangri-la sand sea offers the opportunity to explore dune and interdune areas of distinct composition and texture, thereby providing Dragonfly a representative sampling of Titan’s diverse sedimentary and geomorphic transport processes.

Acknowledgements.

No author declares conflicts of interests. Cassini SAR data are available on the PDS (https://pds-imaging.jpl.nasa.gov/portal/cassini_mission.html). Due to the large file size, NLDSAR data set is available upon request by contacting the corresponding author. Part of the dataset is accessible from Zenodo doi:10.5281/zenodo.528545. The authors thank Jani Radebaugh and Catherine Neish as well as anonymous reviewers and the JGR-Planets associate editor for their constructive evaluation that helped to substantially improve the manuscript. Authors also thank Stéphane Jacquemoud, Sebastien Labarre, Léon Sylvain and Emma Ayçoberry for their fruitful inputs that improved this work. Finally, the authors also thank the Cassini Radar Science Team for constructive feedbacks during this research. The authors thank Y. Callot and G. Steinmetz for providing field pictures from Niger and China, respectively. The authors acknowledge the financial support of the UnivEarthS Labex program at Sorbonne Paris Cité (ANR-10-LABX-0023 and ANR-11-IDEX-0005-02), the ANR EXODUNES, Space Campus from Université Paris Diderot, the Nvidia Academic Research Support program as well as the CNES and Institut Universitaire de France. SMM was supported by NASA Outer Planets Research Program NNX14AR23G/118460.

Appendix

A1. Cassini SAR data summary

The overall incidence coverage for each respective swath used in this work is shown in Fig. A.1.

A2. Microwave backscatter models and their domain of validity

From very smooth to very rough terrain, 1-layer Integral Equation Method (IEM), PO (Physical Optics) and GO (Geometrical Optics) surface scattering models are appropriate and read:

[TABLE]

where $k$ is the wave number, and $W^{n}$ the Fourier transform of the nth power of the surface correlation function, which depends on $\zeta$, and with

[TABLE]

with $f_{hh}$ and $F_{hh}$ being the coefficients of the Kirchhoff and complementary fields (Fung, 1992); and

[TABLE]

with $J_{0}$, $\Gamma_{hh}$ and $\rho(x)=\exp(-x^{2}/\zeta^{2})$ being respectively the zeroth-order Bessel function of the first kind, the Fresnel reflectivity and the Gaussian surface autocorrelation function; and

[TABLE]

with $m=\sqrt{2}\xi/\zeta$ and $\Gamma(\theta=0)$ being the Fresnel reflectivity at normal incidence. Note that the coherent part of the backscatter is correctly simulated with the GO model which is therefore valid at all incidences. The non-coherent surface scattering model is selected as a function of the roughness domain relative to the RADAR wavelength. From smooth to very rough terrain, IEM, POM and GOM are alternatively considered for the modeling of the surface backscattering coefficient $\sigma^{0}_{Shh}$ (e.g., Paillou et al., 2006):

[TABLE]

We checked that the three $\sigma^{0}_{Shh}$ models provide the continuity of the solutions at the limits of each validity domain (Figure A.3). We run 5.105 Monte Carlo simulations accounting for 1-layer models with no model parity and found that most (¿90%) of our reduced $\chi^{2}$-tests best fits are obtained with the GOM for the 3 geomorphic units considered. We hence obtain realistic synthetic data that mimic the actual data from Cassini SAR (Figure A.2). Note that this is compatible with previous works (Paillou et al., 2014).

A3. Surface roughness of natural terrains

As discussed in the main text, statistical descriptors of natural surface roughness depends on the scale and the nature of the considered terrain. Here are combined RMS height ($\xi$) and correlation length ($\zeta$) from literature and photogrammetric surveys on agriculture terrains performed by the corresponding author (Figure A.5). These terrains exemplify the range of surface roughness observed in nature and are to be compared with respect to the values retrieved in this study. Note that there is no direct link between the nature of the terrain and its surface roughness as the latter is controlled by different processes that sculpt the topography. Also, multi-scale effects of the surface roughness are known to play a significant role in the radiometry (Labarre et al., 2017). Further investigation into such multi-scale effects is necessary but beyond the scope of this paper.

A4. Synthetic tests on surface physical properties and subsequent microwave backscatter cross-section

Synthetic tests are done in order to show the effect of spatial variations of the physical properties that control the microwave backscatter signal ($\sigma^{0}$). Synthetic maps are generated (Figure A.6) and the resulting $\sigma^{0}$. Note that shadowing and layer over effects are not considered when synthesizing the SAR image. Multiplicative speckle noise with a gamma distribution has been considered based on noise estimation performed in Lucas et al. (2014a).