Dr Camille Litty
email: camille.litty@univ-grenoble-alpes.fr

Julien Charreau

Pierre-Henri Blard

Raphael Pik

Sébastien Nomade

Spatial variability of Quaternary denudation rates across a volcanic ocean island (Santo Antão, Cape Verde) from cosmogenic 3 He

Denudation of volcanic ocean islands creates remarkable landscapes and contributes to Earth's carbon cycle, since the chemical alteration of basalts is a CO2 sink. Because many volcanic islands have large climate gradients and relatively low variations in lithology and tectonic history, they represent excellent natural laboratories for studying climatic effects on landscape evolution. However, little is known about the control of denudation rates in ocean islands and the respective influences of climatic gradients and morphological parameters. Here, we present new measurements of long-term denudation rates from Santo Antão Island, Cape Verde (17°N). In this 779 km 2 island, mean annual precipitation ranges from 40 mm.yr -1 in the southwest to 1100 mm.yr -1 in the north. To constrain the spatial distribution of denudation rates, we measured the cosmogenic 3 He concentration in river-transported pyroxene grains from 23 river bedload samples. We obtained basin wide denudation rates ranging from 2.7 ± 0.1 to 57.5 ± 0.3 m/Ma. The denudation rates display a significant spatial variability, with the highest rates in catchments located in the northeast side of the island where modern precipitation are the highest and low denudation rates in the southern and western dry basins. Our study shows that precipitation is the main control on denudation and landscape development of the Santo Antão volcanic island. This study provides for the first time the spatial distribution of denudation rates across a volcanic island located in a tropical zone.

Introduction

Climate is thought to be one of the main parameters controlling denudation rates (e.g., [START_REF] Portenga | Understanding Earth's eroding surface with 10 Be[END_REF], which in turn impact many processes such as the building and evolution of mountain ranges (e.g., [START_REF] Willett | Orogeny and orography: The effects of erosion on the structure of mountain belts[END_REF][START_REF] Whipple | The influence of climate on the tectonic evolution of mountain belts[END_REF]. Understanding the effects of climate on denudation rates is therefore thought to be critical, but the modality of its impact remain a matter of long-standing debate (e.g., [START_REF] Langbein | Yield of sediment in relation to mean annual precipitation[END_REF][START_REF] Dibiase | The influence of erosion thresholds and runoff variability on the relationships among topography, climate, and erosion rate[END_REF]. Conversely, denudation may also exert a strong feedback on climate at the geological time scale (e.g., [START_REF] Berner | The carbonate-silicate geochemical cycle and its effect on atmospheric carbon dioxide over the last 100 million years[END_REF]. Silicate weathering on continents consumes CO2 and releases cations (Ca, Mg, K, Na) to the ocean. In the ocean, Ca and Mg are then combined with CO3 2-and precipitated into bicarbonate, resulting

Revised manuscript with no changes marked in a net CO2 sink (e.g., [START_REF] Berner | The carbonate-silicate geochemical cycle and its effect on atmospheric carbon dioxide over the last 100 million years[END_REF]. The weathering of basaltic lithologies are thus supposed to be an efficient CO2 sink at the geological time scale [START_REF] Dessert | Basalt weathering laws and the impact of basalt weathering on the global carbon cycle[END_REF].

In such context, volcanic islands in tropical areas are noteworthy natural laboratories allowing the investigation of links and interactions between denudation and climatic processes for two main reasons. First, the amplitude of the denudation variations can be compared to climate only. Indeed, most of the volcanic islands often exhibit a strong climatic gradient, while tectonic and lithological forcing are comparatively negligible [START_REF] Hildenbrand | Geomorphological study of long-term erosion on a tropical volcanic ocean island: Tahiti-Nui (French Polynesia)[END_REF][START_REF] Salvany | Geomorphic evolution of the Piton des Neiges volcano (Réunion Island, Indian Ocean): competition between volcanic construction and erosion since 1.4 Ma[END_REF][START_REF] Jefferson | Controls on the hydrological and topographic evolution of shield volcanoes and volcanic ocean islands. The Galapagos: A natural laboratory for the[END_REF].

Second, according to geochemical measurements of dissolved load in rivers, the modern weathering of volcanic islands, especially in tropical areas, is especially intense (e.g., Louvat andAllègre, 1997, Jefferson et al., 2014). It may represent up to 35% of the total CO2 consumed from continental silicate weathering [START_REF] Gaillardet | Global silicate weathering and CO2 consumption rates deduced from the chemistry of large rivers[END_REF] and, hence, have a large influence on long-term climate changes [START_REF] Lenton | Enhanced carbonate and silicate weathering accelerates recovery from fossil fuel CO2 perturbations[END_REF].

In this general context, it is of the upmost importance to better constrain denudation rates (i.e., the sum of physical and chemical erosion) and their spatial distribution across these islands. Several studies have estimated the long-term (10 5 -10 6 a) denudation rates of several islands by reconstructing the original volcanic surfaces and then by calculating the volumes removed by erosion [START_REF] Hildenbrand | Geomorphological study of long-term erosion on a tropical volcanic ocean island: Tahiti-Nui (French Polynesia)[END_REF][START_REF] Salvany | Geomorphic evolution of the Piton des Neiges volcano (Réunion Island, Indian Ocean): competition between volcanic construction and erosion since 1.4 Ma[END_REF][START_REF] Menéndez | Fluvial dissection, isostatic uplift, and geomorphological evolution of volcanic islands (Gran Canaria, Canary Islands, Spain)[END_REF][START_REF] Ferrier | Covariation of climate and long-term erosion rates across a steep rainfall gradient on the Hawaiian island of Kaua 'i[END_REF][START_REF] Ricci | Construction and destruction rates of volcanoes within tropical environment: Examples from the Basse-Terre Island (Guadeloupe, Lesser Antilles)[END_REF]. This geomorphological approach relies on the postulate that the construction of the volcanoes precedes their destruction by erosion. However, these two processes can be indeed coeval, which can hamper an accurate estimate of the initial pre-erosion surface. Moreover, it is often hazardous to accurately estimate the age of the original surfaces of the volcanoes. The derived denudation rates are therefore tied to large uncertainties and must be rigorously considered as minimum values [START_REF] Hildenbrand | Geomorphological study of long-term erosion on a tropical volcanic ocean island: Tahiti-Nui (French Polynesia)[END_REF]. Other studies have constrained modern denudation rates from direct geochemical measurements of dissolved load in rivers [START_REF] Louvat | Present denudation rates on the island of Reunion determined by river geochemistry: basalt weathering and mass budget between chemical and mechanical erosions[END_REF][START_REF] Ferrier | Covariation of climate and long-term erosion rates across a steep rainfall gradient on the Hawaiian island of Kaua 'i[END_REF]. However, these estimates remain only strictly valid over short time scales (10 1-10 2 a). Therefore, the mechanisms by which volcanic islands erode at the Quaternary time scale (10 2 -10 4 ) remain largely unknown despite being a critical parameter for a comparison with long-term climate change.

Basin wide denudation rates at the Quaternary time scale can be quantified by measuring the concentration of cosmogenic nuclides in river sediment. In-situ cosmogenic nuclides are produced by the flux of secondary cosmic particles in the top few meters of Earth's surface (e.g., [START_REF] Granger | Spatially averaged long-term erosion rates measured from in situ-produced cosmogenic nuclides in alluvial sediment[END_REF].

However, denudation rate measurements from cosmogenic nuclides have been mainly based on quartzrich lithologies through 10 Be analyses (e.g., [START_REF] Granger | Spatially averaged long-term erosion rates measured from in situ-produced cosmogenic nuclides in alluvial sediment[END_REF]. On basaltic and andesitic islands, quartz is rare or inexistent, implying that 10 Be cannot usually be used (Blard et al., 2008). Cosmogenic 3 He measured in pyroxene and olivine in river bedload sediment represents a good alternative to 10 Be for these mafic terrains (e.g., [START_REF] Puchol | Variability of magmatic and cosmogenic 3He in Ethiopian river sands of detrital pyroxenes: Impact on denudation rate determinations[END_REF]. These minerals have a high helium retentivity (Blard and Pik, 2008) and their 3 He cosmogenic production rate is relatively well established [START_REF] Balco | A complete and easily accessible means of calculating surface exposure ages or erosion rates from 10Be and 26Al measurements[END_REF]. 3 He has been widely used for the past 30 years to date lava flows (e.g., [START_REF] Kurz | Cosmogenic helium in a terrestrial igneous rock[END_REF] or to reconstruct continental paleoclimates (e.g., [START_REF] Blard | Persistence of full glacial conditions in the central Pacific until 15,000 years ago[END_REF]. However, only a few studies have used cosmogenic 3 He to estimate denudation rates [START_REF] Gayer | Spatial variability of erosion rates inferred from the frequency distribution of cosmogenic 3He in olivines from Hawaiian river sediments[END_REF][START_REF] Ferrier | Covariation of climate and long-term erosion rates across a steep rainfall gradient on the Hawaiian island of Kaua 'i[END_REF][START_REF] Puchol | Variability of magmatic and cosmogenic 3He in Ethiopian river sands of detrital pyroxenes: Impact on denudation rate determinations[END_REF]. [START_REF] Sarda | Cosmogenic neon and helium at Réunion: measurement of erosion rate[END_REF] were the first to determine denudation rates from cosmogenic 3 He in a volcanic island, but this study was limited to a depth profile at one single location in La Reunion. To quantify Quaternary denudation rates [START_REF] Gayer | Spatial variability of erosion rates inferred from the frequency distribution of cosmogenic 3He in olivines from Hawaiian river sediments[END_REF] and [START_REF] Ferrier | Covariation of climate and long-term erosion rates across a steep rainfall gradient on the Hawaiian island of Kaua 'i[END_REF] measured the cosmogenic 3 He in river sand in Kaua'I Island (Hawaii), but both studies were limited to a single drainage basin.

In this paper, we analyze the 3 He cosmogenic concentrations in pyroxenes and olivines of 23 sand samples collected in 19 different modern riverbeds across the entire Santo Antão Island (Cape Verde). We hence determine, for the first time across a volcanic island, the near-whole spatial variability of denudation rates averaged over the last 10 3 to 10 5 ka.

Geological and climatic settings of the Cape Verde islands and Santo Antão

The Cape Verde archipelago (4033 km 2 ) includes 10 major and several smaller islands located ~600 km west of the African continent (Fig. 1). The development of the islands likely initiated 19 Ma ago with tholeiitic sills that intruded the Jurassic-Cretaceous sea-floor basement [START_REF] Hayes | Mesozoic magnetic lineations and the magnetic quiet zone off northwest Africa[END_REF][START_REF] Duncan | Geochronology of basaltic rocks recovered by DSDP Leg 41, Eastern Atlantic Ocean[END_REF]. Since 19 Ma, the volcanic activity has not migrated and has been well distributed across the islands until present [START_REF] Mitchell | On dating the magmatism of Maio, Cape Verde islands[END_REF][START_REF] Plesner | 40Ar/39Ar geochronology of Santo Antao, Cape Verde Islands[END_REF]. The active volcanic activity is now limited to Fogo Island with the last eruption in 1995.

The climate of the Cape Verde archipelago is driven by the respective position of the subtropical high pressure of the North Atlantic (Açores anticyclone) and the Intertropical Convergence Zone [START_REF] Mannaerts | Rainfall erosivity in Cape Verde[END_REF]. The wet season lasts from July to October when the ITCZ, usually centered around the equator, shifts toward 10-20°N and brings 80% of the annual precipitation, whereas during the rest of the year, precipitation is scarce [START_REF] Mannaerts | Rainfall erosivity in Cape Verde[END_REF].

Santo Antão is the northernmost and second largest island (770 km 2 ) of the archipelago (Fig. 1).

Its elevation ranges from 0 to 1979 m (Tope de Coroa), with a mean elevation of 362 m. The relief of the island forms a major orographic barrier that creates a large precipitation gradient with annual precipitation that varies by a factor of 20-30 between the wettest northeast and the driest southwest sides (Fig. 3a). Local rainfall is also strongly dependent on elevation, with some low coastal areas having low rainfall of 0 -150 mm.yr -1 whereas the mountain areas above 1000 m elevation receive ~800 -1000 mm.yr -1 (Fig. 3a; [START_REF] Mannaerts | Rainfall erosivity in Cape Verde[END_REF]. This precipitation gradient seems to have induced significant topographic heterogeneities across the island: major valleys on the north and northeast side of the island receive regular rain during the wet season and are deeply incised. In contrast, the southern part of the island receives only very irregular rain and the topography exhibits wide smooth seaward slopes preserving volcanic surfaces that are poorly incised (Fig. 3b). In addition, two deep (>1 km) and large (> 5km) circular depressions are located in Tarrafal and the Cha di Morte regions (Fig. 2 and3).

These two depressions together with the deeply entrenched valley in the north provide numerous outcrops of the lavas that built the island and enable a good reconstruction of its volcanic history (e.g., [START_REF] Plesner | 40Ar/39Ar geochronology of Santo Antao, Cape Verde Islands[END_REF][START_REF] Holm | Sampling the Cape Verde mantle plume: evolution of melt compositions on Santo Antão, Cape Verde Islands[END_REF].

On Santo Antão, the oldest lavas reported are ~7.5 Ma old altered picrites and ankaramites of the Dos Bodes formation. They can be observed in a limited area at the bottom of the Cha di Morte depression [START_REF] Plesner | 40Ar/39Ar geochronology of Santo Antao, Cape Verde Islands[END_REF][START_REF] Holm | Sampling the Cape Verde mantle plume: evolution of melt compositions on Santo Antão, Cape Verde Islands[END_REF]Fig. 2). This formation is overlain by the thick (~800 m) 3-1 Ma old Cha di Morte formation (Fig. 2; [START_REF] Plesner | 40Ar/39Ar geochronology of Santo Antao, Cape Verde Islands[END_REF]. The lava flows of this formation dip toward the coast and point toward sources of magma emission located in the central part of the island. The Dos Bodes and Cha di Morte formations are usually mapped together as "Older Volcanics" by several authors (Fig. 2; [START_REF] Holm | Sampling the Cape Verde mantle plume: evolution of melt compositions on Santo Antão, Cape Verde Islands[END_REF]. After the emplacement of the Cha di Morte formation, between ~1.4 Ma and 0.3 Ma, other lavas flows were emplaced in other areas of the island, indicating that the recent volcanic activity was less focused and more distributed (Plenser et al., 2002[START_REF] Holm | Sampling the Cape Verde mantle plume: evolution of melt compositions on Santo Antão, Cape Verde Islands[END_REF]). The Tope de Coroa volcano experienced a large flank collapse that is not yet dated and yields to the large Tarrafal circular depression. This depression was later refilled by the development of the modern Tope de Coroa stratovolcano ~0.1 Ma ago. Since ~0.3 Ma, the volcanic activity of the island is mainly characterized by explosive magmatism including numerous individual scoria cones, phreatomagmatic to strombolian and several pluri-meter-thick phonolitic to tephri-phonolotic tephra deposits referenced as the Cao Grande formation [START_REF] Plesner | 40Ar/39Ar geochronology of Santo Antao, Cape Verde Islands[END_REF]Eisele et al., 2015). During this recent explosive period, several lava flows are reported in the Tope de Coroa and in the Porto Novo region.

Material and methods

Basin average denudation rates

Assuming isotopic steady state, meaning that the cosmogenic production of 3 He is balanced by the erosional loss [START_REF] Lal | Cosmic ray labeling of erosion surfaces: in situ nuclide production rates and erosion models[END_REF], approximating (similarly to the spallogenic production) the two muogenic production rates at depth with two different exponential laws, then the basin average denudation rate (cm.yr -1 ) is given by:

𝜀 = ∑ 𝑃 𝑥,3 Λ 𝑥 𝜌3 𝐻𝑒 𝑐 𝑥=𝑠𝑝,𝜇 𝑠𝑚 ,𝜇 𝑓𝑚 (1)
where the subscripts sp, 𝜇 𝑠𝑚 and 𝜇 𝑓𝑚 correspond to the production by spallation, slow muon capture and fast muon processes, respectively. 𝑃 𝑥,3 is the basin average in situ cosmogenic production rates at the surface (at.g -1 .a -1 ). ρ and Λx are the rock density (g.cm -3 ) and the attenuation length in the rocks (g.cm - 2 ) of each production pathway x, respectively. 3 Hec (at/g) is the mean concentration of the cosmogenic 3 He in the river sediment. We use a basalt density of 3 g/cm 3 and assumed constant attenuation lengths of 160 g/cm 2 , 4320 g/cm 2 and 1500 g/cm 2 for spallation, fast muons capture and slow muons processes (e.g. Braucher et al., 2011). The basin average production rates were computed using a cell-by-cell approach and the SIG tool Basinga [START_REF] Charreau | Basinga: A cell-by-cell GIS toolbox for computing basin average scaling factors, cosmogenic production rates and denudation rates[END_REF] with a 30 m digital elevation model of the Santo Antão Island. The production rates of each of the cells of the DEM were scaled using the time dependent Lal-Stone scheme (Stone, 2000), and a reference surface production at sea level and high latitude of 122±12 at.g -1 .a -1 has been used (world-wide mean of all calibration studies computed with the CREp calculator, Martin et al., 2017). Following DiBiase (2018), we did not apply any topographic shielding correction.

Determination of cosmogenic 3 He in pyroxenes and olivines

The total 3 He concentration ( 3 Hetot) in both pyroxenes and olivine is a three-component system (Blard and Farley, 2008) that can be written as follows:

3 Hetot = 3 Hec + 3 Hemag + 3 Henucl (2) where indices c, mag and nucl refer to the concentration in matrix-sited cosmogenic 3 He, the magmatic 3 He located in fluid and melt inclusions and the nucleogenic 3 He produced by neutron capture on Li nuclei since the eruption of the lavas, respectively.

Similarly, the total 4 He budget is a two-component system consisting of magmatic 4 He inherited in fluid and melt inclusions and matrix-sited radiogenic 4 He mainly produced by the decay of 238 U, 235 U and 232 Th since the eruption: 4 Hetot= 4 Herad + 4 Hemag (3) Therefore, the cosmogenic 3 Hec can be estimated from Eq (2) if the magmatic, the nucleogenic and the total 3 He can be estimated independently with accuracy and precision. The 3 Hetot can be measured by melting the sample. The nucleogenic contribution 3 Henucl can be calculated using the mean neutron flux, the Li concentration and the eruption age of the minerals [START_REF] Andrews | Natural production of tritium in permeable rocks[END_REF]. In practice, we used here the major and trace elements and the Li-concentrations of the phenocrysts measured at the SARM in Nancy to compute this nucleogenic component (Table A.1).

Because it can be quite variable from one aliquot to another (e.g.s [START_REF] Puchol | Variability of magmatic and cosmogenic 3He in Ethiopian river sands of detrital pyroxenes: Impact on denudation rate determinations[END_REF], the magmatic 3 Hemag component is difficult to estimate. Magmatic He can be preferentially released and analyzed by in vacuo mineral crushing [START_REF] Kurz | Cosmogenic helium in a terrestrial igneous rock[END_REF]. The crushing must, however, be sufficiently long to release and purge most of the trapped magmatic He. On the other hand, this crushing should not be too intense nor too long to avoid releasing the matrix sited 3 Hec and/or the 4 Herad [START_REF] Blard | Cosmogenic 3He production rates revisited from evidences of grain size dependent release of matrix sited helium[END_REF](Blard et al., , 2008)), nor contaminate the sample with atmospheric He [START_REF] Protin | Irreversible adsorption of atmospheric helium on olivine: A lobster pot analogy[END_REF].

In the present study, we crushed all the pyroxenes/olivines and analyzed the 3 He and 4 He concentrations extracted during this extraction step ("crush" indice). We fused in a high vacuum furnace additional uncrushed pyroxene/olivine aliquots (above 1400°C, to ensure a complete He release), and measured the 3 He and 4 He concentrations release ("melt" indice). Then, we describe and apply two approaches to estimate the magmatic component and to calculate the 3 He cosmogenic concentration.

Estimation of 3 Henucl is the same in both approaches.

The first approach is derived from the initial method proposed by [START_REF] Kurz | Cosmogenic helium in a terrestrial igneous rock[END_REF], and includes the latter improvements by [START_REF] Dunai | Scaling factors for production rates of in situ produced cosmogenic nuclides: A critical reevaluation[END_REF] and Blard and Farley (2008), to account for the presence of 4 Herad. In this method, the cosmogenic 3 He concentration is calculated as follow:

3 Hec = 3 Hetot -3 Hemag -3 Henucl 3 Hec = 3 Hetot -[ 4 Hemag × ( 3 Hemag / 4 Hemag)] -3 Henucl from Eq (3): 3 Hec = 3 Hetot -[( 4 Hetot -4 Herad) × ( 3 Hemag / 4 Hemag) -3 Henucl andthen: 3 Hec = 3 He melt -[( 4 Hemelt -4 Herad) × ( 3 He/ 4 He)crush] -3 Henucl (4)

The ( 3 He/ 4 He)crush is the isotopic ratio measured during the initial crushing of the pyroxenes and olivines, which is considered as representative of the magmatic isotopic ratio.

Because pedogenesis and sediment transport induces break up of phenocrysts, we consider that the pyroxenes and olivines sampled in the riverbed loads do not bear any implanted 4 Herad from the (U-Th) richer lava (Blard and Farley, 2008). This is justified because the implanted component is only present in the 20 µm outer rim (e.g., [START_REF] Dunai | Scaling factors for production rates of in situ produced cosmogenic nuclides: A critical reevaluation[END_REF]. Consequently, we computed the 4 Herad with the U-Th-Sm concentrations measured in the sole phenocrysts (Tables 1 andA.2). To estimate the mean ages of the lava flow, we used the geological map of [START_REF] List | Cape Verde Islands: Santo Antão, São Vicente, Sal, Fogo, Thematic Maps[END_REF], and weighted each eruption age with the respective contributive size of each studied basin. Additionally, if the measured 4 He concentration during melting is lower than the calculated 4 Herad, then we assume that the 4 Hemag is null in these aliquots and we only apply a correction for 3 Henucl.

The second approach assumes that the in vacuo crushing step was efficient enough to release all the magmatic 3 He from the aliquots and that the measured concentrations in the crushed aliquots are representative of those of the melted aliquot. In this case, Eq. ( 2) can be rewritten as follows:

3 Hec = 3 He melt -3 Hecrush -3 Henucl (5) with the indices melt and crush referring to the 3 He measured after melting and crushing experiments of the sample. The calculated cosmogenic 3 He concentrations from both approaches; Eqs. ( 4) and ( 5) are then used to estimate the average denudation rates in the different sampled river catchments on Santo Antão Island applying Eq. ( 1) (e.g., [START_REF] Granger | Spatially averaged long-term erosion rates measured from in situ-produced cosmogenic nuclides in alluvial sediment[END_REF].

Sampling, treatment and analyses

Twenty-three samples were collected at the mouth of the main rivers all around the island (Fig. 3,Table A.1). Approximately 3 kg of river bedload sediment was collected from each sampling site.

Bulk sands were sieved under flowing water into the 0.15-0.5 mm granulometric fraction. A hand magnet was first used to remove a large fraction of the magnetic sediment. Then the pyroxene and olivine grains were concentrated from the sand by two magnetic separations using a Frantz device.

Nearly pure pyroxene and olivine grains were then obtained by gravimetric separation in di-iodomethane (density ~ 3.3 g cm -3 ). After washing in diethyl ether and drying, the grains were cleaned in an ultrasonic bath of deionized water, and then acetone, for 10 min each, to remove any impurities resulting from dust. After drying, all grains with attached lava fragments were removed and discarded by handpicking under a binocular microscope. We chose to work mainly on pyroxenes because they were more abundant in our samples and thus easier to separate. Nevertheless, when olivines were present in some samples, we also analyzed mixed olivine/pyroxene aliquots. To measure the total 3 He and 4 He concentrations, grain aliquots were wrapped in copper foils and fused at 1350°C for 15 min. High temperature extractions were performed with the full metal induction furnace designed at CRPG [START_REF] Zimmermann | A new single vacuum furnace design for cosmogenic 3He dating[END_REF]. After extraction, the gas was purified using activated charcoal at -196 °C and with hot (400°C and 600°C) and ambient temperature (21°C) titanium sponge getters. A complete description of the furnace and the purification line can be found in [START_REF] Zimmermann | A new single vacuum furnace design for cosmogenic 3He dating[END_REF]. Helium was analyzed on a GV HELIX Split Flight Tube dual collector mass spectrometer. The instrument sensitivity was determined using the HESJ standard [START_REF] Matsuda | The 3He/4He ratio of the new internal He Standard of Japan (HESJ)[END_REF] whose certified 3 He/ 4 He ratio is 20.63 ± 0.10 Ra (Ra = 1.384  10 -6 being the atmospheric 3 He/ 4 He ratio). To measure the Hemag concentration, invacuo crushing extractions were performed on additional aliquots with a soft iron slug activated by external solenoids (5 min at 100 strokes/min). After a first measurement, to optimize the release of all the Hemag, second crushing extractions of the same aliquot were performed. In the first approach, the mean value of the two crushing experiments was used to calculate the ( 3 He/ 4 He)crush ratios (i.e., the magmatic 3 He/ 4 He ratio). In the second approach, the sum of the two crushing experiments has been used to compute the 3 Hec as this approach assumes that the crushing extracts the whole magmatic He.

Finally, to test the reproducibility of the methods, nine samples have been replicated and analyzed after high temperature extraction. For these replicates, magmatic corrections were performed using the same He results from in-vacuo crushing extraction. 1.

Results

Information about samples and catchment metrics are given in supplementary

Mean 4 He and 3 He blanks during heating extraction were (7.3 ± 0.3)  10 -15 mol and (1.24 ± 0.2)  10 -19 mol, respectively. Repeated measures after sample heating extraction always yielded He concentrations in the blank range, indicating that the extraction was complete after the first heating step.

Mean 4 He and 3 He sample values during heating extraction before blank correction were (1.10 ± 0.03)  10 -12 mol and (4.2 ± 0.2) 10 -18 mol, respectively, one to three orders of magnitude higher than the measured blank concentrations. This implies that blanks do not represent a significant source of uncertainty. Mean 4 He and 3 He blanks measured from the crushing device (with the slug moving) were (8.6 ± 0.4)  10 -15 mol and (1.4 ± 0.1)  10 -19 mol, respectively. Mean 4 He and 3 He sample values for the crushing experiments before blank correction were (5.8 ± 1.5)  10 -14 mol and (6.5 ± 1.7)  10 -19 mol, respectively. Blanks therefore represent 10 to 20% of the crushing concentrations. The calculated concentrations have been then corrected with the mean blank values. This corresponds to total 4 He and 3 He concentrations measured after heating experiments ranging between (13 ± 1)  10 10 and (27 ± 1)  10 11 at/g and between (22 ± 3)  10 5 and (29 ± 1.3)  10 6 at/g, respectively (Table 1). The magmatic 4 He and 3 He concentrations measured after the first in-vacuo crushing extractions range between (6.3 ± 0.1)  10 9 and (2.6 ± 0.2)  10 11 at/g and between (1.9 ± 0.1)  10 4 and (8.8 ± 0.4)  10 6 at/g, respectively (Table A.1). Corresponding 3 He/ 4 He ratios measured by crushing range between 2 ± 14 Ra and 15 ± 10 with a weighted mean of 12.2 ± 0.3 Ra (Ra = 1.384×10 -6 being the atmospheric 3 He/ 4 He ratio). Second in-vacuo crushing extractions show similar results and, when summed with the first crushing, returned total magmatic values that range between (6.3 ± 0.1)  10 9 and (5.8 ± 0.2)  10 11 at/g and between (1.9 ± 0.1) ± 10 4 and (1.1 ± 0.06)  10 7 at/g for 4 He and 3 He, respectively. The isotopic ratios of the two crushing together range between 2 ± 14 Ra and 14.1 ± 0.9 Ra, with a weighted mean of 11.8±0.4 Ra. However, some samples have not been measured for a second in-vacuo crushing extraction.

With the first approach, calculated cosmogenic 3 He concentrations range between (2.1 ± 1.4)  10 6 at/g and (2.0 ± 0.17)  10 7 at/g (Table 1). Eight aliquots show 3 Hemag > 3 Hemelt therefore yielding null cosmogenic 3 Hec concentrations. The second approach leads to 3 Hec concentrations ranging between (2.0 ± 0.3)  10 6 and (2.3 ± 0.1) 10 7 at/g. Importantly, the two approaches yield 3 Hec concentrations and denudation rates that are close to the 1:1 line for most of the samples (Fig. 4a andb), indicating that the choice of the correction procedure for the non-cosmogenic 3 He components is not a significant source of uncertainty. Based on these cosmogenic 3 He concentrations, we determined the denudation rates of 23 samples in 19 main catchments all across the Santo Antão Island. Considering the two approaches, and ignoring the few samples that yielded null-cosmogenic 3 He concentrations, denudation rates show a variability across the island, from 2.7 ± 0.1 m/Ma (CV18-01) to 58 ± 4 m/Ma (CV18-20), with a mean denudation rate over the studied catchments of 18.0 ± 5.2 m/Ma. Finally, the nine replicated measurements are compatible within analytical uncertainties (Fig. 5a,b,c). For denudation rates lower than 15m/Ma, the correlation between replicates is robust and is not dependent on the approach (Fig. 5d; R 2 = 0.76, p-value = 0.002).

Discussion

Reliability of the calculated cosmogenic 3 He concentrations and derived denudation rates

The reconstruction of the 3 Hec concentrations relies on several assumptions and corrections that must be considered. First, the two approaches give similar results, which suggest that the concentration in cosmogenic 3 He has been accurately quantified. However, when looking in detail, the 3 Hec concentrations obtained with the second approach are lower by ~10 to 20% compared to those of the first approach, for most of the samples (Fig. 4a), leading to ~10 to 20% higher denudation rates (Fig. 4b). The first approach depends on the accuracy of the 4 He radiogenic and the magmatic ( 3 He/ 4 He) estimates. The assumption that the implanted radiogenic 4 He from the (U-Th) richer lava can be neglected is supported by the fact that, if it is considered, computed radiogenic 4 He are larger than the total 4 He measured in almost 90% of the analyzed aliquots. On the other hand, the lower 3 Hec concentrations obtained with the first approach could suggest an underestimate of the 4 Herad. This could result from an underestimate of the eruption age of the bedrock material. Rivers may actually erode more efficiently the older lithologies than the younger ones. Alternatively, older rocks are porphyritic and probably contain more pyroxenes than the younger ones. However, we carefully estimated the mean bedrock ages in each basin, propagating high uncertainties considering the age distributions within a given basin (from 10 to 93%, average = 60%; Table A.1). For that reason, the first approach yields higher uncertainties (mean 1σ error of 20%) than the second one (mean 1σ error of 10%). Second, the lower cosmogenic 3 He concentrations calculated with the first approach compared to the second may also suggest an overestimation of ( 3 He/ 4 He)mag ratio from the crushing experiments. Too intense crushing could have partially released the 3 Hec trapped in the crystal lattice [START_REF] Blard | Cosmogenic 3He production rates revisited from evidences of grain size dependent release of matrix sited helium[END_REF](Blard et al., , 2008)).

Moreover, the good reproducibility between the nine replicates suggests that merging hundreds of grains in each aliquot seems sufficient to average the 3 Hec concentration, despite possible spatial variability of the denudation in the watershed [START_REF] Puchol | Variability of magmatic and cosmogenic 3He in Ethiopian river sands of detrital pyroxenes: Impact on denudation rate determinations[END_REF]. Samples yielding denudation rates above ~15 m/Ma display an increase of the inter-aliquot scatter, with up to 57% of differences between the replicates (Fig. 5c). This observation suggests that, in Santo Antão Island, ~15 m/Ma seems to be a threshold value, above which denudation rates can only be measured with lower accuracy and precision. This threshold is directly controlled by the cosmogenic/magmatic ratio, as well as by the amount of radiogenic 4 He. It may thus vary from one geological setting to another and it may represent a significant limitation, particularly in catchments where denudation rates are high (>> 10 m/Ma; [START_REF] Puchol | Variability of magmatic and cosmogenic 3He in Ethiopian river sands of detrital pyroxenes: Impact on denudation rate determinations[END_REF]. [START_REF] Foeken | Long-term low latitude cosmogenic 3He production rate determined from a 126 ka basalt from Fogo, Cape Verdes[END_REF] have calibrated the local 3 He cosmogenic production rate in the nearby Fogo Island to 75.6±7.0 at.g -1 .yr -1 in olivine and 70.1±2.9 at.g -1 .yr -1 in pyroxene. These values are significantly lower than other 3 He production rates found elsewhere in the world and the reference surface production at sea level and high latitude of 122±12 at.g -1 .yr -1 (Martin et al., 2017) that we use in our calculations. [START_REF] Foeken | Long-term low latitude cosmogenic 3He production rate determined from a 126 ka basalt from Fogo, Cape Verdes[END_REF] explained this discrepancy by geomagnetic field strength spatial variability. Considering these local production rates would decrease all our Cape Verde denudation rate values but would not change their relative spatial variations.

In summary, we conclude that we accurately measured and estimated the different He components and that the two approaches provide reliable estimates of the cosmogenic 3 Hec concentration. Because the uncertainties associated to the second approach are lower than those of the first one, and because both approaches produce results that are identical within uncertainties, we restrain the following discussion to the denudation rates derived from the second approach.

Denudation rates variability in Santo Antão

The studied hydrological basins of Santo Antão are characterized by a broad range of basinaveraged denudation rates, ranging between ~3 and ~60 m/Ma and therefore varying by a factor of 20 across the island (Fig. 6; Table 1). The basins in the south and west sides of the island display the lower calculated denudation rates ranging from ~3 and ~20 m/Ma. One sample (CV18-26) located in the south of the island but shows a relatively high denudation rate (26.1 ± 5.4 m/Ma and 36.6 ± 4.4 m/Ma with the first and second approach respectively) and compared to the other southern basins seems to be an outlier that is not representative of this subregion. The northern basins are characterized by a gradient of the denudation rates from the southeast to the northwest. Finally, the region of the Cha di Morte displays the highest denudation rates of the island, with values reaching ~60 m/Ma (Fig. 6). These denudation rates inferred from 3 Hec are averaged over characteristic timescales (Tint = (Λ//ɛ). In Santo Antão, the integrated timescale of these denudation rates ranges from ~10 ka for the fastest eroded basins in the north, to ~200 ka for the southern basins, with a mean value of ~50 ka for the island (Table 1).

However, in the south side of the island, the calculated integration time is in some places higher than the age of the rocks. This suggests that the southern side could exhibit even lower denudation rates than the calculated rates. This variability of long-term denudation rates variability across Santo Antão seem to correlate with the modern precipitation pattern on the island. Basins located in the northeast and receiving relatively high precipitation are characterized by the highest denudation rates. Basins located in the south receiving the lowest annual precipitation rates of the island yield the lowest denudation rates (Fig. 3a and6). This observation is consistent with the study of [START_REF] Ferrier | Covariation of climate and long-term erosion rates across a steep rainfall gradient on the Hawaiian island of Kaua 'i[END_REF]: by reconstructing the eroded volume in 24 basins across Kaua'I in Hawaii, these authors found that long-term (Ma) denudation rates on the island are correlated with modern annual precipitation.

However, given the typical integration time of our cosmogenic 3 He-derived denudation rates of 53 ka, spatial patterns in precipitation rates may have been different in the past, meaning that the denudation rates we have measured could result from superimposed successive climatic settings. In Santo Antão, the trade winds have been stable since at least the last glacial maximum, 20 ka [START_REF] Rognon | Holocene Saharan dust deposition on the Cape Verde Islands: sedimentological and Nd-Sr isotopic evidence[END_REF].

Another study even suggests that trade winds have been present over the last 140 ka [START_REF] Hooghiemstra | Variations of the NW African trade wind regime during the last 140 000 years: changes in pollen flux evidenced by marine sediment records[END_REF]. Such timescales are similar to, or longer than, the integration time of the cosmogenic 3 He-derived denudation rates. This created a stable mean gradient of precipitation that caused one side of the island to erode at higher rates than the other. Moreover, during this time span the volcanic activity that could have reshaped and reconstructed the surface of the islands was relatively limited [START_REF] Plesner | 40Ar/39Ar geochronology of Santo Antao, Cape Verde Islands[END_REF].

Therefore, climatic conditions can explain the spatial variations of the long-term denudation rates across the island.

On the other hand, the topographic patterns also strongly vary across the island and is also in apparent correlation with the variability of denudation rates. Steeply sloping valleys are located in the northeast side of the island, in the Tarrafal and Cha di Morte depressions where denudation rates are high (> 40 m/Ma). The southern side of the island presents relatively low denudation rates, a relatively smooth morphology and low slope angles (Fig. 7). [START_REF] Gayer | Spatial variability of erosion rates inferred from the frequency distribution of cosmogenic 3He in olivines from Hawaiian river sediments[END_REF] also explained variations in cosmogenic 3 He measured in 26 samples across the Waimea basin in Kaua'I (Hawaii) by a non-linear link between denudation rates and hillslope gradient. Yet the topography of the island could be inherited from its volcanic construction [START_REF] Lahitte | DEM-based reconstruction of southern Basse-Terre volcanoes (Guadeloupe archipelago, FWI): Contribution to the Lesser Antilles Arc construction rates and magma production[END_REF] independently from the precipitation pattern.

To better explore and separate the correlation existing between denudation rates, precipitation rates and slopes, we performed a principal component analysis (PCA) using XLSTAT software [START_REF] Jolliffe | Principal Component Analysis[END_REF].

The PCA for the mean annual precipitation, mean and maximum slopes in the catchments resulted in a distance biplot of the first two principal axes shown in Fig. 8. The first and second axis explained 74.6% and 18.6% of the variance, respectively. Samples were grouped in five clusters according to their denudation rates (Fig. 8). The striking feature of this PCA is that basins with low denudation rates (< 10 m/Ma) present both low mean precipitation and low mean and maximum slopes at the same time. Basins with denudation rates between 10 and 40 m/Ma display both relatively high mean precipitation and relatively high slopes in the catchments. The sample CV18-26, which seemed to be an outlier from a visual point of view, is also an outlier here as it does not belong to the right part of the graph (Fig. 8).

Past studies have suggested strong coupling and feedbacks between climate, denudation and landscape evolution in active mountain ranges (e.g., Willet, 1999;[START_REF] Whipple | The influence of climate on the tectonic evolution of mountain belts[END_REF]. This study reinforces this idea and moreover suggests that these coupling and feedback processes also exist on volcanic islands, and that these forcing are not necessarily linear. In Santo Antão, the precipitation patterns over time control the potential for physical erosion and weathering [START_REF] Jefferson | Controls on the hydrological and topographic evolution of shield volcanoes and volcanic ocean islands. The Galapagos: A natural laboratory for the[END_REF]. The spatial variability of denudation rates influences the differential evolution of the morphology over the island, initiating a probable positive feedback: high precipitation led to higher denudation rates on the northeast side of the island, allowing for the creation of deeply incised basins with high slope angles, high drainage density and fluvial valleys that penetrate inland to the central part of the island, amplifying the local denudation.

In contrast, low precipitation leading to relatively low denudation rates on the south side of the island causes the formation of weakly dissected basins with low slopes, little tributary development and lack of deep valleys. Over the long term, the morphology of the island may thus have enhanced the denudation variability through positive feedbacks with the slopes where high slopes contributing to increase the potential for fluvial erosion. Therefore, topography on the island seems to be controlled by denudation rates rather than the contrary, at least at the initial stage of the island history. Moreover, climatic conditions seem to be the predominant process that has shaped the morphology of Santo Antão.

However, the highest denudation rates of the island (> 40 m/Ma) are correlated with high precipitation rates, but are not necessarily associated with the highest slopes, suggesting that other parameters, such as geological properties, may also play a role. Similarly, denudation rates estimated for CV18-19 and CV18-20, which are two of the highest values of the island, but located in the lee side of the island, can also result from other causes. (i) The valley bottom of the Cha di Morte is filled by alluvial conglomeratic deposits that are incised by modern streams. Sediments that have been sampled can possibly be recycled material from paleo conglomeratic deposits. They therefore record paleo-denudation rates that can differ from the modern ones.

Additionally, these conglomerates are also probably characterized by a higher erodibility than pristine lava flow basalts. (ii) The Cha di Morte depression has been interpreted as being close to the source of the older volcanic magma and numerous dike and sills are intruding into the Dos Bodes and Cha di Morte formations [START_REF] Holm | Sampling the Cape Verde mantle plume: evolution of melt compositions on Santo Antão, Cape Verde Islands[END_REF]. Dikes can significantly influence basin dissection by compartmentalizing the volcanic structure into individual blocks and increase water infiltration (Hildenbrand et al., 2005). Thus, in the Cha di Morte depression, the relative high denudation rates can also possibly be caused by an increase of the chemical alteration or be associated with physical erosion through mass-wasting of individual blocks.

Denudation rates in other volcanic islands in the world

In Cape Verde, our denudation rates results obtained for Santo Antão, are ranging between ~3 and 60 m/Ma and are comparable to those calculated on other volcanic islands characterized by different climatic patterns. In Kaua 'I Island, Hawaii, Ferrier et al (2013) computed denudation rates based on eroded volume that range between 3 and 112 m/Ma. For the same island, [START_REF] Gayer | Spatial variability of erosion rates inferred from the frequency distribution of cosmogenic 3He in olivines from Hawaiian river sediments[END_REF] proposed a denudation rate of 56 m/Ma from cosmogenic 3 He measured in 26 samples across the Waimea basin. In Gran Canaria, Canary Islands, denudation rates based on eroded volume range between 0.1 and 240 m/Ma [START_REF] Menéndez | Fluvial dissection, isostatic uplift, and geomorphological evolution of volcanic islands (Gran Canaria, Canary Islands, Spain)[END_REF]. One extreme value reaches 975 m/Ma, however, it has been measured in a small basin of 2.3 km 2 . We decided to discard this value as well as 14 other denudation rates that have been computed in basins smaller than 5 km 2 . This avoids most of the outliers that seem to be related to local effects related to small eroded volumes in modest catchments. Finally, in Tahiti-Nui, in French Polynesia, Hildebrand et al. (2008) calculated denudation rates between 72 and 2700 m/Ma. Results from the different volcanic islands show significant differences in the range and absolute value of the observed denudation rates. Despite the different techniques used to infer denudation rates, previous studies and our study can be directly compared as all techniques measure both the mechanical and chemical erosion [START_REF] Kirchner | Cosmogenic nuclide methods for measuring longterm rates of physical erosion and chemical weathering[END_REF][START_REF] Ferrier | Covariation of climate and long-term erosion rates across a steep rainfall gradient on the Hawaiian island of Kaua 'i[END_REF]. Remarkably, we observe a logarithmic correlation between the Quaternary denudation rates and the mean annual precipitation rates when considering the data of the four islands (R 2 = 0.63, p-value < 0.001; Fig. 9). Although these different islands have variable volcanic morphologies and lithology, resulting from different geodynamic contexts, climatic conditions seem to be a predominant control on the denudation rates, and therefore, on landscape evolution not only in Santo Antão but also in other volcanic islands in general. Our results demonstrate the strong coupling through which climate influences the variability of denudation and show how climate acts to shape the landscape of volcanic islands.

Conclusion

We accurately measured long-term (over 10 ka to 200 ka timescale) denudation rates from cosmogenic 3 He in pyroxene from 23 river bedload samples in Santo Antão Island. The calculation of the long-term denudation rates allows the reconstruction of a volcanic island landscape and landform history in response to a contrasted climate. Exposed to the trade winds and associated relatively heavy rainfall, the northeast part of the island exhibits high denudation rates and deeply incised basins. The southern part of the island that currently receives low precipitation displays low denudation rates and gently southward dipping slopes. Accordingly, the denudation rates and the morphological evolution of the Santo Antão Island is probably primarily controlled by climatic processes. After comparison with previously published denudation rates from other islands, our data shows that the large climatic gradients existing on volcanic ocean islands have a predominant role on their landscape evolution. Our data represent additional inputs for construction of the global database of long-term denudation rates affecting volcanoes within tropical settings. The development of such database is key for better quantifying chemical and mechanical erosion processes in tropical environments. 3 Hec calculated with the first approach (Eq. 4) versus the 3 Hec calculated with the second approach (Eq. 5). b. Denudation rates calculation from the first approach versus the denudation rates calculation from the second approach. The grey band represents the single prediction band of the one standard deviation of the correlation trend in red.

Figure captions

Figure 5: Comparison of the replicate measures of the denudation rates. a. Comparison of the replicates of the denudation rates calculated with the first approach. b. Comparison of the replicates of the denudation rates calculated with the second approach. c. Comparison of the replicates of the denudation rates calculated with the two approaches together. D. Zoom of the Fig. 5c showing the best reproducibility of the measurements for denudation rates lower than ~15 m/Ma. The square symbols correspond to computed denudation rates based on eroded volume, whereas the triangles correspond to denudation rates based on 3 He cosmogenic concentration.

Figure 1 .

 1 Figure 1. Overview of the Cape Verde islands and ASTER satellite images of Santo Antão Island.

Figure 2 .

 2 Figure 2. Simplified geological map of Santo Antão modified from List et al. (2007).

Fig 3

 3 Fig 3. a. Precipitation and elevation maps over Santo Antão Island from TRMM. The black arrows show the main direction of the wind. b. Photograph of the dry region of Tarrafal on the southwest side of the island. c. Photograph of the wetter valley of the Ribeira on the northeast side of the island.

Fig

  Fig 4: a.3 Hec calculated with the first approach (Eq. 4) versus the 3 Hec calculated with the second approach (Eq. 5). b. Denudation rates calculation from the first approach versus the denudation rates calculation from the second approach. The grey band represents the single prediction band of the one standard deviation of the correlation trend in red.

Figure 6 :

 6 Figure 6: Spatial variation of the denudation rates on Santo Antão Island.

Figure 7 :

 7 Figure 7: Map of the slope distribution on Santo Antão and hillshade map on the left top side. The north and northeast side of the island as well show high slopes corresponding to the dissected valleys. The edges of the Cha di Morte also show high slope angles, whereas the south of the island shows relatively low slope angles corresponding to a smooth morphology.

Figure 8 :

 8 Figure 8: Log-ratio transformed compositional biplots derived from principle component analysis. The length of the black rays is proportional to the amount of variation explained by the corresponding controlling factors. Samples taken from different basins are clustered according to their calculated denudation rates with the second approach.

Figure 9 :

 9 Figure9: Denudation rates measured on different volcanic islands in the world using logarithmic scales.

Figure

  Figure 1 (Color)

  Fig 3b.

  Figure 3 (Color)

  Table A.1. Annual precipitation rates have been estimated from the TRMM model with 20% of uncertainties and gradient information have been extracted from a 30 m resolution DEM (Table A.1). Helium results of the pyroxenes/olivines crushing are presented in Table A.1. Major and trace element concentrations used to calculate the 4 Herad and 3 Henucl components are given in supplementary Table A.2. Results of high temperature extractions and calculated 3 Hec concentration and denudation rates are presented in Table

Acknowledgments

The manuscript has been read by all authors, who have approved the content and agree to its publication.

This study was supported by the INSU-CNRS SYSTER. We thank Emmanuel Davy for the mineralogical separation of the samples and Laurent Zimmerman for his assistance and help during He measurements.

Denudation rates (m/Ma) Mean annual precipitation rates (mm/yr) 500 5000 Tahiti-Nui, French polynesia (Hildebrand et al., 2008) Kaua'i, Hawaii [START_REF] Ferrier | Covariation of climate and long-term erosion rates across a steep rainfall gradient on the Hawaiian island of Kaua 'i[END_REF] Kaua'i, Hawaii [START_REF] Gayer | Spatial variability of erosion rates inferred from the frequency distribution of cosmogenic 3He in olivines from Hawaiian river sediments[END_REF] Santo Antao, Cape Verde (This study) Gran Canaria, Canary Islands (Menendez et al., 2008) 50 1


Bin Cheng

Yang Yu

Erik Asphaug

Patrick Michel

Derek C Richardson

Masatoshi Hirabayashi

Makoto Yoshikawa

Hexi Baoyin
email: baoyin@tsinghua.edu.cn

Reconstructing the formation history of top-shaped asteroids by the surface boulders distribution

Finding the basic mechanism governing the surface history of asteroids of various shapes is essential for understanding their origin and evolution. In particular, the currently visited asteroids (162173) Ryugu 1 and (101955) Bennu 2 by Hayabusa2 and OSIRIS-REx appear to be top-shaped. This distinctive shape, characterized by a raised equatorial bulge, is shared

by other similarly sized asteroids, including Didymos A 3 , 2008 EV5 [START_REF] Busch | Radar observations and the shape of near-Earth asteroid 2008 EV5[END_REF] and 1999 KW4 Alpha [START_REF] Ostro | Radar imaging of binary near-Earth asteroid (66391) 1999 KW4[END_REF] . However, the possibly common formation mechanism that causes the top-like shape is still under debate. One clue may lie in the boulders on their surfaces. The distribution of these boulders, that was precisely measured in unprecedented detail by the two spacecrafts [START_REF] Watanabe | Hayabusa2 arrives at the carbonaceous asteroid 162173 Ryugu-A spinning top-shaped rubble pile[END_REF][START_REF] Lauretta | The unexpected surface of asteroid (101955) Bennu[END_REF] , constitutes a record of the geological evolution of the surface regolith since the origin of these asteroids. Here, we show that during the regolith migration driven by YORP spin-up [START_REF] Walsh | Rotational breakup as the origin of small binary asteroids[END_REF][START_REF] Hirabayashi | Internal structure of asteroids having surface shedding due to rotational instability[END_REF][START_REF] Sánchez | Rotational evolution of self-gravitating aggregates with cores of variable strength[END_REF][START_REF] Sánchez | Cohesive regolith on fast rotating asteroids[END_REF] , the surface boulders coevolve with the underlying regolith and exhibit diverse dynamical behaviors: they can remain undisturbed, sink into the regolith layer and become tilted, or totally be buried by the downslope deposition, depending on their latitudes. The predominant geological features commonly observed on top-shaped asteroids, including the boulder-rich region near the pole [START_REF] Watanabe | Hayabusa2 arrives at the carbonaceous asteroid 162173 Ryugu-A spinning top-shaped rubble pile[END_REF][START_REF] Michikami | Boulder size and shape distributions on asteroid Ryugu[END_REF] , the deficiency of large boulders at the equatorial area [START_REF] Michikami | Boulder size and shape distributions on asteroid Ryugu[END_REF][START_REF] Dellagiustina | Properties of rubble-pile asteroid (101955) Bennu from OSIRIS-REx imaging and thermal analysis[END_REF] and partially buried, oblique boulders exposed on the regolith surface [START_REF] Sugita | The geomorphology, color, and thermal properties of Ryugu: Implications for parent-body processes[END_REF][START_REF] Walsh | Craters, boulders and regolith of (101955) Bennu indicative of an old and dynamic surface[END_REF] , are commensurate with this coevolution scenario. The surface regolith migration thus is the prevalent mechanism for the formation history of the top-shaped asteroids with stiffer cores.

Our investigation of the quasi-static process of boulder-regolith evolution by YORP spinup was performed using numerical code based on the soft-sphere discrete element method [START_REF] Cundall | A discrete numerical model for granular assemblies[END_REF][START_REF] Sánchez | Simulating asteroid rubble piles with a self-gravitating softsphere distinct element method model[END_REF][START_REF] Schwartz | An implementation of the soft-sphere discrete element method in a high-performance parallel gravity tree-code[END_REF] (SSDEM). The considered structure and surface of a top-shaped asteroid consist of a "two-layer structure" with a shallow mantle of loose granular regolith covering a stiffer interior [START_REF] Hirabayashi | Internal structure of asteroids having surface shedding due to rotational instability[END_REF][START_REF] Sánchez | Rotational evolution of self-gravitating aggregates with cores of variable strength[END_REF] . This model is supported by multiple observations of the surface geologic features on Ryugu and Bennu, e.g., the longitudinal ridges and long linear grooves denoting internal stiffness and structural rigidity with a weak strength of less than a few Pa 17 ; the measured thermal inertia and substantial material flows evidencing loose and unconsolidated surface regolith [START_REF] Watanabe | Hayabusa2 arrives at the carbonaceous asteroid 162173 Ryugu-A spinning top-shaped rubble pile[END_REF][START_REF] Walsh | Craters, boulders and regolith of (101955) Bennu indicative of an old and dynamic surface[END_REF] . The pit feature of the artificial impact crater on Ryugu provides further evidence of a stronger substrate covered by cohesionless regolith [START_REF] Arakawa | An artificial impact on the asteroid (162173) Ryugu formed a crater in the gravity-dominated regime[END_REF] . This heterogeneity could appear as a result of stronger interparticle cohesion, grains of different materials, higher packing fraction, a big inner monolith, etc., or a combination of those particular characteristics [START_REF] Sánchez | Rotational evolution of self-gravitating aggregates with cores of variable strength[END_REF] . The "granular shell" model is shown in Fig. 1a, i.e., the asteroid interior is modelled as a rigid spheroid (∼ 400 m in radius), and a longitudinal slice of the regolith layer (∼ 30 m in thickness) is considered for the SSDEM simulation, consisting of 2,063,044 pebble-size particles. The boulders of various sizes, modelled as massive polyhedra, are located at four latitudes {15 • , 45 • , 75 • , 87 • } as representatives of low, mid-low, mid-high, and high latitudes, respectively. We assumed these boulders are initially on top of the regolith layer, which could be the direct result of the re-accumulation that formed the asteroids after the catastrophic disruption of their parent bodies [START_REF] Sugita | The geomorphology, color, and thermal properties of Ryugu: Implications for parent-body processes[END_REF][START_REF] Walsh | Craters, boulders and regolith of (101955) Bennu indicative of an old and dynamic surface[END_REF] , or much later by ballistic sorting [START_REF] Wright | Boulder stranding in ejecta launched by an impact generated seismic pulse[END_REF] or Brazil-nut effect [START_REF] Matsumura | The brazil nut effect and its application to asteroids[END_REF] during impactinduced seismic events. Two different types of regolith grains were used: one with moderate friction (internal friction angle of ∼ 25 • ) and one with gravel-like friction (internal friction angle of ∼ 30 • ). We then used two types of regolith grains for each of the four latitudes, and each of the three boulders, giving a total of 24 simulations (all simulation parameters are summarized in the supplementary information).

The time evolution of this regolith migration during the spin-up process is shown in Fig. 1c.

The regolith remains stable at first. When the progenitor reaches the critical spin period T c = 4.35 hr, a localized landslide first occurs at mid-latitudes [START_REF] Sánchez | Cohesive regolith on fast rotating asteroids[END_REF][START_REF] Scheeres | The dynamic geophysical environment of (101955) Bennu based on OSIRIS-REx measurements[END_REF] . Note that T c here is for surface landslide failure, which is different from the shedding limits and internal failure limits [START_REF] Hirabayashi | Internal structure of asteroids having surface shedding due to rotational instability[END_REF][START_REF] Sánchez | Rotational evolution of self-gravitating aggregates with cores of variable strength[END_REF] . The regolith mi-gration forms a narrow zone with intense shear strain (Fig. 1c, first panel). Within the shear band, particle contacts break and remake. Such rearrangement events soften the shear strength of the granular material [START_REF] Reddy | Evidence of mechanically activated processes in slow granular flows[END_REF][START_REF] Germanovich | Dynamic growth of slip surfaces in catastrophic landslides[END_REF] , i.e., the particles in this zone rotate and slip, transporting the regolith above this weakened band downslope along the failure surface like a frozen ice block, with an average creeping speed of less than 1 mm/s (Fig. 1b). Regolith above the slip region remains stable and undisturbed due to the low surface slope near the pole [START_REF] Sánchez | Cohesive regolith on fast rotating asteroids[END_REF] . With increasing spin rate Ω, the scarp at the head of the slip zone retreats towards high latitudes, leading to the shrinking of the stable region, as shown by the kinetic energy distribution change in Supplementary Fig. 3, in parallel with the massive mass wasting towards low latitudes. The retreat of the headscarp to high latitudes and the expansion of the depletion region at mid-latitudes are both evident from the growth of the failure surface (Fig. 1c, second panel). In the later stage of the migration, due to the existence of a stable region near the equator (low surface slope), the barely moved grains here increase the resisting force along the failure surface in the landslide foot, blocking the material in the head from moving towards the foot. Thus, the creeping debris moves uphill and then overlies the pre-existing materials, forming a distinct deposit bulge near the equator. The original shear band propagates to low latitudes and subsequently creates a rupture surface inside the deposition (Fig. 1c, third panel).

It is notable that the elevated geologic upheaval that marks the foot of the landslides is not at the equator at first, i.e., the collapse is limited to the mid-latitudes due to the resistance force between regolith grains. Then, as the accumulation of downslope materials from the depletion zone goes on, this bulge swells upward, leading to a steeper local slope that eventually collapses towards the equator. The raised, pronounced equatorial ridge denotes a fully sped-up progenitor for which the top-shaped morphology emerges (Fig. 2a). The simulated topography, with low elevation at mid-latitudes, mean elevation at the polar region and high elevation at the equator, resembles the observed shapes of typical top-shaped asteroids (Fig. 2b), which is consistent with previous findings that regolith migrations are a major force that has reshaped these asteroids if the interior is stronger [START_REF] Hirabayashi | Internal structure of asteroids having surface shedding due to rotational instability[END_REF][START_REF] Sánchez | Rotational evolution of self-gravitating aggregates with cores of variable strength[END_REF][START_REF] Sánchez | Cohesive regolith on fast rotating asteroids[END_REF] . Simulations with higher friction give a similar picture, except that the stable region near the pole is larger and that the equatorial ridge is steeper (Supplementary Video 6).

The motion of boulders across the surface shows three distinct dynamical behaviours (see Supplementary Videos 2-5). The 15 • and 45 • boulders are located in the early migration zone, thus they start to creep downslope in synchronization with the first landslide of granular regolith.

During the creeping-down process of the 45 • boulder, the strength of its supporting force chain network is weakened due to the grain agitations induced by the shear deformation [START_REF] Nichol | Flow-induced agitations create a granular fluid[END_REF] . These weak grain bonds generate a fluid-like behaviour; thus this boulder sinks into the regolith layer by a few metres and becomes tilted, reminiscent of the abundant partially buried boulders on Ryugu and Bennu [START_REF] Sugita | The geomorphology, color, and thermal properties of Ryugu: Implications for parent-body processes[END_REF][START_REF] Walsh | Craters, boulders and regolith of (101955) Bennu indicative of an old and dynamic surface[END_REF] . This flow-induced submergence phenomenon has also been observed in granular experiments [START_REF] Reddy | Evidence of mechanically activated processes in slow granular flows[END_REF][START_REF] Nichol | Flow-induced agitations create a granular fluid[END_REF] . At lower latitudes, the 15 • boulder is first elevated by the extruded underlying sediment, subsequently transferred downslope to the foot of the deposit upheaval, and eventually buried entirely by the accumulated debris from the collapsed bulge hill. This mechanism, i.e., boulders initially located at low latitudes being totally buried by the downslope deposition, may have caused the deficiency of large boulders at the equator as observed on Ryugu 10 and Bennu [START_REF] Walsh | Geology of Bennu's equatorial ridge[END_REF] .

At higher latitudes, the 75 • boulder remains static until the regressive headscarp of the depletion zone reaches 75 • latitude. The growth of this scarp removes support from adjacent grains, forming a fragile force-chain network downstream (Fig. 3). The weak support increases the equivalent local slope and thus triggers a downslope migration of the boulder. Finally, this boulder sinks into the granular mantle and leaves only a partially exposed body above the regolith layer, generating a distinct outcrop-like landform (Fig. 4). The 87 • boulder, however, is located in a stable region (see Fig. 2a) and thus does not move throughout the reshaping process of the top-shaped asteroids. This result is consistent with the geologic indication of a stable pole deduced from the disorganized boulder structure near Bennu's pole and the non-slipping polar boulder Otohime on Ryugu [START_REF] Watanabe | Hayabusa2 arrives at the carbonaceous asteroid 162173 Ryugu-A spinning top-shaped rubble pile[END_REF][START_REF] Sugita | The geomorphology, color, and thermal properties of Ryugu: Implications for parent-body processes[END_REF] .

We infer that Otohime could be a primordial relic of the original disruption/reaccumulation of Ryugu's parent body. Future studies on the geologic age of Otohime from thermal or fracture properties would provide significant constraints on when the catastrophic impact event that formed Ryugu occurred. Note that we do not rule out the possibility that other exhumation processes could contribute to the subsequent evolution of these boulders. Some of the buried boulders in the equatorial region on Ryugu could be exhumed again by migrating grains from the ridge towards the current geopotential lows at higher latitudes [START_REF] Sugita | The geomorphology, color, and thermal properties of Ryugu: Implications for parent-body processes[END_REF] . Other mechanisms including Brazil nut effect [START_REF] Matsumura | The brazil nut effect and its application to asteroids[END_REF] and ballistic sorting [START_REF] Wright | Boulder stranding in ejecta launched by an impact generated seismic pulse[END_REF] could also further modify these landforms.

Our results show that the top shapes of rubble pile asteroids can be evolved via regolith migrations from mid-latitudes towards the equator driven by the YORP spin-up effect. If the interior is stronger than the outer regolith, then miniature surface landslides could be frequent as the spin rate approaches the critical spin limit [START_REF] Sánchez | Cohesive regolith on fast rotating asteroids[END_REF][START_REF] Scheeres | The dynamic geophysical environment of (101955) Bennu based on OSIRIS-REx measurements[END_REF] . These landslides are to be distinguished from their counterpart on large-gravity bodies like Mars, the Moon or Earth. On such planets the large boulders usually bounce and deposit near the toe of the avalanche sediment [START_REF] Gray | Particle segregation in dense granular flows[END_REF] . By contrast, regolith migra-tion on low gravity rubble piles, driven by YORP, once near the spin limit, creeps imperceptibly at glacial speeds. A detailed look into these creep-like landslides shows that the boulders are in a "co-evolving" relationship with the underlying granular regolith, which may account for the unique geological landforms observed on top-shaped asteroids, i.e., the stable boulders near the pole, the abundant partially buried boulders and the deficiency of large boulders at the equator. We find consistency of the derived shapes and boulder features between the observations and the numerical results, which lends support to our model as a way to decipher the evolution history of top-shaped asteroids. Further simulations show that the large craters that overlay the equator of Ryugu would have been eroded by the migrating regolith if they were formed before the equatorial ridges (Supplementary Fig. 7), implying an early-stage formation of the overall shapes [START_REF] Watanabe | Hayabusa2 arrives at the carbonaceous asteroid 162173 Ryugu-A spinning top-shaped rubble pile[END_REF] . This would appear to be in tension with the long timescales over which YORP effect could significantly modify the spins [START_REF] Michel | Collisional formation of top-shaped asteroids and implications for the origins of Ryugu and Bennu[END_REF] . A plausible explanation is that the progenitor bodies that emerged from the catastrophic disruptions are already oblate shapes with relatively high spin rates, which is a possible result of the direct reaccumulation of fragments in the case of large starting angular momentum [START_REF] Michel | Collisional formation of top-shaped asteroids and implications for the origins of Ryugu and Bennu[END_REF] . These early-formed "diamond-shaped" rubble piles with high spin rates could serve as a midway starting point that facilitates the formation of top shapes, i.e., the YORP-induced migrations further reshaped the bodies and remodeled their surface features like large boulders to the current geomorphology. Another possibility is these craters could be created by detachment of a small chunk from the equatorial ridge at a more rapid rotation state [START_REF] Tardivel | Equatorial cavities on asteroids, an evidence of fission events[END_REF] . This mechanism seems to be plausible given that the past spin period of Ryugu could be up to around 3 hours as deduced from its east-west dichotomy [START_REF] Hirabayashi | The western bulge of 162173 Ryugu formed as a result of a rotationally driven deformation process[END_REF] . In such a scenario, the circum-equatorial ridge and overlying craters could both be landscapes sculpted by YORP spin-up.

Methods

We performed discrete element simulations [START_REF] Cundall | A discrete numerical model for granular assemblies[END_REF][START_REF] Sánchez | Simulating asteroid rubble piles with a self-gravitating softsphere distinct element method model[END_REF][START_REF] Schwartz | An implementation of the soft-sphere discrete element method in a high-performance parallel gravity tree-code[END_REF] using an original soft-particle N-body code DEMBody [START_REF] Cheng | Numerical simulations of the controlled motion of a hopping asteroid lander on the regolith surface[END_REF] , which is capable of tracing the quasi-static behaviour of particle-particle and particle-polyhedron interactions. In this code, the nonlinear contact forces between two contacting objects (spheresphere, sphere-polyhedron, sphere-wall) are calculated based on Hertz-Mindlin contact theory [START_REF] Somfai | Elastic wave propagation in confined granular systems[END_REF] according to the mutual overlap (typically <1% of their radii). A dimensionless coefficient µ is used to control the stick-slip friction between colliding particles. A physically based rotational resistance model incorporating rolling and twisting friction was implemented into this code, parameterized by a quantity β that represents a statistical measure of real particle nonsphericity [START_REF] Cheng | Numerical simulations of the controlled motion of a hopping asteroid lander on the regolith surface[END_REF][START_REF] Zhang | Rotational failure of rubble-pile bodies: Influences of shear and cohesive strengths[END_REF] .

To calculate the contact between non-spherical boulders and grains, the boulder is represented as an assembly of movable walls with the same geometry and inertia tensor [START_REF] Cheng | Numerical simulations of the controlled motion of a hopping asteroid lander on the regolith surface[END_REF] , which allows it to behave as a rigid body. A second-order leapfrog integrator is used for integration of the equations of motion of both particles and boulders.

Our simulation scheme is as follows. First, we prepared a polydisperse granular medium consisting of 2,063,044 cohesionless particles with properties chosen according to Ryugu's regolith 1 , i.e., material density of 2.4 g/cc. A slight dispersion in the grain radius was considered, with a random distribution ranging from 0.4 m to 0.6 m to avoid any crystallization of the packing. The particles were dropped onto a rigid, frictional shell with a radius of 400 m and width of 100 m confined by two vertical frictionless walls [START_REF] Sánchez | Cohesive regolith on fast rotating asteroids[END_REF] , as shown in Fig. 1a, forming a longitudinal slice of surrounding walls has only a minor influence on the boulder dynamics for our simulation settings. By using a granular bed of both hemispheres, which consists of 4,126,088 particles, we tested the symmetry of the regolith migration. Simulations show the sliding debris meet and accumulate symmetrically at the equator (Supplementary Video 9). The overall geomorphology is similar to that in hemisphere cases.

Equatorial crater. We investigated the evolution of equatorial craters during the YORP spin-up.

A crater with radius of 80 m was generated at the equator before rotational acceleration (Supplementary Fig. 7). Simulations show that the migrating debris would have filled the equatorial craters if these craters formed before the equatorial ridges, which implies an early-stage formation of the overall morphology.

Boulders configuration. We investigated the dependency of the boulder evolution on their configurations. Simulations with boulders ten meters below the regolith surface show these initially underground boulders keep buried throughout the creeping evolution (Supplementary Videos 10-11). We also tested the interaction between multiple boulders. Given that the number density of boulders 10 larger than 15 m on Ryugu is ∼ 100 km -2 , we randomly placed eight boulders, either on the surface or embedded in regolith, at middle and low latitudes of the simulated longitudinal slice. Simulations give results similar to those described above, in which surface boulders sink into or are buried by the creeping regolith, while underground boulders remain buried (Supplementary Video 12).

Figure 1

 1 Figure 1 Asteroid morphological change during the YORP spin-up process. a, "Two-layer structure" with a shallow mantle of loose granular regolith covering a stiffer interior. A large boulder lies on the regolith surface. b, Spin-up path and consequent creep speed of global regolith migrations for a 75 • boulder with moderate-friction regolith case. c, Four selected frames of the regolith and boulder at different spin periods showing the regolith migration under the YORP effect. The cuboid represents the boulder and is located at 75 • latitude initially. The regolith has moderate friction. When the progenitor reaches the critical spin period T c = 4.35 hr, the centrifugal tension induces a narrow zone with intense shear strain (dashed line in first panel; see Supplementary Fig. 2 for the shear heat of regolith particles). This shear band expands towards high latitudes, and creeping debris moves to low latitudes along the rupture surface in the deposition (dashed line in third panel). At the end of the simulation, deposit sediment accumulates at low-latitudes, forming a pronounced ridge along the equator (fourth panel). The failure surface acts as a dynamical divider separating the static zone from the depletion zone, whose growth controls the landslides across the regolith bed. Animation of this figure: Supplementary Video 1.

Figure 2

 2 Figure 2 The simulated geomorphology. a, The topography changes during the YORP spin-up process, showing the barely moved pole and migrating debris from midlatitudes overlying the pre-existing materials near the equator. Here the elevation is the radial height of particles above the interior shell. Data are for a moderate-friction case.b, Simulated morphological profile (black line) compared with Ryugu's (blue) 1 , Bennu's (red)[START_REF] Scheeres | The dynamic geophysical environment of (101955) Bennu based on OSIRIS-REx measurements[END_REF] , 1999 KW4's (green)[START_REF] Ostro | Radar imaging of binary near-Earth asteroid (66391) 1999 KW4[END_REF] and 2008 EV5's (yellow)[START_REF] Busch | Radar observations and the shape of near-Earth asteroid 2008 EV5[END_REF] average shapes. All profiles are normalized by the corresponding equivalent radius. Shaded areas stand for the 1σ error bands.

Figure 3

 3 Figure 3 Visualization of the force networks excited by the creeping boulder. Colourindicates the elastic energy, E, stored in each particle relative to the average elastic energy, E 0 . Note that only those particles with E > E 0 are considered part of the force chains[START_REF] Peters | Characterization of force chains in granular material[END_REF] . The white-curve indicates the free surface of the creeping regolith. Red arrows show the sinking and tilting directions of the boulder. The chain structure downstream is significantly fragmented and spatially sparse due to the removal of adjacent grains by the retrogressive landslide, leading to subduction migration of the boulder. Data are for a 75 • boulder with moderate-friction regolith case.

Figure 4

 4 Figure 4 Typical landforms of boulders sculpted by YORP evolution. The 87 • boulder remains stable, consistent with the non-slipping polar boulder Otohime on Ryugu 12 . The boulders at mid-latitudes sink into the regolith layer by a few metres and become tilted. These landforms generated by regolith migrations are likely the geophysical origins of the partially buried/oblique boulders on Ryugu 12 and Bennu 13 . (panel d, hyb2 onc 20180831 101059 panel e, hyb2 onc 20180801 160642 tnf l2c and panel f, hyb2 onc 20180720 075208 tvf l2b are from Hayabusa2 ONC dataset 34 .)

Acknowledgements This work is funded by the National Science Fund for Distinguished Young Scholars of China (No. 11525208). Y.Y. acknowledges support from the Natural Science Foundation of China (No. 11702009). P.M. acknowledges funding from the French space agency CNES, from the European Union's Horizon 2020 research and innovation programme under grant agreement No. 870377 (project NEO-MAPP) and from Academies of Excellence: Complex systems and Space, environment, risk, and resilience, part of the IDEX JEDI of the Université Côte d'Azur. The authors acknowledge supports from JSPS Core-to-Core program "International Network of Planetary Sciences". We are grateful to the Hayabusa2 ONC team for providing images which are available at the JAXA Data Archives and Transmission System (DARTS) at http://darts.isas.jaxa.jp/pub/hayabusa2/onc_bundle.

Data availability. The data that support the plots within this paper and other findings of this study are available from the corresponding authors upon reasonable request.

Code availability. The code used to generate the datasets is available from the corresponding author on reasonable request.

the asteroid surface with a thickness of ∼ 30 m. A cuboid boulder, the size of which is similar to that of typical boulders on Ryugu, was initially located above the interface of the granular media.

The whole system is under the two-component gravity [START_REF] Sánchez | Cohesive regolith on fast rotating asteroids[END_REF] exerted by the rocky interior (bulk density of 1.2 g/cc) 1 and the regolith shell. After a settling-down phase of the granular-polyhedron system, we started to accelerate the asteroid linearly following the prescribed spin-up path shown in Fig. 1b. The acceleration was performed until massive debris shedding occurred at the equator.

The angular acceleration used in this study is 2.5 × 10 -10 rad/s 2 , which is much smaller than the reference values used in previous works [START_REF] Sánchez | Rotational evolution of self-gravitating aggregates with cores of variable strength[END_REF][START_REF] Zhang | Rotational failure of rubble-pile bodies: Influences of shear and cohesive strengths[END_REF] . By using this slow spin-up acceleration, we confirm that the motion of the granular grains would immediately stop at any point in time when the spinup is stopped, which means that the regolith layer maintains quasi-static equilibrium during the spin-up process (see Supplementary Information).

Regolith bed geometry. We investigated the dependency of the boulder evolution on the geometry of the regolith bed. By using granular beds with depth of 45 m and of 20 m, which consists of 3,246,524 particles and 1,452,565 particles, respectively, we tested the influence of regolith bed depth (Supplementary Figure 45). Simulations give results similar to those described above, in which tilted boulders and buried boulders emerge with the surface migration driven by YORP spin-up. By using a granular bed with width of 160 m, which consists of 3,332,716 particles, we tested the influence of confining sidewalls (Supplementary Video 7). Simulations indicate the size of surrounding walls is large enough to eliminate the boundary effect. By using a spherical wedgeshaped bed with angular width of 20 • , which consists of 2,063,045 particles, we tested the influence of the boundary shape (Supplementary Figure 6). Simulations indicate the shrinkage effect from