Mohsen Thameur Hajlaoui
email: thameur.hajlaoui@emt.inrs.ca

Elain Hajlaoui

Michaël Josse

Essebti Dhahri

Alain Pignolet

Mohsen Elain Hajlaoui

Tetragonal tungsten bronze/barium hexaferrite room-temperature multiferroic composite ceramics

Keywords: composite, room temperature multiferroic, electromechanical properties, polar domains, magnetic domains

come

Introduction

Recently, multiferroic/multifunctional materials have attracted the interest of scientists due to their rich properties and their potential use in a wide range of advanced applications [START_REF] Spaldin | Multiferroics : Past , present , and future feature[END_REF][START_REF] Lu | Multiferroic oxide thin films and heterostructures[END_REF][START_REF] Ramesh | Multiferroics: progress and prospects in thin films[END_REF][START_REF] Eerenstein | Multiferroic and magnetoelectric materials[END_REF][START_REF] Nan | Multiferroic magnetoelectric composites: Historical perspective, status, and future directions[END_REF]. In particular, the coexistence of ferroelectric and ferromagnetic properties in these materials makes them promising candidates for data storage devices [START_REF] Spaldin | Advances in magnetoelectric multiferroics[END_REF][START_REF] Baek | Ferroelastic switching for nanoscale non-volatile magnetoelectric devices[END_REF], including new generations of non-volatile integrated memories [START_REF] Shen | A multilevel nonvolatile magnetoelectric memory[END_REF][START_REF] Roy | Multiferroic memories[END_REF][START_REF] Scott | Multiferroic memories[END_REF]. In order to overcome the issue of the scarcity of single phase multiferroic materials, especially at the ambient temperature and above, several strategies have been adopted in order to develop new multiferroic materials at room temperature [START_REF] Nechache | Epitaxial thin films of multiferroic Bi2FeCrO6 with B-site cationic order[END_REF][START_REF] Sanchez | Room-temperature single phase multiferroic magnetoelectrics: Pb(Fe, M) x(Zr,Ti)(1-x)O3 [M = Ta, Nb[END_REF][START_REF] Zhao | Near room-temperature multiferroic materials with tunable ferromagnetic and electrical properties[END_REF][START_REF] Scott | Room-temperature multiferroic magnetoelectrics[END_REF]. One of the most efficient is to synthesize composites materials formed by ferroelectric and ferromagnetic phases, which allows the use of a large variety of ferroelectric/piezoelectric and magnetic material couples, having both good ferroelectric and good ferromagnetic properties [START_REF] Nan | Multiferroic magnetoelectric composites: Historical perspective, status, and future directions[END_REF][START_REF] Jayachandran | Solutions for maximum coupling in multiferroic magnetoelectric composites by material design[END_REF][START_REF] Ma | Recent progress in multiferroic magnetoelectric composites: From bulk to thin films[END_REF][START_REF] Ai | Synthesis and magnetoelectric properties of multiferroic composites of lead lanthanum zirconate titanate and mesoporous cobalt ferrite[END_REF][START_REF] Gao | A comparative study of the dielectric, ferroelectric and anomalous magnetic properties of Mn0.5Mg0.5Fe2O4/Ba0.8Sr0.2Ti0.9Zr0.1O3 composite ceramics[END_REF]. An important issue with composites is the chemical compatibility of the targeted phases, and the formation of secondary phases (or of interphases) can be complex and challenging [START_REF] Laguta | Crystal structure transformations induced by surface stresses in BaTiO3 and BaTiO3@SiO2 nanoparticles and ceramics[END_REF][START_REF] Chung | Evidence of diffusion at BaTiO3/silicon interfaces[END_REF].

Tetragonal tungsten bronzes (TTB) are a class of materials having a crystal structure consisting of a framework of corner-sharing octahedra connected in such a way that different channels with different shapes (triangular, square and pentagonal channels) are formed. These channels are all oriented along the short crystallographic c-axis, and can be filled, partially filled or empty. Such a crystal structure with a large number of cationic and anionic sites allows for the synthesis of large variety of composites and dopings, resulting in a fine control and a large enhancement of the functional properties of the synthesized material [START_REF] Josse | The Ba2LnFeNb4O15 "tetragonal tungsten bronze": Towards RT composite multiferroics[END_REF][START_REF] Aamlid | The Effect of Cation Disorder on Ferroelectric Properties of SrxBa1-xNb2O6 Tungsten Bronzes[END_REF][START_REF] Josse | Original crystal-chemical behaviors in (Ba,Sr)2Ln(Fe,Nb,Ta)5O15 tetragonal tungsten bronze: Anion-driven properties evidenced by cationic substitutions[END_REF][START_REF] Yang | Effects of A-site cations on the electrical behaviors in (Sr1-xCax)2.1Na0.8Nb5O15 tungsten bronze ferroelectrics[END_REF]. Crystalizing in this particular structure, the ferroelectric tetragonal tungsten bronze Ba 2 LnFeNb 4 O 15 (TTB-Ln), with Ln being a Lanthanide, have been shown to exhibit a particular behavior, namely the generation in-situ of ferromagnetic barium hexaferrite (BaFe 12 O 19 : BaFO), during the proper sintering of the ferroelectric matrix [START_REF] Josse | The Ba2LnFeNb4O15 "tetragonal tungsten bronze": Towards RT composite multiferroics[END_REF][START_REF] Roulland | Influence of ceramic process and Eu content on the composite multiferroic properties of the Ba6-2xLn2xFe1+xNb9-xO30 TTB system[END_REF][START_REF] Castel | In-situ formation of barium ferrite in iron-doped "tetragonal tungsten bronze": Elaboration of room temperature multiferroic composites[END_REF]. This processes therefore produces a spontaneously forming multiferroic ceramic composite and overcomes any chemical compatibility issue that are common with conventional composites.

For complete understanding of the multiferroic properties and their microscopic distribution, this contribution reports a complementary investigation of the ferroelectric and ferromagnetic properties of TTB-Sm/BaFO composite ceramics synthesized by solid state reaction [START_REF] Josse | The Ba2LnFeNb4O15 "tetragonal tungsten bronze": Towards RT composite multiferroics[END_REF]. By studying the dielectric coefficients and the magnetization versus the temperature, we determined the ferroelectric and the ferromagnetic transition temperatures, which confirms the room temperature multiferroic nature of the composite. In addition, we used microelectromechanical investigations to study the polar domains at the ceramics surface and determine their piezoelectric coefficients. We also studied the magnetic response at the microscopic level, and we determined the spatial distribution of the magnetic secondary phase.

Materials and methods

Tetragonal tungsten bronze multiferroic composite ceramics were prepared by solid state reaction method.

The precursors BaCO 3 , Fe 2 O 3 , Nb 2 O 5 and Sm 2 O 3 (having a purity higher than 99%) were mixed by mechanical grinding and pressed under uni-axial pressure of 4 ton/cm² into cylindrical pellets with a diameter of 8 mm and a thickness around 1.2 mm before being sintered at 1300 °C during 4 hours in air. The ceramics crystallization and phase identification were characterized at room temperature by X-ray diffraction collected using a 4-circle X-ray diffractometer (PANalytical X-Pert PRO MRD). In order to prepare the samples for ferroelectric and dielectric measurements, we deposited, on each side of the pellet, a 200 nm thick layer of gold by using thermal evaporation under vacuum giving disk-shaped electrodes of 6 mm diameter each. Dielectric characteristics were performed with Wayne-Kerr 6425 component analyzer under dry helium at frequencies ranging from 10 3 to 10 6

Hz. The electric field dependence of the polarization was measured using an analogical Sawyer-Tower circuit with compensation of stray capacitances and resistances. The microstructure, the local electromechanical properties and the microscopic magnetic response were characterized by an atomic force microscope (AFM) using a DI-enviroscope AFM (Bruker, Santa Barbara, CA). To image and manipulate the ferroelectric polarization of the samples at the nanoscale, we used piezoresponse force microscopy (PFM) and we applied a voltage of 3 V at 20 kHz between the conductive tip and a silver past bottom electrode. The surface-induced piezoelectric vibrations were detected using a lock-in amplifier from Signal Recovery (model 7265, Wokingham, UK). Magnetic force microscopy (MFM) was performed to investigate the local magnetic response of the samples. The microscopic measurements were carried out with conductive/magnetic tips, having a Co/Cr coating (MESP from Bruker, Santa Barbara, CA), thus suitable for both electric and magnetic measurements. At the macroscopic level, the variation of the magnetization as a function of temperature and as a function of an applied magnetic field were studied using an EV9 Vibrating Sample Magnetometer (VSM) produced by Microsense (Lowell, MA, USA). Figure 1 shows the X-ray diffractogram of the composite ceramics. The diffraction peaks shown with the Miller indices are attributed to the principal phase TTB-Sm according to the reference file JCPDS 00-059-0425, confirming that the principal phase crystallizes in the tetragonal system with the tungsten bronze structure belonging to with the space group P4/mbm. The coherent domain size is determined to be around 35 nm using Scherrer's formula [START_REF] Scherrer | Bestimmung der Größe und der inneren Struktur von Kolloidteilchen mittels Röntgenstrahlen[END_REF], which is in accordance with [START_REF] Hajlaoui | Highly oriented multiferroic Ba 2 NdFeNb 4 O 15 -based composite thin films with tetragonal tungsten bronze structure on silicon substrates[END_REF]. A sketch of the TTB-Sm crystal structure seen along the crystallographic c-axis is drawn using the software VESTA and the crystallographic data provided in [START_REF] Josse | The Ba2LnFeNb4O15 "tetragonal tungsten bronze": Towards RT composite multiferroics[END_REF][START_REF] Bovtun | Relaxor-like behavior of lead-free Sr2LaTi2Nb 3O15 ceramics with tetragonal tungsten bronze structure[END_REF], is shown in the inset of Figure 1. In this structure, the barium ions occupy the pentagonal sites, the samarium ions are contained in the square sites, while the triangular sites are kept empty. The Fe 3+ and Nb 5+ metal ions are statistically distributed within oxygen octahedral, which are responsible to the appearance of the ferroelectric polarization [START_REF] Josse | The Ba2LnFeNb4O15 "tetragonal tungsten bronze": Towards RT composite multiferroics[END_REF][START_REF] Castel | In-situ formation of barium ferrite in iron-doped "tetragonal tungsten bronze": Elaboration of room temperature multiferroic composites[END_REF]. Moreover, SmNbO 4 phase (#) is always detected, due to the limited structural distortions of the TTB phase, which results only in a partial accommodation of the samarium ions [START_REF] Josse | The Ba2LnFeNb4O15 "tetragonal tungsten bronze": Towards RT composite multiferroics[END_REF][START_REF] Castel | In-situ formation of barium ferrite in iron-doped "tetragonal tungsten bronze": Elaboration of room temperature multiferroic composites[END_REF]. The remaining ions interact with a small amount of niobium and oxygen to form the Fergusonite phase. In order to maintain electrical neutrality, iron ions can only be partially incorporated in the TTB framework and the excess portion of iron oxide reacts with a small amount of barium oxide in order to form the BaFO phase (•) [START_REF] Josse | The Ba2LnFeNb4O15 "tetragonal tungsten bronze": Towards RT composite multiferroics[END_REF][START_REF] Roulland | Influence of ceramic process and Eu content on the composite multiferroic properties of the Ba6-2xLn2xFe1+xNb9-xO30 TTB system[END_REF][START_REF] Josse | Des fluorures aux ferroïques, l'empire de la cristallochimie[END_REF].

Results and discussion

Although the fergusonite phase does not present ferroic properties at room temperature, its presence is required to trigger the formation of the magnetic phase of hexaferrite. In addition, the formation of the magnetic BaFO phase is very desirable because of the large use of this phase in several applications, such magnetic recording and microwave applications [START_REF] Chen | Barium hexa-ferrite technology for MAMIR and advanced magnetic recording applications[END_REF][START_REF] Sharrock | Application of barium ferrite particles in advanced recording media[END_REF][START_REF] Ozah | Nanosized barium hexaferrite in novolac phenolic resin as microwave absorber for X-band application[END_REF]. We found that Curie temperature determined here by dielectric measurements is slightly higher than that determined in [START_REF] Josse | The Ba2LnFeNb4O15 "tetragonal tungsten bronze": Towards RT composite multiferroics[END_REF] (T C ≈ 405 K). In the paraelectric region (T>T C ), the real part of the permittivity does not obey the linear behavior expected from the Curie-Weiss law in the paraelectric state, as a change of slope is obvious around 430 K. This can be related to the behavior observed in the Ba 2 NbFe(Nb 1-x Ta x ) 4 O 15 solid solution, and to the existence of a metastable ferroelectric phase that appear as a general feature in this TTB series [START_REF] Josse | Original crystal-chemical behaviors in (Ba,Sr)2Ln(Fe,Nb,Ta)5O15 tetragonal tungsten bronze: Anion-driven properties evidenced by cationic substitutions[END_REF][START_REF] Albino | Metastable ferroelectric phase and crossover in the Ba2NdFeNb4-xTaxO15 TTB solid solution[END_REF].

The temperature dependence of the imaginary part of the dielectric permittivity ( ) at various frequencies has also been measured and is depicted in Figure 2(d). A peak having a maximum position independent to the frequency at T = T C is obtained, which is clearly related to the paraelectric-ferroelectric phase transition.

Moreover, a second dielectric anomaly is observed at lower temperature of around 250 K, suggesting a second phase transition at reduced temperature that can be related to the recurrent ferroelectric to relaxor crossovers observed in this series [START_REF] Josse | The Ba2LnFeNb4O15 "tetragonal tungsten bronze": Towards RT composite multiferroics[END_REF][START_REF] Josse | Original crystal-chemical behaviors in (Ba,Sr)2Ln(Fe,Nb,Ta)5O15 tetragonal tungsten bronze: Anion-driven properties evidenced by cationic substitutions[END_REF][START_REF] Castel | Flexible relaxor materials: Ba2PrxNd1-xFeNb4O15 tetragonal tungsten bronze solid solution[END_REF]. The presence of a very wide peak which is extending over a large frequency range indeed indicates that the low temperature transition is a relaxor-like ferroelectric transition [START_REF] Peláiz-Barranco | Relaxor Behaviour in Ferroelectric Ceramics Chapter[END_REF][START_REF] Cross | Ferroelectrics Relaxorferroelectrics : An overview[END_REF].

Finally, the imaginary part of the permittivity increases drastically at high temperatures (T> 432 K), which reflects the presence of extrinsic contributions related to the establishment of a significant electrical conductivity [START_REF] Zhang | Oxygen-vacancy-related relaxation and conduction behavior in (Pb1-xBax)(Zr0.95Ti0.05)O3 ceramics[END_REF]. This dielectric investigation, however, confirms that the studied ceramics are ferroelectric at room temperature. In order to get more information about the local ferroelectric and piezoelectric properties of the studied system at room temperature, we evaluated the variation of the longitudinal (d ZZ ) and the transversal (d XZ ) piezoelectric coefficients [START_REF] Harnagea | Local ferroelectric switching properties in BiFeO3 microstructures and their piezomagnetic response[END_REF][START_REF] Harnagea | Towards ferroelectric and multiferroic nanostructures and their characterisation[END_REF] by cycling an external voltage via the PFM tip fixed above selected locations on the sample surface. The dependence of d ZZ and d XZ on the applied bias pulse is described by a hysteresis loop as shown in Figure 3(c) and 3(d), respectively. The well-saturated character of the hysteresis curves proves that the polarization can be reversibly switched between two opposite stable states, which clearly confirms the ferroelectric nature of the material. This result suggests that this material could possibly be used for memory devices that require the polarization to switch between two states with an angle of 180° [START_REF] Palneedi | Status and Perspectives of Multiferroic Magnetoelectric Composite Materials and Applications[END_REF]. The maximum value of piezoelectric coefficients d ZZ and d XZ are estimated at 1.3 pm/V and 1 pm/V, respectively. The difference of coercive voltages in both positive and negative directions is related to the fact that the top electrode (the AFM/PFM tip) has both a very different geometry and a different composition than the continuous bottom electrode [START_REF] Hajlaoui | Enhanced ferroelectric properties in multiferroic epitaxial Ba2EuFeNb4O15 thin films grown by pulsed laser deposition[END_REF]. The average coercive voltage is similar in both cases of piezoelectric coefficients, which is estimated at ≈ 4.5 V. In order to determine the nature and the spatial distribution of the magnetic phase embedded in the ferroelectric matrix, the magnetic properties have been studied at macroscopic and microscopic levels, as shown in Figure 4.

First, the hysteresis loop describing the variation of the macroscopic magnetization as a function of an applied magnetic field was studied at room temperature as depicted in that these magnetic domains are probably formed by several grains of BaFO. This proves that magnetic domains are indeed spontaneously formed within the ferroelectric matrix during its thermal treatment, and demonstrates that their magnetic properties can be experimentally observed, giving rise to room temperature multiferroic composite in the form of ceramics.

Conclusion

Multiferroic composite ceramics TTB-Sm/BaFO were synthesized and their multifunctional properties were investigated. The ferroelectric and the magnetic phase transition temperatures were determined to be 417 K and 590 K, respectively. In addition, ferroelectric and ferromagnetic hysteresis loops were measured at room temperature. The microelectromechanical characterization confirmed the presence of ferroelectric properties down to nanoscale and allowed the determination of the piezoelectric coefficients. The formation of the magnetic phase in the form of microdomains was confirmed by studying the magnetic response of the ceramic at the microscopic level. In this study, not only the multiferroic properties of the TTB-Sm/BaFO composite ceramics have been investigated, but the distribution of the polar and magnetic properties was locally studied as well. Finally, the coexistence of the polar and magnetic domains in the same material suggests the investigation of possible couplings between these domains (magnetoelectric coupling) as a prospective for this study.

Fig. 1 :

 1 Fig. 1: X-ray diffractogram of the synthesized TTB-Sm/BaFO ceramic composite: in addition to the diffraction peaks of the TTB-Sm (shown with their Miller indices), the peaks of the secondary phases BaFe 12 O 19 and SmNbO 4 are shown with the symbols (•) and (#), respectively. The inset displays a unit cell of TTB-Sm, indicating the different channels with pentagonal (filled with Ba 2+ ), square (filled with Sm 2+ ) and triangular (empty) shapes and the positions of the different ions.

Fig. 2 :Figure 2 (

 22 Fig. 2: Room temperature ferroelectric properties of the synthesized ceramics: (a) and (b) represent the variation of the polarization and switching current as a function of an applied voltage, respectively, for various electric field amplitudes ranging from 12kV/cm to 26 kV/cm. The inset of (a) shows the variation of the coercive field as a function of the maximum electric field applied. Variation of the dielectric permittivity as a function of temperature and frequency: (c) real part and (d) imaginary part of the dielectric permittivity.

Fig. 3 :

 3 Fig. 3: Microscopic electromechanical properties of TTB-Sm ceramics investigated at room temperature: (a) and (b) show respectively the out-of-plane and the in-plane piezoelectric responses detected at the ceramics surface.

  (c) and (d) represent the applied voltage dependence of the out-of-plane and the in-plane piezoelectric coefficients, respectively. In order to investigate the microelectromechanical properties of the studied ceramics and get space-resolved information about their ferroelectric properties at room temperature, piezoelectric force microscopy was used. The out-of-plane and the in-plane piezoelectric responses (mixed signal, i.e. amplitude cos(phase)) are shown in Figure 3(a) and 3(b), respectively. First, the bright and dark contrasts in Figure 3(a) are respectively attributed to polar regions where the out-of-plane component of the polarization is oriented upward ( ) and downward ( ) perpendicular to the plane of the image, i.e. the surface of the ceramic sample. Second, the bright and dark contrasts observed in Figure 3(b) indicate the presence of polar regions where the polarization is lying in the plane of the sample surface and oriented in opposite directions, as depicted by the arrows (↑↓). These contrasts evidence the out-of-plane and the in-plane piezoelectric vibrations, thus the piezoelectric properties of the studied material.

Fig. 4 :

 4 Fig. 4: Investigation, at the macroscopic and microscopic scales, of the magnetic properties of TTB-Sm/BaFO ceramic composite. (a) Represents the room temperature variation of the magnetization as a function of an applied magnetic field, while (b) depicts the variation of the magnetization as a function of temperature. (c) and (b) represent respectively the topography and the microscopic magnetic response of the sample surface measured at room temperature, showing the magnetic domains distribution at the ceramic surface.

Figure 4 (

 4 a), evidencing of the presence of a room temperature ferromagnetic ordering in the studied TTB-Sm/BaFO ceramic composite. The magnetic coercive field (H C ) measured to be around 2200 Oe (0.22 Tesla) corresponds to the magnetic coercive field of the in-situ forming 'hard' magnetic phase within the ferroelectric matrix TTB-Sm[START_REF] Josse | The Ba2LnFeNb4O15 "tetragonal tungsten bronze": Towards RT composite multiferroics[END_REF][START_REF] Castel | In-situ formation of barium ferrite in iron-doped "tetragonal tungsten bronze": Elaboration of room temperature multiferroic composites[END_REF][START_REF] Josse | Des fluorures aux ferroïques, l'empire de la cristallochimie[END_REF]. The magnetization as a function of temperature is studied and is depicted in Figure4(b), which shows a decrease of magnetization and a transition from a ferromagnetic to a paramagnetic phase when increasing the temperature. For accurate determination of the transition temperature, the derivation of the magnetization is plotted as a function of the temperature (Figure4(b)), and its minimum is obtained at 590 K, which is ascribed to the magnetic transition temperature of the TTB-Sm/BaFO samples. The high transition temperature together with the large magnetic coercive field are in concordance with the magnetic signature of BaFO phase. This result, together with the ferroelectric behavior discussed above, prove that these ceramics are multiferroic at room temperature.To determine the spatial distribution of the magnetic phase at the sample surface, magnetic force microscopy was used. The surface topography and the magnetic response measured at room temperature are shown in Figure4(c) and 4(d), respectively. The image of the topography shows that the microstructure of the ceramic is formed by a grain structure of lateral dimensions varying from 0.8 μm to 3 μm, giving rise to a roughness of about RMS ≈ 250 nm. The presence of domains with dark contrast in the image of the magnetic response reveals the presence of magnetic domains within a non-magnetic matrix (the ferroelectric matrix TTB-Sm), at locations indicated by the red circles in Figure4(c) and 4(d). A careful comparison between these two images indicates the absence of any correlation between the topography and the magnetic response, confirming that the contrasts observed in Figure4(c) are solely due to the presence of magnetic domains at or near the sample surface. The average size of the magnetic domains (which we associate with the BaFO phase) is around 1.38 μm, indicating

Funding: this research was funded by the Natural Sciences and Engineering Research Council of Canada (NSERC) Individual Discovery Grant # RGPIN 261662. Conflicts of interest/Competing interests: no conflict related to this work.


Laure Guerit

Xiao-Ping Yuan

Sebastien Carretier

Stephane Bonnet

Sébastien Rohais

Jean Braun

Delphine Rouby

Sébastien Carretier

Stéphane Bonnet

Fluvial landscape evolution controlled by the sediment deposition coefficient: Estimation from experimental and natural landscapes

de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

INTRODUCTION

Continental landscapes constantly evolve in response to tectonics and climate via erosion, sediment transport and deposition. Based on a simple mass balance, the evolution of a topography through time can be expressed as:

where is the topographic elevation (m), is the time (yr), is the uplift rate (m/yr), is the erosion rate (m/yr), and is the deposition rate (m/yr).

According to the stream power model [START_REF] Howard | Channel changes in badlands[END_REF][START_REF] Howard | Modeling fluvial erosion on regional to continental scales[END_REF][START_REF] Whipple | Dynamics of the stream-power river incision model: Implications for height limits of mountain ranges, landscape response timescales, and research needs[END_REF][START_REF] Lague | The stream power river incision model: evidence, theory and beyond: Earth Surface Processes and Landforms[END_REF], the erosion rate at steady state along river channels can be expressed as:

where is the erodibility (m 1-3m yr m-1 ), is the effective precipitation rate (m/yr), is the drainage area (m 2 ), is the slope, and and are two coefficients. [START_REF] Kooi | Escarpment evolution on high-elevation rifted margins: Insights derived from a surface processes model that combines diffusion, advection, and reaction[END_REF] suggested that the deposition term is equal to the sediment flux (m 2 /yr) divided by a transport length (m). The transport length can be seen

G G = 0 G > 0.4 dh dt = U -E + D, (1) h t U E D E = K f (PA) m S n , ( 2
)
K f P A S m n D q s L
manuscript submitted to Geology as the distance a grain can travel within the flow before it is deposited on the river bed, and it exerts a strong control on river morpho-dynamics [START_REF] Davy | Fluvial erosion/transport equation of landscape evolution models revisited[END_REF][START_REF] Bradley | Measuring gravel transport and dispersion in a mountain river using passive radio tracers: Earth Surface Processes and Landforms[END_REF][START_REF] Ganti | Quantitative bounds on morphodynamics and implications for reading the sedimentary record[END_REF]Kaspark et al., 2015;[START_REF] An | Morphodynamic model of the lower Yellow River: flux or entrainment form for sediment mass conservation?[END_REF]. [START_REF] Davy | Fluvial erosion/transport equation of landscape evolution models revisited[END_REF] improved this formulation by relating to the sediment concentration in the water column and to the net settling velocity of the grains which depends on the flow and particles characteristics. The transport length has been implemented as proportional to the water discharge (m 2 /yr) in several models to give:

where is a scaling factor (yr/m). Equations (1-3) are used to model small range and foreland up to mega-fans [START_REF] Carretier | Modelling sediment clasts transport during landscape evolution: Earth Surface Dynamics[END_REF][START_REF] Mouchené | Autogenic versus allogenic controls on the evolution of a coupled fluvial megafan-mountainous catchment system: numerical modelling and comparison with the Lannemezan megafan system (northern Pyrenees, France)[END_REF][START_REF] Shobe | The SPACE 1.0 model: a Landlab component for 2-D calculation of sediment transport, bedrock erosion, and landscape evolution[END_REF][START_REF] Langston | Developing and exploring a theory for the lateral erosion of bedrock channels for use in landscape evolution models[END_REF].

For a given water discharge, when is low, it is the ability of river to transport sediments that controls the rate of landscape evolution (transport-limited system, TL), whereas it is controlled by the ability of rivers to erode and provide sediments to the system when is high (detachment-limited system, DL) [START_REF] Davy | Fluvial erosion/transport equation of landscape evolution models revisited[END_REF]. These two endmember modes lead to similar steady-state landscapes but their transient behavior can be quite distinct [START_REF] Tucker | Topographic outcomes predicted by stream erosion models: Sensitivity analysis and intermodel comparison[END_REF][START_REF] Davy | Fluvial erosion/transport equation of landscape evolution models revisited[END_REF]. For example, the response of a landscape to a climatic change will be faster for a TL system than for a DL one [START_REF] Armitage | Numerical modelling of landscape and sediment flux response to precipitation rate change[END_REF]. Accurate estimation of is thus required for landscape evolution models and to understand the response of continental landscapes to perturbations in tectonics or climate.

However, the scaling factor is defined from physical parameters related to the particles in transport and to characteristics of the flow that are difficult to estimate in natural systems.

L q

w D = q s L = q s q w ζ , (3) ζ ζ ζ ζ ζ manuscript submitted to Geology

APPROACH

We develop here a method allowing a simple estimation of the deposition coefficient from the landscape morphology. Considering local equilibrium, the sediment flux can be expressed as the integrated upstream net erosion rate. Assuming uniform precipitation over the upstream catchment of area , the water flux is equal to the product of and We thus rewrite the deposition term as:

where is a dimensionless coefficient which controls the local deposition rate.

This coefficient has a strong impact on the drainage density and landscape morphology [START_REF] Carretier | Colluvial deposits as a possible weathering reservoir in uplifting mountains: Earth Surface Dynamics[END_REF]their Ndepo), and [START_REF] Davy | Fluvial erosion/transport equation of landscape evolution models revisited[END_REF] proposed that (their equation 16) determines the DL-or TL-like nature of the topographic evolution ( or for ) . Equation (1) thus reads [START_REF] Yuan | A new efficient method to solve the stream power law model taking into account sediment deposition[END_REF]:

Consider an uplifting area surrounded by a sedimentary basin where alluvial fans can develop.

We define as the vertical displacement relative to the basin's basement (i.e., the sedimentary basin is fixed). At steady state, the rate of topographic change is equal to zero. Accordingly, in the uplifting area, Equation (5) reduces to:

and along the fans in the sedimentary basins, where there is no uplift, it reduces to:

q s P A q w A P . D = 1 APζ ∫ A ( U - ∂h ∂t ) d A = G A ∫ A ( U - ∂h ∂t ) d A, (4) G = 1/ ( Pζ ) G G > 1 G < 1 n = 1 ∂h ∂t = U -K f (PA) m S n + G A ∫ A ( U - ∂h ∂t ) d A . ( 5
) U ∂h / ∂t (1 + G)U = K f (AP) m S R n , ( 6
)
where SR and SF refer to the river slopes in the range and along the fan, respectively. At the coupling point of these two domains, combining Equations ( 6) and ( 7) leads to:

We can thus estimate the deposition coefficient from the slopes of the rivers near the coupling point between a catchment and its fan. This equation is derived at the coupling point of catchment-fan system and therefore, only local steady-state river profile is required to estimate from equation ( 8). Diffusive processes along hillslopes could be added to Equation ( 5). However, we focus here on river profiles at the catchment outlets and we thus neglect this process for the sake of simplicity.

In the following, we apply this methodology to experimental and natural catchmentfan system to estimate and the corresponding scaling factors . We then discuss implications in terms of sediment transport and catchment dynamics.

FROM EXPERIMENTAL LANDSCAPES

We selected two sets of laboratory experiments where an uplifted area submitted to rainfall erosion is coupled to a depositional zone [START_REF] Babault | Influence of piedmont sedimentation on erosion dynamics of an uplifting landscape: an experimental approach[END_REF][START_REF] Rohais | Architecture stratigraphique et flux sédimentaires sur la marge sud du golfe de Corinthe (Grèce) : Analyse de terrain, modélisations expérimentales et numériques[END_REF][START_REF] Rohais | Sedimentary record of tectonic and climatic erosional perturbations in an experimental coupled catchment-fan system[END_REF]. Details about the experiments are presented in the GSA Data Repository Text 1 and are summarized here. [START_REF] Rohais | Architecture stratigraphique et flux sédimentaires sur la marge sud du golfe de Corinthe (Grèce) : Analyse de terrain, modélisations expérimentales et numériques[END_REF] and [START_REF] Rohais | Sedimentary record of tectonic and climatic erosional perturbations in an experimental coupled catchment-fan system[END_REF] used a rectangular erosion box with a single outlet coupled to a depositional plateau to simulate an uplifting catchment and its alluvial fan (Fig. 1A, Fig. DR1). The catchment only is submitted to rainfall. [START_REF] Babault | Influence of piedmont sedimentation on erosion dynamics of an uplifting landscape: an experimental approach[END_REF] used a

GU = K f (AP) m S F n , (7) G = S F n S R n -S F n . (8) G G G ζ G manuscript submitted to Geology
rectangular box surrounded by a depositional plateau to simulate an uplifted range with multiple outlets and associated fans (Fig. 1B, Fig. DR2). The whole landscape is submitted to rainfall. In the 14 experiments analyzed, the erosion box is filled with similar silica powder, uplift is uniform and applied at a constant rate in the catchment areas, and erosion is induced by rainfall from sprinklers located above. Optical stereo data were acquired to build DEMs of the experiments (horizontal resolution: 0.5-1 mm/pixel, vertical resolution <0.5 mm/pixel).

Based on a large dataset of experiments [START_REF] Moussirou | Modulation of the erosion rate of an uplifting landscape by long-term climate change: An experimental investigation[END_REF], [START_REF] Moussirou | Modélisation expérimentale de la dynamique long terme des reliefs en réponse aux forçages climatiques et tectoniques[END_REF] shows that the time to reach steady-state non-linearly depends on . According to the stream power model, this implies an exponant in Equation ( 2), with a value of for the material used here. This is in agreement with estimations from similar experiments [START_REF] Viaplana-Muzas | Modeling of drainage dynamics influence on sediment routing system in a foldand-thrust belt[END_REF].

With each experiment, we used landscapes at steady-state (constant mean elevation through time in the uplifted zone) and determined the drainage network with Topotoolbox Matlab subroutines (Schwanghart and Khun, 2010). We extracted the river profiles along a few centimeters on both sides of the river outlets and calculated the river slope in the range and along the fan by linear regressions (Figs. 1C-D, DR1, DR2, Table DR1). In total, 29 catchment-fan systems were analyzed. The experiments are performed with the same lithology and erosion processes. We therefore expect to be similar in all the experiments.

Estimates of are consistent whatever the experimental designs or the applied forcings: values range between 0.4 and 1.9 (Figs. 2 and DR3, Table DR1). A whisker plot suggests that is an anomaly (Fig. 2B). Excluding this point, we observe that the experiments are characterized by with a median value of 0.7 (Fig. 2B).

U n > 1

n = 1.7 ± 0.2 S R S F G G G = 1.9 G ∈ [ 0.4 -1.2 ]
manuscript submitted to Geology

FROM NATURAL LANDSCAPES

We applied the same approach to natural landscapes and selected 68 systems of coupled catchment-fan systems worldwide based on the visual identification of alluvial bodies in sedimentary basins. We excluded fans that are cut or deformed by obvious tectonic features to avoid out-of-equilibrium systems. Two well-individualized catchment-fan systems and six systems of adjacent fans were selected as an equivalent to the two sets of experiments. We also analyzed fan deltas to document extensional settings. The data set cover a variety of tectonic and climatic settings and is composed of small (debris-flow dominated) to very large (mega-fans) systems (Table DR2, Figs. DR4-13).

For each site, we mapped the sedimentary system on GoogleEarth images (Figs. 3A and DR4-13) and we used Topotoolbox Matlab subroutines (Schwanghart and Khun, 2010) to extract the main river path close to the apex using ASTER DEMs (horizontal resolution: 30 m, vertical resolution: 10 m, Fig. 3B). A river profile at steady state can be converted into a straight line by a slope transformation ( plots introduced by Perron and Royden, 2012) .

Such plot are less influenced by noise associated to DEM resolution than simple slope plots.

We therefore used plots to identify the longest linear segment of the river upstream of the apex (typically a few km) and we assumed that it is at steady state. We then measured the slopes by linear fit along this segment and along the fan, a few km downstream of the apex (Figs. 3C, DR4-13, Table DR2). Several studies suggest for natural landscapes and a median of 2 emerges from worldwide database (DiBiase and [START_REF] Dibiase | The influence of erosion thresholds and runoff variability on the relationships among topography, climate, and erosion rate[END_REF][START_REF] Whittaker | Tectonic and climatic controls on knickpoint retreat rates and landscape response times[END_REF][START_REF] Lague | The stream power river incision model: evidence, theory and beyond: Earth Surface Processes and Landforms[END_REF][START_REF] Harel | Global analysis of the stream power law parameters based on worldwide 10Be denudation rates[END_REF]. In the following, we thus use for all sites.

Estimates of are more scattered than the experimental ones. Three systems (5% of the data set) show a slope continuity from the catchment to the fan, implying .

G

χ -h χ -h n > 1 n = 2 ± 0.5 G G → ∞ manuscript submitted to Geology
Excluding these systems, ranges from ∼0 to 3.1 (Fig. 4 and DR3, Table DR2). A more linear relationship between slope and erosion would lead to higher values of but in the same order of magnitude (Fig. DR14).

DISCUSSION

To calibrate , we assume that and are uniform in space at the scale of the catchment-fan system. These assumptions are correct in the experiments made with a unique material evolving under uniform rainfall. For these systems, we suggest that the variability observed around the median value (~, Fig. 2) is related to local variations in material properties (i.e. particle size) and in precipitation rate, which can not be perfectly controlled in such setups. This could be considered as an intrinsic variability associated to .

The variability observed for natural systems is higher than the experimental one (Fig. 4). We propose that it partly comes from violation of our assumptions (e.g., steady state, uniform and ). In natural settings, and may vary at the scale of the catchment-fan system. However, most channels are covered by sediments (e.g., [START_REF] Turowski | Hydraulic geometry, river sediment and the definition of bedrock channels[END_REF], and within-catchment variations should be averaged at the river outlets where we estimate . At the scale of a catchment-fan system, it is thus reasonable to assume locally uniform and that is the catchment-averaged precipitation rate. Within a given study site, different lithologies between catchments could induce differences in and in grain sizes which potentially influence sediment transport and deposition.

We also assume landscapes at steady-state, which is again correct for the experiments.

For the natural systems, the plots of the river profiles close to the apexes suggest that these segments are in equilibrium. However, alluvial fans are typically sensitive to temporal

G G G K f P ± 50% G K f P K f P G K f P K f χ -h manuscript submitted to Geology
variations at the scale of 10 2 -10 5 yrs [START_REF] Hinderer | From gullies to mountain belts: A review of sediment budgets at various scales[END_REF] and systems responding to such variations (out-of-equilibrium systems, ) could explain part of the variability observed on .

Beyond those uncertainties, we observe a correlation between and the slope of the river in the range, and between and the drainage area (Fig. DR15A-B). Steep rivers and small catchments tend to be associated with low values of and vice versa. This may reveal a difference in fluvial dynamics. In fact, based on previous studies, we propose that experimental landscapes tend toward a TL mode. Using a similar laboratory setting, [START_REF] Lague | Laboratory experiments simulating the geomorphic response to tectonic uplift[END_REF] showed that their experiments are best reproduced by a TL model. Based on similar experiments, [START_REF] Turowski | Experimental channel response to tectonic uplift[END_REF] concluded that a DL model is unable to reproduce their observations. We thus argue that the coefficient characterizes the erosional and transport behaviors of a landscape for a large range of spatial and temporal scales. Based on our estimates, we propose that indicates a system tending toward a TL mode.

About 1/3 of the natural systems studied here show a lower than 0.4 (Fig. 4). Those systems should thus be close to a pure DL mode, well-described by the stream power model.

The remaining 2/3 are characterized by a large range of (from 0.4 to ) which should indicate systems transitioning from DL to TL mode, up to pure TL mode ( ). Natural landscapes thus seem to describe a continuum between the two modes with a preference for TL mode (as suggested by [START_REF] Davy | Fluvial erosion/transport equation of landscape evolution models revisited[END_REF], supporting the use of landscape models coupling erosion and deposition in a single equation.

Finally, we can estimate the scaling factor from (Equation ( 4)). In the experiments, we observe similar estimated with a median of 1.7×10 -3 yr/m (Fig. DR16A).

For the natural systems, the estimated show a large dispersion with a median of 17 yr/m

∂h / ∂t ≠ 0 G G G G G G > 0.4 G G ∞ G → ∞ ζ G ζ ζ
manuscript submitted to Geology (Fig. DR16B). [START_REF] Davy | Fluvial erosion/transport equation of landscape evolution models revisited[END_REF] defined the scaling factor as , where

is the ratio of the sediment concentration near the bed and within the flow, and is the effective settling velocity of the particles. The estimated values of for the natural systems lead to unrealistic low settling velocities of maximum 10 -8 m/s. When measured from natural topography, the apparent scaling factor thus appears to be overestimated. This could be related to an underestimation of the morphogenic (as we only have access to the present-day, average ) or to the intermittency of natural rainfall that induces intermittent transport. Such transport can lead to an overestimation of (GSA Data Repository Text).

CONCLUSIONS

We propose a simple formulation for fluvial landscape evolution taking into account sediment deposition. The deposition term is characterized by a dimensionless coefficient .

This coefficient is related to a scaling factor that controls river bedforms and fluvial dynamics from a detachment-limited (DL) mode to a transport-limited (TL) mode. Although this factor is key to describe fluvial dynamics and landscape evolution, is difficult to measure on natural landscapes. On the contrary, is easy to estimate from topographic data.

We measured this coefficient on experimental and natural landscapes across a large variety of

Figure 1 .

 1 FIGURE CAPTIONS

Figure 2 .

 2 Figure 2. A) Estimations of the deposition coefficient from the 14 experiments (gray dot s:

Figure 3 .

 3 Figure 3. Extraction of the average river slope in the range and in the fan for natural

Figure 4 .

 4 Figure 4. Estimations of the deposition coefficient from several natural systems

  1 GSA Data Repository item 201Xxxx, with supplementary method, tables and figures, is available online at www.geosociety.org/pubs/ft20XX.htm, or on request from editing@geosociety.org or Documents Secretary, GSA, P.O. Box 9140, Boulder, CO 80301, USA.

	∞
	G
	n
	G ∈ [ 0.4 -1.2 ]

ACKNOWLEDGMENTS

This work is part of the COLORS project funded by Total. We thank John Armitage, Katy Barnhart, Jeffrey Kwang, an anonymous reviewer and the journal editor for their comments that greatly improved the manuscript.


Shannon M Mackenzie

Jason W Barnes

Jason D Hofgartner

Samuel P D Birch

Matthew M Hedman

Antoine Lucas

Sébastien Rodriguez

Elizabeth P Turtle

Christophe Sotin

Sebastien Rodriguez

December

The case for seasonal surface changes at Titan's lake district

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Titan, Saturn's largest moon, hosts lakes and seas of liquid hydrocarbons at its poles 1 . General circulation models demonstrate that regional evaporation and precipitation rates of methane likely change with the seasons (Titan's year ∼ 29.5 Earth years) and evolve on a geological timescale (∼ 10 5 Earth years) 2,3,4 . Cassini observations suggest shoreline recession at a few south polar lakes during local summer 3 , but similar seasonal changes have yet to be observed at the north pole where lakes are larger and more numerous 6,3,7 . We present three "phantom lakes" that appear to be north polar surface liquids in winter observations by Cassini RADAR but which are inconsistent with lakes in infrared images obtained up to 7 years later, after vernal equinox, suggesting that the liquids were removed in between. If this were the case, the phantom lakes could be interpreted as shallow ponds, with either a pure methane composition or a regolith porous enough to remove the less volatile ethane. These phantom lakes provide observational constraints on removal timescales for surface liquids at Titan's north pole. The location, size, and longevity of surface liquid reservoirs affect sediment processing 7 , seasonal weather 8 , climate evolution 9 , and even, perhaps, their habitability 10 . As solubility of the possible non-polar mixtures is generally low, short-lived lakes might be nutrient-poor 10 and thus have low astrobiological potential.

Complete emptying of small or shallow lakes may be the easiest shoreline change to identify in Cassini data, given that the largest predicted evaporation and infiltration rates peak at meters per year 11,6,12 . For example, though observations of Ontario Lacus, Titan's largest south polar lake, were interpreted as evidence of kilometers of shoreline recession over four years 8,3 , later analysis 14 contested that such change was within the uncertainty of spatial sampling. The difficulty in identifying even the most dramatic liquid level change lies in establishing a long enough baseline between observations with sufficient spatial resolution to discern changes. Data from a single instrument are rarely able to produce suitable comparisons for this purpose. But, by combining data from RADAR, the Visual and Infrared Mapping Spectrometer (VIMS) and the Imaging Science Subystem (ISS) we can take advantage of 13 years of Titan-targeted flybys at the cost of added complexity.

Each instrument's dataset is unique in spatial and temporal resolution and coverage, making precise comparison between them nontrivial. Furthermore, the instruments probe different properties of the surface. VIMS and ISS are sensitive to composition and grain size on the scale of a few microns while Synthetic Aperture RADAR (SAR) images convey information about dielectric properties and surface roughness on the centimeter scale.

Lakes are identifiable due to their lacustrine morphology and characteristic darkness in each dataset. Light specularly reflects off the surface of the smooth lakes 15,21 , often away from the direction of Cassini. In the infrared (IR), diffuse scattering off haze in the atmosphere also provides some illumination, making lakes appear dark grey to ISS (0.93 µm) and dark blue to VIMS (0.9 µm -5 µm) in the color scheme of Figures 1 and2. At 2.17 cm, the observed RADAR signal is unaffected by the atmosphere 17 and is a combination of quasispecular backscatter and diffuse backscatter 15 . For lakes, SAR observes low signal return resulting from little to no quasispecular contribution (due to the smooth lake surface); weak diffuse contribution from volume scattering within the liquid if it is nearly homogenous liquid; and reduced scattering off the lake floor boundary (due to the small Fresnel reflection coefficient between liquid and solid organics) 18 . Attenuation within a liquid column begins to affect the signal when path lengths are greater than tens of meters 19,20 (Supplemental Figure 1). Each instrument is sensitive to liquid depths of at least a few wavelengths.

The location and extent of Titan's lakes and seas are thus some of the strongest cases of agreement between RADAR, VIMS, and ISS datasets. Three surface features, however, defy that agreement. These "phantom lakes" appear liquid-filled and lacustrine in shape to RADAR but not to the IR instruments in later observations. Observed by RADAR in SAR mode on July 22 and October 9 of 2006 1 (top panel, Figure 2), the phantom lakes are dark in both the denoised 2 images (Figure 1) and raw images (Figure 2). Their normalized radar cross sections (σ 0 ) are more consistent with nearby lakes than dry surfaces. In Figure 3 we show the σ 0 of the phantom lakes relative to empirical fits of empty lakes (red) and partially-filled lakes (blue); shading represents 95% confidence intervals 3 . A second observation on May 23, 2013 (Figure 2 inset) shows that phantom lake B remained dark enough to agree with the lake fit, though is brighter than predicted.

Of course, RADAR darkness does not a lake make. We do not definitively rule out that some combination of smoothness or reduced volumetric scattering 22 makes these features RADAR-dark, dry surfaces. However, the lacustrine morphology of each candidate makes the presence of liquids a more likely explanation, especially given the phantom lakes' locations within the lake district and far from, for example, the dark dunes of the equatorial region. While phantom lake A lacks the sharp edges seen at some lakes 23 , it is locally confined at a scale similar to neighboring lakes and has diffuse boundaries like those seen at Jingpo and MacKay Lacus; A is therefore unlikely to be a dark plain (examples of which are plotted as grey circles in Figure 3).

Phantom lake B is a dark area along the southern shoreline of Oneida Lacus, a lake confirmed to have liquids by echoes in the April 2017 altimetry pass (which unfortunately ends before reaching phantom lake B). Overlap between SAR swaths T19 and T91 indicates that the rim extends partially along the southern shore, separating half of the phantom lake terrain from the liquid of Oneida. Thus, phantom lake B's morphology suggests a plain in contact with a lake shoreline, potentially a flood plain. Finally, C's boundaries resemble those of neighboring sharp-edge depressions 7 .

Because the ISS and VIMS cameras could not observe Titan's north pole until the Sun rose with northern spring (August 2009) and the winter clouds dissipated 24 . ISS first glimpsed the lake district on October 13, 2013. In that image, phantom lake A was not dark like a lake while phantom lakes B and C were not resolved. A month later (Figure 2 middle), phantom lakes A and C appear bright, inconsistent with the presence of surface liquids. The observation of phantom lake B is less conclusive due to its proximity to Oneida Lacus.

VIMS followed up on December 1, 2013 (Figure 2 bottom). At as fine as 4 km/pixel, at least four spectral units-surface materials with similar reflectance properties in the IR-are evident. Evaporites 25,26 appear bright pink/orange in our infrared color images because of their diagnostic brightness at 5 µm; they surround some of the dark blue lakes. (The blue and green materials are spectrally distinct but not yet geomorphologically interpreted.) VIMS pixels covering phantom lakes A and B do not behave as a single unit as they do in RADAR.

In Figure 4 we show that the southern pixels of A (black square) are darker than the northern pixels (grey square). Similarly, the eastern half of B (grey star) is darker than the western (black star). These bright western pixels of B parallel the raised rim of Oneida Lacus. Though the outline of the phantom lakes is unmistakeable in RADAR, A and C do not show these boundaries in VIMS. Phantom lake B's two VIMS units do correspond to the RADAR boundaries, providing morphological evidence that this RADAR-dark patch is a unique and isolated feature. Phantom lake C's spectrum most closely resembles the evaporites of Ypoa and Suwa Lacus, nearby lakes that appear liquid filled in all datasets.

Subsequent VIMS observations are too coarse for in-depth analysis but appear consistent with a lack of surface liquids. ISS images from February 2014 and July 2017 show that the phantom lakes are bright even in the final Cassini glimpses of the lake district.

Thus, the phantom lakes scatter like nearby solid surfaces in the infrared. They cannot be open bodies of liquid at the time of VIMS or ISS observations. How then can we reconcile the liquid radar-darkness with an IR bright surface? We put forth three possible explanations for the phantom lake observations: (1) constant liquid level with a RADAR-invisible topmost layer, (2) an IR scattering layer developing between observations, and (3) liquid removal from at least the top few microns of the surface. In all cases, the radar-darkness arises from a liquids on or saturating the surface at the time of the RADAR observations.

In scenario (1), no surface change occurs. Rather, the difference between the RADAR and IR observations stems from the difference in wavelengths observed by the RADAR mapper and the IR imagers. For example, the diffuse scattering off suspended Rayleigh particles much smaller than a RADAR wavelength with dielectric constants near that of liquid methane could contribute a small enough σ 0 that the lake would appear radar-dark 18 . In this case, the scattering layer must be less than a RADAR wavelength (2.17 cm) to be effectively transparent to the instrument.

In scenario (2), the scattering layer could develop in the time between the RADAR and IR observations. In this case, the IR bright layer need only be thick enough to obscure the the liquid below.

The scattering "scum" layer could be bubbles, haze fallout, or suspended sediment atop a liquid layer. But, as even a few mm of liquid would show significantly decreased reflectance to VIMS 27 , the overall IR brightness of the phantom lakes requires that the scattering layer be optically thick (a few mm). Furthermore, as no ethane absorption features are present in any of the phantom lakes' spectra, ethane cannot be present in the scattering layer (Supplemental Figure 3). This, along with the fact that the phantom lakes do not exhibit the same boundaries in VIMS as they do in RADAR, suggests that bubbles, haze fallout, and suspended sediment are unlikely explanations for the IR brightness of the phantom lakes.

Alternatively, a dry rind sitting atop a saturated regolith would not necessarily show unique spectral boundaries and is thus the more likely no-change scenario. Huygens observations 28 suggest even Titan's equatorial desert has liquid methane close to the surface. Modeling Titan's climate with extensive surface liquid reservoirs better reproduces the observed distribution of seas and cloud activity 9,12,4,[START_REF] Turtle | Titan's meteorology over the Cassini mission: Evidence for extensive subsurface methane reservoirs[END_REF] . The phantom lakes could therefore be highly localized areas of near surface wetted regolith distinct from the surrounding features, an explanation perhaps more likely for phantom lake B given its morphological similarity to a flood plain of Oneida Lacus and less likely for A and C given their lacustrine morphologies and distance from other lakes. If the phantom lakes are the first Cassini identification of regolith saturated up to the top few mm, it is curious that such saturated regolith is limited to a few isolated locales given that lakes are widely distributed throughout the polar region.

We therefore currently favor (3) in which liquids are responsible for the dark σ 0 and are removed to at least a few mm below the surface by the time VIMS and ISS observe the phantom lakes. A freshly exposed solid surface could explain both the IR brightness and the similarity between the phantom lakes' VIMS units and those of the surrounding terrain. Processes for emptying a lake on Titan include draining by channels, infiltration into a porous regolith, and evaporation. Cassini could not resolve any channels feeding into the phantom lakes; if any such features exist they must be subpixel in size (< 300 m).

For features A and C, the emptying process would have to take place within the seven years between observations, a timescale that easily fits within the estimates for emptying via a combination of evaporation and infiltration. If the phantom lakes are shallow, evaporation would empty them more quickly than a nearby deep lake-creating more easily resolvable shoreline change-and should pick up with increased insolation during northern summer, in line with the sequence of our observations. Infiltration can happen on relatively short timescales: even the least permeable regoliths hypothesized on Titan can drain lakes of 10 m to a methane table a few centimeters below the surface within several years 3 . Modeled evaporation rates for Titan range from ∼ 0.01 mm/Earth year 12 to meters/Earth year 11 . The most conservative evaporation estimates (0.01 mm/Earth year) constrain the evaporated volumes of phantom lakes A and C to as small as 10 6 m 3 ; corresponding to depths too shallow to be sufficiently dark to RADAR. Thus, infiltration of about the same volume would have to have played a role.

Feature B, on the other hand, would have a much shorter timescale for change. Though brief, the six months between observations may be enough time to remove liquid from the shallow depths of the surface that the IR instruments sense. Phantom lake B appears brighter to RADAR at higher incidence angles than would be expected for a shallow lake. Scattering off suspended particles can increase σ 0 , but so can looking through less liquid, i.e. a decrease in the amount of liquid at phantom lake B. Furthermore, B is situated along the shoreline of Oneida Lacus, enhancing the likelihood of a shallow deposition of liquid (e.g. lake overflow) or draining into a nearby reservoir (e.g. deeper into the subsurface). As a raised rim separates the western half of phantom lake B from Oneida Lacus, it is more likely for liquid to overflow on the eastern half, consistent with the difference in VIMS units.

Together, the phantom lakes suggest that observable change may have occurred at Titan's north pole during the Cassini mission at timescales (6 months -7 years) comparable to that observed for south polar lakes 3 of similar size. This seasonal timescale would also be consistent with the range of removal rates predicted for evaporation 11,12 and infiltration 6 when operating in conjunction. Barring substantial local (less than a few km) differences in evaporation and/or infiltration properties, we expect that larger and deeper lakes are exchanging their liquids between the surface and atmosphere at similar rates, but their greater depths make their shoreline charges undetectable by Cassini. Because our temporal coverage relies on three datasets that probe different surface properties, we consider both liquid level change and physical phenomena as possible explanations.

Figure 3: Normalized radar cross sections (σ 0 ) from SAR swaths as a function of incidence angle for the phantom lakes (A,B,C); lakes (blue circles); empirical fits for partially-filled lakes (blue) and empty lakes (red) with 95% confidence intervals shaded 3 ; and non-lacustrine, uniform, SAR-dark plains 7 (dark grey circles). O and U represent the σ 0 of Oneida Lacus and an unnamed adjacent lake, respectively; each have altimetry returns consistent with a specular reflection off a smooth, presumably liquid surface. σ 0 is a measure of the microwave reflectivity of a material. The standard errors for each feature are smaller than the plotted symbols. The flyby number for each of the two B and U observations are indicated in subscripts.

Figure 4: Behavior of phantom lakes relative to other VIMS spectral units within the atmospheric windows. Because VIMS can only see down to Titan's surface in seven windows (0.93, 1.0, 1.3, 1.6, 2.0 , 2.7, 2.8, and 5 µm), we group surface materials into spectral units defined by the relative reflectance at these wavelengths (b,c). In the T96 VIMS image (CM 1764550739 16; a), we identify four units (lake in purple, evaporite in pink, green, and blue) as endmembers to which we compare phantom lake pixels (examples labeled in top image). Phantom lakes A and B have two groups of pixels (black, and grey), suggesting two spectral units within the area that RADAR sees as a dark unit. Atmospheric effects have been mitigated with a radiative transfer model.

Methods

After calibration (described below), we imported images into QGIS, an open-source GIS, for projection and mapping.

In Synthetic Aperture Radar (SAR) mode, the Cassini RADAR instrument mapped the surface at resolutions of ∼ 350 m/pixel at 2.17 cm 1 . We use incidence angle, noise-subtracted, and NLDSAR 2 swaths from T16, T19, and T91. Incidence and σ 0 values were extracted via polygon mapping and statistical tools within QGIS. The fits shown in Figure 3 are taken directly from a previously published empirical fit to SAR data 3 .

ISS imaged Titan at 0.938 µm 4 . At this wavelength, haze particles in the atmosphere scatter most of the light before it reaches the surface. However, the shape of the haze phase scattering function 5 actually facilitates observation of the surface because some of that light forward scatters to the ISS detector. The resulting ISS maps are generally 1 km per pixel. Thus, ISS data have been used to identify the presence of liquid on Titan's surface 6,7,8 as well as monitor cloud activity in the atmosphere 6,9 . For this work, we used the global mosaic map available from the USGS, as well as data from October 13, 2013 to February 2, 2014 available on the PDS. The latter were calibrated and mosaicked using USGS Integrated Software for Imagers and Spectrometers (ISIS).

VIMS observed Titan in the near infrared windows at which methane absorption (and that of trace gases) is weak enough to allow light to penetrate through the atmosphere and interact with the surface 10,11 : 0.93 µm, 1.0 µm, 1.3 µm, 1.6 µm, 2.0 µm, 2.7µm, 2.8 µm, and 5 µm. Though specific chemical composition is difficult to identify 12,13 , the diverse surface can be divided by spectral reflectance properties into spectral units 14 .

VIMS data were calibrated and colored 14 for initial identification of spectral units, shown in Figures 1 and2. We identified four spectral units in the high resolution cube-lake, evaporite, blue, and green-based on the unique behavior of the band ratios. These endmember representative spectra were created by averaging over pixels within the identified units. To rule out atmospheric phenomena, we applied an empirical correction 15 and also modeled the spectra with the radiative transfer model 16,17,18 . In both data sets, the spectral behavior of the phantom lakes shown in Figure 4 is consistent and agrees with the raw data. The radiative transfer model, however, allowed us to significantly reduce the standard deviation of the spectral ratios with respect to the raw and empirically corrected data, thus reinforcing the conclusions that the phantom lakes are not liquid filled in the VIMS data.

The radiative transfer model used is that of Hirtzig et al. ( 2013) 16 with a look-up-table scheme 17 as described in detail by Brossier et al. (2018) 18 . Though the model relies on Huygens ground truth from the equatorial region, the decreased standard deviation of the spectral ratios of all units indicates that the radiative transfer model fairly mitigated the atmospheric contributions, one of the major sources of data variability.

There are many IR bright phenomena in VIMS observations of Titan. Specular reflections off smooth surfaces 19,20,21 are bright at all wavelengths (especially 5 µm) but are only observed when Cassini images the surface at the same angle at which the light is incident. None of our phantom lakes were imaged at this geometry. Broad specular reflections from wavy surfaces can occur further from the specular point but are also much brighter than the reflectance of the phantom lakes 22 . Furthermore, no specular reflections have been identified in ISS data to date, as the the high optical depth at 0.93 µm makes direct surface scattering unlikely even at favorable geometries. Thus, neither direct nor broad specular reflections can explain the brightness of the phantom lake features to ISS. Clouds are also bright but have a distinct, ephemeral spectral signature 23,24,25,26 that is inconsistent with that of the phantom lakes (Supplemental Figure 6). As the phantom lakes do not appear to drastically change between observations (October 2013-July 2017) with VIMS and ISS, clouds are unlikely to hide liquids. The adjacency effect can artificially brighten dark surfaces when photons reflected from nearby bright surfaces get scattered into the cone of acceptance of the pixel imaging the dark material. However, this effect dominates at shorter distances with increasing optical depth 27 and we observe the opposite: the phantom lakes appear brighter at the longer wavelengths. The brightness of the phantom lakes is therefore most likely due to compositional and grain size properties of the surface and not some observational or atmospheric phenomenon.

Figure 1 :

 1 Figure 1: Cassini images of Titan's north pole. (a) Azimuthal stereographic projection of Titan's north pole as observed by Cassini VIMS on December 13, 2013. Titan's hydrocarbon seas (maria) lie opposite the lake district and at lower elevation 30,7 . (b-e) Eastern edge of the lake district as seen by VIMS (up to 4 km/pixel; b), ISS (map scale at 4 km/pixel; c), and Non-Localized Denoised 2 SAR (350 m/pixel; d). These instruments observed the eastern half of the lake district at different wavelengths and seasons (northern winter and northern spring). A,B, and C labels indicate the phantom lake features on the NLDSAR image. Relevant lake names are given on the ISS map.

Figure 2 :

 2 Figure 2: Sequence of observations of the phantom lakes, outlined in each dataset. SAR coverage from T16 and T19 (a) and T91 (b) show the phantom lakes as dark. But later, the same areas

Acknowledgements

This work was supported under the NASA Earth and Space Science Fellowship Program grant NNX14AO30H to SMM. JWB acknowledges support from NASA Cassini Data Analysis Program NNX15AI77G. AL and SR acknowledges the financial support of the UnivEarthS Labex program at Sorbonne Paris Cité (ANR-10-LABX-0023 and ANR-11-IDEX-0005-02). SR is also supported by the French National Research Agency (ANR-APOSTIC-11-BS56-002, ANR-12-BS05-001-3/EXO-DUNES). Part of this work was carried out at the Jet Propulsion Laboratory, California Institute of Technology, under contract with NASA. Copyright 2018. Government sponsorship is acknowledged.

Data Availability

All data presented here are available from the NASA Planetary Data System (https://pds-imaging.jpl.nasa.gov/portal/cassini mission.html), except the NLDSAR swaths which are maintained by coauthor Lucas (http://cssnldsar.geophysx.org/). The data that support the plots within this paper are also available from the corresponding author upon reasonable request.

Author Contributions SMM led the procuring and analysis of the data and wrote the manuscript. JWB assisted in all aspects of the analysis and writing. JH and SPDB contributed to the analysis of RADAR data. MMH helped develop scattering models. AL wrote the code to produce NLDSAR and provided these data. All authors contributed to the discussion of these results. SR provided the radiative transfer model and contributed to their analysis. EPT assisted in the calibration of ISS data and planned the ISS observations. CS planned the VIMS Titan observations.

Competing Interest

The authors declare no competing interest. *Correspondence or request for materials should be directed to SMM (shannon.mackenzie@jhuapl.edu).