Neogene fluvial landscape evolution in the hyperarid core of the Atacama Desert

Dating of extensive alluvial fan surfaces and fluvial features in the hyperarid core of the Atacama Desert, Chile, using cosmogenic nuclides provides unrivalled insights about the onset and variability of aridity. The predominantly hyperarid conditions help to preserve the traces of episodic climatic and/or slow tectonic change. Utilizing single clast exposure dating with cosmogenic 10Be and 21Ne, we determine the termination of episodes of enhanced fluvial erosion and deposition occurring at ~19, ~14, ~9.5 Ma; large scale fluvial modification of the landscape had ceased by ~2–3 Ma. The presence of clasts that record pre-Miocene exposure ages (~28 Ma and ~34 Ma) require stagnant landscape development during the Oligocene. Our data implies an early onset of (hyper-) aridity in the core region of the Atacama Desert, interrupted by wetter but probably still arid periods. The apparent conflict with interpretation that favour a later onset of (hyper-) aridity can be reconciled when the climatic gradients within the Atacama Desert are considered.

Cosmogenic nuclide exposure dating Samples were ground and sieved to 250-1000 µm and subsequently purified by sequential HFleaching 2 . Purified quartz separates were investigated under a microscope. Those with a high abundance of visible fluid inclusions and/or fragments resembling chalcedony were further etched or excluded for 21 Ne. Inductively coupled plasma optical emission spectrometry (ICP-OES) was used to verify the purity of the quartz before dissolution. 10 Be target for accelerator mass spectrometry (AMS) were prepared using the stacked column approach outlined in Binnie et al. 3 . A reagent blank was prepared in tandem with the sample batch 10 Be/ 9 Be AMS measurements were made on CologneAMS 4 normalized to the ICN standard dilutions prepared by Nishiizumi 5 . Concentrations of 10 Be are reported following blank subtraction, which was less than 1% of the total number of nuclides measured in each case. The 1 standard deviation analytical precision of the nuclide concentrations was estimated by summing in quadrature the relative uncertainties on the AMS measurements of both the samples and the blank, along with a 1% (1 s.d.) estimate for the precision of the mass of 9 Be added during spiking.
For cosmogenic 21 Ne, cleaned and purified samples were sieved to <125 µm and packed into aluminium foil cups. The samples were measured with a noble-gas mass spectrometer at Scottish Universities Environmental Research Centre (SUERC) applying the standard procedure [6][7][8] .
For age calculation, we used a 10 Be half-life of 1.36 ± 0.07 Ma 5 and a 'nuclide dependent scaling' after Lifton, et al. 9 , calculated with "The online exposure age calculator formerly known as the CRONUS-Earth online exposure age calculator." Version 3, http://hess.ess.washington.edu/math/v3/v3_age_in.html; 10 . This employed a high-latitude, sealevel spallogenic 10 Be production rate of 3.92 ± 0.31 atoms g -1 a -1 11 and for 21 Ne 18.7 ± 2.3 atoms g -1 a -1 12 and ERA40 atmospheric model. We applied a time-integrated inverse modelling of cosmogenic nuclide production rates by subsequent modelled subsidence to paleo elevations. We used a linear uplift model of 40 m/Ma 13,14 with uplift commencing at 20 Ma (assuming that uplift of the alluvial fans started from a vertical position near sea-level). To apply uplift correction using LSDn scaling as implemented the CRONUS-Earth online calculator (V.3), we calculated the uplift-dependant change by using an erosion rate (the reduction of the mass air above the sample by uplift is mathematically equivalent to erosion).
The density of air decreases with altitude (rising from sea-level to 1100 m, our mean sample elevation) the density of air is reduced by~10% (U.S. Standard Atmosphere 1976, NASA-TM-X-74335). For the purpose of calculating the effect of atmospheric mass removal by uplift we utilize the mean density of air at 500 m.a.s.l. (1.167 kg/m 3 , NASA-TM-X-74335). The ρair 500m/ρQtz multiplied by the uplift rate provides the 'erosion rate' (in our case 1,762 x*10 -6 cm/yr) for input into the online calculator. Since we use the density of air at 500 m at all elevations, we assign a 5% internal uncertainty for our uplift correction (caused by the utilization of a mean density of air).
For samples with concentrations implying an age >20 Ma, we use the scaling calculated for a virtual 20 Ma sample at the sampling location, to calculate the exposure prior to 20 Ma (from the concentration that is in excess of that produced since 20 Ma); which is added to 20 Ma to obtain the exposure age.
The uplift correction is relatively insensitive to deviation of the assumption of linear uplift (within the bounds of the same long-term average). Ages change by ±10% if the uplift rate decreased/increased fourfold after half the exposure time (Dunai et al. 2005). We assign a ±10% external uncertainty to the uplift-corrected ages. Ages without uplift correction would be younger. Early Miocene ages would be younger by ~30%, differences decrease towards younger ages (~2% at 2 Ma).  The zircon U/Pb age record is represented by histograms and associated probability density function curve for zircons ages between 0 and 10Ma, determined for the sampled tephra T4. Note that the scale of age axis is interrupted and changes between 3 and 5 Ma. Eruption ages (Ar/Ar) of adjacent volcanic complexes are marked as vertical bars 15 . Orange bar indicates potential pre-eruptive crystallization duration of zircons within the Purico magma chamber 16 . U/Pb zircon ages could therefore be related to the Purico or the Tatio eruptions. (B) Rank order plots (ROP) and zircon U/Pb ages as probability density function (PDF) curves. The shaded PDF curve represents the distribution for the ROP displayed. The normal distribution is interpreted to represent a single population of crystals from an individual eruption. Horizontal lines display 2σ analytical uncertainty of data from which cumulative probability curves were calculated. The age of 0.98 ± 0.04 is the weighted mean and 2σ error for the dominant peak in the distribution of 35 zircons of T4.

Additional Site Data
The easternmost fan (Fan East) has a maximum width of 1.5 km and length of 4 km, moderately dipping (1.46°) to the north. The lower 1.5 km are buried under lacustrine sediments, related to lake stages in the Quillagua-Llamara basin 17 . Diatom-rich sediments indicate deposits of the Quillagua Formation 5.5-4. 5 Ma,17,18 outcrop at the surface. Deflation removed material to the bed-rock level adjacent to the fans resulting that the fan surface is elevated approximately 5-10 m above the adjacent surface, resembling an inverted landscape topography in the lower part of the fan.
The central fan (Fan Center) has a maximum length of 3 km, including the final 1 km covered by lacustrine sediments, and a maximum width of 2.4 km. Similar like the eastern fan, the center fan is surrounded by an inverted landscape-like topography. The fan surface dips moderately to the north (1.29°). A feeder channel (RL15-005) is well preserved and has incised into blocks of the Adamito Fault System, indicating synchronous evolution/activity.
The determination of the outline of the western-most fan (Fan West) is difficult; is due to significant surface modification by fluvial erosion in the distal portions, and dust cover approaching its apex. To the west, the fan is bound by the Atacama Fault. The east is limited by deflation sinks, and the western margin of the Central Fan in its lower section. The surface dips moderately to the north (0.93°), with a maximum width of 2 km and maximum length of 4.5 km. Fan sediments are buried under lacustrine sediments in the distal part.   To increase the visibility of boundaries between alluvial fans and inverted landscape, the Pléiades 1B Multi-Spectral Image was modified using contrast enhancement. C) Pléiades 1B Multi-Spectral Image including calculated slope map using ArcGIS 10.5.1 ® Slope angles >8° are displayed in colour and indicate boundaries between alluvial fan sediments and deflation hollows. D) Enlarge image of the eastern alluvial fan surface north of the Adamito fault (Pléiades 1B Multi-Spectral Image). Visible are former topographic lows (channels with gravel infill), which are now topographic highs due landscape inversion at the eastern limit of the eastern fan. Cross-cutting on the center fan towards the eastern fan is marked.