Mass flux decay timescales of volcanic particles due to aeolian processes in the Argentinian Patagonia steppe

We investigate the timescales of the horizontal mass flux decay of wind remobilised volcanic particles in Argentina, associated with the tephra-fallout deposit produced by the 2011–2012 Cordón Caulle (Chile) eruption. Particle removal processes are controlled by complex interactions of meteorological conditions, surface properties and particle depletion with time. We find that ash remobilisation follows a two-phase exponential decay with specific timescales for the initial input of fresh ash (1–74 days) and the following soil stabilisation processes (3–52 months). The characteristic timescales as a function of particle size shows two minimum values, identified for sizes around 2 and 19–37 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\upmu$$\end{document}μm, suggesting that these size-range particles are remobilised more easily, due to the interaction between saltation and suspension-induced processes. We find that in volcanic regions, characterised by a sudden release and a subsequent depletion of particles, the availability of wind-erodible particles plays a major role due to compaction and removal of fine particles. We propose, therefore, a simple and reproducible empirical model to describe the mass flux decay of remobilised ash in a supply-limited environment. This methodology represents an innovative approach to link field measurements of multi-sized and supply-limited deposits with saltation erosion theory.

. Global distribution of major wind-remobilisation events associated with pyroclastic deposits reported in literature. Source volcanoes and eruption years are indicated. Colours show time elapsed between the occurrence of the eruption and the last remobilisation event reported. Major events are associated with several kilometres of extension, usually detected by satellite imagery. Even though they can represent a major concern for exposed communities, small events are not reported here as they are either not recorded or recorded in local or national media and are understudied. References: (a) Novarupta 1,2,11,12  Due to the prevailing winds, a vast area in the Argentinian Andes and the Patagonian steppe was covered by ∼ 1 km 3 of rhyolitic tephra [49][50][51][52] . The primary tephra-fallout deposit consists of 4 units of which the first three were deposited during the 6-15 June, 2011 and the Unit IV was deposited only in very proximal areas in Chile after 15 June, 2011. Unit III, being the most widely dispersed, upper and fine-grained, is the unit that was first remobilised 23 (Fig. 2). However, after Unit III was completely eroded, Unit I and II also started being remobilised, particularly in unprotected areas (e.g., low vegetation cover). We analyse (i) the temporal distribution of the remobilisation processes associated with the Unit III of the CC tephra-fallout deposit, (ii) the mass flux decay timescale as a function of particle size, (iii) the role of meteorological and soil conditions based on the comparison of the measured (Q msd ) and theoretical (Q th ) streamwise mass fluxes, and (iv) the influence of particle supply in volcanic erosion systems. A complementary analysis of the spatial distribution of remobilisation is available in the Electronic Supplementary Material (ESM-Section 1). We show that meteorological fluctuations, material supply, compaction of primary deposit and particle size are key parameters affecting the temporal evolution of the mass flux. Consequently, the current theoretical saltation schemes, which are based on assumptions of unlimited supply and uniform particle size, should be carefully applied to fresh volcanic deposits.

Results
Wind-remobilisation of the Cordón Caulle tephra-fallout deposit: temporal distribution. In order to study soil erosion processes, a systematic network of 4 horizontal collectors setup by INTA (National Agricultural Technology Institute of Argentina) in strategic erosion sites of the Patagonia steppe, was operated from 2011 until 2016 (Fig. 2). Each collector consists of 3 dedicated samplers at fixed heights above the ground (Fig. 3a). We calculated mass fluxes from the mass of collected material, the duration of accumulation and the cross-sectional area of the collector (See "Materials and methods" section). Temporal evolution of the mass flux is the result of a complex interaction between the surface and meteorological conditions, and is intrinsically associated with the eruptive evolution and the primary tephra sedimentation. Whilst the climactic eruptive phase started on 4 June 2011 and lasted approximately 24-30 hours, the Unit III was produced and sedimented some days after, between the 6 and 15 of June 52 (Fig. 2). There are evidences of syn-and post-eruptive remobilisation immediately after the eruption 23 ; however, it is impossible to reconstruct all the remobilisation events.
Since the largest amount of mass was collected on the site S4 and it represents the most complete sub-dataset, we focus on this single site to illustrate the main features of the temporal mass flux variation. Figure 3b shows that the mass flux decays non-monotonically over time for collector S4. The first peak at this site is recorded on Horizontal mass flux decay of volcanic particles. Total mass flux decay. Aeolian processes are strongly influenced by the supply of erodible particles, which has been shown to decrease exponentially with time 41,46,47,53 . Consequently, mass flux q as a function of time t for each height, can be described as, where τ is the decay timescale and a is a constant. A mass flux peak significantly above the background erosion was recorded after the emplacement of the tephra-fallout deposit (streamwise mass flux before the eruption for the same location was estimated as 0.09 kg m -1 day −121 ). Following this peak, we find that the evolution of the mass flux can be described by two successive phases of exponential decay with distinct timescales (Fig. 3c-e). This is true for each collector height (0.15 m, 0.50 m, 1.50 m), with inflection points between the two decay stages at 139 ± 25, 127 ± 26 and 183 ± 5 days, respectively. Decay timescale, τ (Eq. 1), is also shown for each exponential curve. The two-phase processes represent, therefore, two different timescales with the first, faster phase having shorter timescales (59 ± 16, 55 ± 16, 77 ± 16 days) than the second, slower phase (775 ± 342, 516 ± 186, 445 ± 175 days) for the 3 heights, 0.15 m, 0.50 m and 1.50 m, respectively. These results also show that the slowest decay is observed for the collector at 1.50 m for the first phase, and at 0.15 m for the second phase ( Fig. 3c-e). www.nature.com/scientificreports/ Mass flux decay as a function of particle size. Considering that tephra-fallout deposits have disperse grainsize distributions, we quantify the evolution of mass flux of different size fractions and investigate if different fractions exhibit the same two-phase mass flux decay as the total remobilised material. To this end, the mass collected at site S4 has been analysed for grainsize using half-phi-size classes. This sub-dataset corresponds to 14 collection periods at 3 different heights (see "Materials and methods" section and ESM- Table S3). We fit the mass flux of each size class at each height to a double exponential decay of the form, where τ 1 and τ 2 are the fast and slow decay timescales corresponding to phase I and phase II, respectively; b and c are constants giving the peak mass flux at t=0. Mass flux decay timescales, expressed by τ 1 and τ 2 , provide an insightful description of erosion phenomena when displayed as a function of particle size (Fig. 4). Phase I decay timescales are all within the range of ∽1-74 days (Fig. 4a), whereas phase II decay timescales are typically one-two orders of magnitude slower ( ∽3-52 months, Fig. 4b). Similarly to the whole mass (Fig. 3b), the slowest decay for the phase I and the fastest decay for the phase II are associated with the 1.50 m height for all the size fractions. Whilst mass flux decay timescales show a clear trend for the phase II at all the heights, this trend is less notable for phase I, for the 0.15 and 1.50 m measurement heights. Consequently, τ 2 shows a robust and distinct U-shaped curve at all heights with two minima values: at 2 µ m and more prominently at 19-37 µ m (Fig. 4b). In contrast, τ 1 exhibits a weak minimum from 19 to 37 µ m slightly varying with height (Fig. 4a). Following sections present how the meteorological, surface parameters and material supply can influence this temporal variation.
Role of meteorological conditions and material supply on the streamwise mass flux. The frequency and intensity of ash-remobilisation events are related to specific meteorological variables (e.g., surface wind velocity, precipitation, low-level and atmospheric humidity) and surface features (e.g., vegetation cover, soil moisture, grainsize and density of particles) 38,39,41,42 . In the Argentinian Patagonia, recurrent wind remobilisation is associated with the months of austral spring and austral summer 17,22 , due to the strong winds and gusts Error bars indicate the 68% confidence interval; uncertainties associated with the first phase are notably larger than those of the second phase. Phi and metric scales are shown; the equivalence between both scales, as well as the representative size values (used hereafter) for each bin are presented in the ESM- Table S2. Fitting procedure details are presented in the "Materials and methods" section.  (Fig. 5a). Local observations (informal interviews in Ingeniero Jacobacci, Fig. 2) suggest that events are more frequent during the afternoon, particularly after midday, which coincides with the strongest winds (> 36 km h −1 , Fig. 5a) and high surface temperatures on the ground. Maximum gust values are mainly recorded during the summer ( ∼ 90 km h −1 , Fig. 5a). The prevailing surface wind direction at surface is W to E 23 . A strong seasonality in soil moisture is evidenced by lower values of gravimetric soil moisture during spring and summer (19-26% ) than in autumn and winter (21-35% ) periods (Fig. 5b). There was a notable increase in soil moisture since the winter in 2013, peaking in spring 2014 (Fig. 5b). In addition to these meteorological conditions, the measured streamwise mass flux Q msd calculated as the vertical integral of the mass flux q is shown in the Fig. 5c. The four mass flux peaks in August 2011, December 2011, January 2013 and September 2013, that were recorded by the individual collectors (Fig. 3b) are also present in the integrated signal. In order to better understand the influence of each meteorological and soil parameter, and the importance of the availability of particles on the mass flux decays, we analyse the temporal evolution by comparing the measured Q msd and the theoretical Q th streamwise mass fluxes.
Comparison of measured and theoretical streamwise mass fluxes: impact of a finite supply of particles. In order to calculate Q th , we first use the model of Shao and Lu 44 to estimate the threshold friction velocity u * t , considering a γ of 1.65 × 10 −4 kg s −2 , which fits best the airborne material analysed in this study (ESM- Fig. S2). Secondly, we account for soil moisture through the parameterisation of Fecan et al. 45 . We combine then the resulting u * t with the wind friction velocity, u * , to calculate Q th using the saltation model of Owen 39 for all particle sizes. Two meteorological conditions have been considered in order to constrain the theoretical streamwise mass flux: (i) the precipitation rate must not be higher than 0.01 mm h −1 , according to the UK Met Office strategy 30 ; and (ii) u * must exceed u * t for each particle size (See "Materials and methods" section). www.nature.com/scientificreports/ This theoretical framework assumes that saltation occurs under an infinite supply of particles 41,47 . However, in volcanic regions characterised by an almost instantaneous release of a large amount of particles and its subsequently depletion, the supply of material is limited. This depletion creates the condition known as supplylimited saltation or source-limited saltation 41 . We account then for a correction factor on the supply-limited streamwise mass flux Q slm based on the temporal evolution of the erodibility, σ , or the potential of the surface to be eroded 41 , according to, We first assume a constant erodibility value of σ = 0.21 following the numerical simulations of Mingari 54 for the outstanding event of ash resuspension in October 2011 associated with the 2011-2012 CC eruption. Figure 6a shows the comparison of the measured Q msd for the 14 collection periods, and the theoretical Q th for a 10-day average period, considering an infinite-supply of particles. A substantial disagreement between Q msd and Q th is evident. Whilst Q th remains relatively constant with slight variations due to the meteorological fluctuations (i.e., wind speed and precipitation), Q msd decreases with time by more than two orders of magnitude from period 1 to 14. These results demonstrate that the decay of Q msd might be intrinsically related to a decreasing supply of particles, and, therefore, σ cannot be constant through time. To account for this depletion, we propose the following empirical model to estimate the temporal variability of the erodibility, σ (t) , according to, where σ 0 corresponds to the initial erodibility (here, 0.21) and α is an empirical constant parameter with units of (m kg −1 ). The Eqs. (3) and (4) were solved using an iterative approach in order to reconstruct the time series for the supply-limited streamwise mass flux Q slm . In order to determine the optimal value of α , the root mean square error between Q msd and Q slm for different values of α was computed (Fig. 6b). The best agreement between the empirical model and the measured streamwise fluxes was obtained for α = 1.3 × 10 −2 m kg −1 . The resulting Q slm , considering a time-dependent erodibility and the optimal α , for a 10-day average and for the sampling period average are shown in the Fig. 6a. There is a very good agreement between the measured and the modelled mass fluxes, except for the last sampling period where the predicted Q slm is overestimating the measured Q msd . Interestingly, this model can predict the peaks in August 2011, December 2011, January 2013 and September 2013, previously discussed (Figs. 3b, 5c and 6a).

Discussion
The temporal variation of ash-remobilisation associated with the 2011-2012 CC eruption highlights important aspects concerning aeolian processes of fresh pyroclastic deposits. Analysis of the mass flux peak values (Fig. 3b) shows that seasonal variability had a strong influence on the frequency of ash remobilisation events and, hence, on the recorded mass flux peaks. Major remobilisation (first and maximum peak) occurred some time after the emplacement of the tephra-fallout deposit associated with the main eruptive phases (4-15 June 2011). The eruption started in the middle of austral winter when high values of precipitation, relative humidity and soil moisture, inhibit aeolian erosion. An important turning point towards a monotonic depletion of mass fluxes since 2014 is www.nature.com/scientificreports/ associated with a strong precipitation event on 2-8 April 2014 (Fig. 5c). Forte et al. 22 also reported a decrease on the frequency of remobilisation events since 2014 based on in-person interviews. This precipitation event most likely contributed to erode the primary deposit by surface water flow and also to increase the compaction due to increased soil moisture (Fig. 5b). Whilst the tephra-fallout deposit was produced and deposited within a few days during the eruption (4-15 June 2011), remobilisation processes continue up to the time of writing (2020). The temporal decay of mass flux is unequivocally not linear. A two-phase decay process is described by a double exponential fitting of the measured mass fluxes, with two strongly contrasting decay timescales (Eq. 2, Figs. 3c-e and 4). The significant mass flux peak due to the sudden large injection of volcanic ash in the system would be neglected by using only a single exponential fit. Indeed, this behaviour is of particular importance in volcanic regions since they are characterised by large amounts of material entering "instantaneously" into the erosional system once a deposit is emplaced. As a consequence, there is a sharp distinction between the first phase, clearly associated with the input of fresh volcanic ash, and the second phase, where deposit stabilisation processes may be developing. Despite the large uncertainty associated with the timescales (see parameters in the Fig. 3c-e), we find that the transition point between the two-phases based on measurement of trapped sediment at 0.15 m and 0.50 m above ground level occurs approximately 5 to 7 months after the eruption onset, whilst this transition occurs 8 months after the eruption for the 1.50 m sampling height. These periods coincide with the end of austral spring and the beginning of austral summer (i.e., October-January), characterised by strong winds and low soil moisture, that might have contributed to a rapid depletion of material due to intense remobilisation. Interestingly, a similar period of about 4-6 months of intense remobilisation was reported by Wilson et al. 14 after the eruption of the Hudson volcano in August 1991. This faster shift of phases at levels below 0.50 m may correspond to the height of the saltation layer where major particle transport occurs 41 . Additionally, we find that, for the total mass at the three sampling heights, the two timescales are orders of magnitude different: a faster τ 1 of 59 ± 16, 55 ± 16, 77 ± 16 days, and a slower τ 2 of 775 ± 342, 516 ± 186, 445 ± 175 days, at 0.15 m, 0.50 m and 1.50 m, respectively. We noticed that, whilst a slower decay is associated with the 1.50 m in the first phase, it is opposite in the second phase (slower at 0.15 m). This trend, as a function of height for the total mass (Fig. 3), is also confirmed by individual particle size classes (Fig. 4). Although we cannot identify a dominant transport mechanism (such as saltation or suspension 38 ), these results suggest that a large fraction of the particles in transport occurred at low levels (< 0.50 m) during the first phase ( ∽June 2011-January 2012), while a larger fraction of the particles in transport shifted to a higher level (> 0.50 m) during the second phase ( ∽February 2012-December 2016).
Mass flux decay timescales are strongly dependent on particle size. Figure 4 shows that the mass fluxes of particles of different sizes decay at different timescales. We expect that decay timescales are related to the particle threshold friction velocity u * t since this physical quantity is the minimum friction velocity at which particle motion starts. It has previously been demonstrated that u * t as a function of particle size is described by an U-shaped curve with a minimum for particle sizes around 75-100 µm 55 , where inter-particle forces are approximately equal to the gravitational force 42 . For particles smaller than approximately 75 µ m, u * t rapidly increases as particle size decreases due to the increasing strength of cohesive forces; in contrast, for particles greater than approximately 100 µ m, u * t increases due to their larger weight 44,55,56 (ESM- Fig. S2). Interestingly, we find the same U-shaped trend between the timescale of decay in mass flux and particle size suggesting that this methodology may be extremely useful in linking the saltation theory of particles with field measurements in natural environments. The complexity of natural environments, however, requires a careful interpretation. The contribution of volcanic particles dominant during the first phase is mainly controlled by the size classes present in the primary volcanic deposit, whilst the second phase is likely controlled by a more general trend of erosion in this region with a mixture of the CC volcanic particles, old sediments and soil. In particular, particles of 19 to 37 µ m diameter, corresponding to the most abundant classes in volume in the primary deposit (ESM -Table S4), show a faster mass flux decay in both phases, although this trend is much more clear for the phase II. Indeed, τ 2 as a function of particle size displays a distinct U-shaped trend, very similar to the theoretical u * t trend, but shifted towards a minimum for a particle size between 19 and 37 µ m instead of 75-100 µ m particles. This discrepancy could be due to the fact that theoretical u * t curve is based on experiments of individual particle sizes and, thus, it does not account for the interaction between particles in a multi-sized distribution, as it is the case of volcanic deposits. Additionally, the lack of reliable experimental data for particles smaller than 50 µ m 41 means that the theoretical value of u * t is poorly constrained for this size range. Given the high complexity, the effect of grainsize wind segregation (resulting in good sorting) and the interaction between saltation and suspension-induced processes (i.e., splashing generating short or long-term suspension of fine particles 41,42 ), are not considered in the theoretical u * t . In other words, the timescales reflected from field measurements are related to a dynamic threshold friction velocity, whilst the theoretical u * t is estimated for particles initially at rest (static threshold friction velocity). Indeed, the distinct trend of τ 2 with 2 minima values at 2 µ m and 19-37 µ m might be associated with the effects of complex saltation and suspension-induced transport mechanisms in a supply-limited system where the particles with shorter decay timescales (i.e., low τ ) can be remobilised more easily until the total depletion. Therefore, a crucial parameter in aeolian processes of volcanic deposits is the removal rate of particles. The comparison of Q th with Q msd clearly shows that measured mass fluxes become much lower than the theoretical estimates with time (Fig. 6a). Since the theoretical saltation models assume an infinite supply of particles, the variations in Q th are only associated with meteorological fluctuations (e.g., seasonality). Even though parameters that also modulate particle flux have not been considered in this study due to the complexity of their numerical analysis (e.g., vegetation, surface crust development, gusts or snow), their influence likely follow seasonal variations as observed by Shinoda et al. 57 in temperate grasslands. However, it has been demonstrated that, in this study, Q msd decays exponentially with time (Figs. 3, 4, 5c and 6a). These results suggest that, although meteorological conditions play an important role in aeolian processes, the exponential variation of mass transport Scientific RepoRtS | (2020) 10:14456 | https://doi.org/10.1038/s41598-020-71022-w www.nature.com/scientificreports/ of volcanic particles with time is substantially governed by the amount of available material in the long-term (i.e., months to years), as it has also been observed in supply-limited studies for loess deposits 47 , numerical simulations 48 , and experiments 53 . The limitation on material supply might be mainly related to the compaction of the primary volcanic deposit, the assimilation of ash by vegetation, and the loss of fine material with time, all processes being very difficult to numerically quantify. One possibility to account for compaction is through the cohesive forces described by the parameter γ of the Shao and Lu model 44 . However, this scheme assumes that cohesive forces are linearly proportional to the particle diameter, but the effect of moisture and chemical bindings are ignored. An alternative approach to consider the effect of soil compaction is the soil moisture correction given by Fecan et al. 45 parameterisation. However, this has not been possible to explore since current soil reanalysis models (e.g., ERA-Interim Reanalysis, from which our dataset comes) do not account for a real and updated volumetric soil moisture once a volcanic deposit is emplaced. Multiple dust emission schemes have been successfully applied to model the removal of fine ash particles suspended in the atmosphere 17 . However, modelling resuspension of fine particles for such a long time series (2011 to 2016) and a widespread area ( ∼ 15 km 2 ) is currently not possible due to computational constraint. Notwithstanding, even when we cannot quantify the compaction of the parental deposit and the loss of fine particles, we propose a simple and reproducible model that allows for the reconstruction of the mass flux temporal decay accounting for all the key processes involved in the depletion of material in a real case. This empirical model is based on the temporal variation of the erodibility σ (t) , which decreases with time from an initial erodibility σ 0 (in this case, 0.21) with a rate proportional to the integral of the mass fluxes in the preceding period [t 0 − t] (i.e., t t 0 Q slm dt , where t 0 is a reference time, Eqs. 3 and 4). The constant of proportionality α can be then empirically determined (in this case, its optimal value is 1.3 x 10 −2 m kg −1 , Fig. 6b). This model provides a very good agreement between the measured Q msd and the theoretical Q slm streamwise mass fluxes, considering a supply-limited system (Fig. 6a). In addition, the model is able to reproduce the mass flux peaks for specific sampling periods (i.e., August 2011, December 2011, January 2013 and September 2013, Fig. 6a grey curve), which are not recognizable when considering only meteorological fluctuations ( Fig. 6a green curve).
This study provides significant insights into the temporal variation of remobilised ash, characterised by a twophase process that can be empirically modelled. We can expect the most impactful first phase (associated with high mass fluxes of remobilised volcanic ash) to occur within a few months after the eruption, and a long-lasting second phase with exponentially decreasing mass fluxes to occur up to years and decades later. These results have important implications for risk reduction of active volcanic regions. Mitigation measures at short, medium and long term should be considered, particularly for regions that could be potentially affected by significant deposition of fine ash and strong winds, and characterised by large flat areas not protected by topography or vegetation, as it is the case of the Argentinian Patagonia.

Conclusions
Aeolian erosion of volcanic particles is a stochastic, complex and extremely dynamic process, whose timescales can vary from days to millennia after eruptions. This study provides a comprehensive analysis of supply-limited processes which are very common in natural environments (e.g., loess deposits, sand dunes, volcanic deposits). Understanding these processes and their effect over the long-term on the environment and the communities exposed to volcanic eruptions is extremely important. It is difficult to constrain the interactions between all the involved parameters and accurately determine which are the most relevant for controlling the mass flux decay over time. Our results suggest that compaction processes and the loss of fine material with time, which ultimately limits the supply of particles, play dominant roles on the timescales of aeolian transport of fresh volcanic deposits. Since these processes are difficult to quantify numerically, we propose an exponential model based on the temporal evolution of the erodibility, the factor accounting for the potential of the surface to be eroded, and a proportionality constant (Eqs. 3 and 4). Such a simple and reproducible model for supply-limited systems has never been developed in mineral or volcanic erosion studies. This proposed empirical model is able to describe the depletion of particles due to aeolian activity, under the conditions of the studied volcanic area. Further investigations are needed to test the suitability of this model under different conditions. We have also shown that, in the case of volcanic eruptions where a large volume of material is suddenly injected into the aeolian-transport system, the timescales of mass flux decay, as a function of particle size, follow the characteristic U-shaped curve; such a trend relates the threshold friction velocity to particle size. The proposed methodology represents an innovative approach to link field measurements with saltation theory. However, further studies of supply-limited volcanic systems are necessary to better understand which parameters contribute the most to remobilisation of volcanic ash. Finally, it is important to stress that the field measurements and techniques necessary to characterise ash-remobilisation deposits and processes are not the same as those used to investigate primary pyroclastic deposits. To address this, a detailed uncertainty analysis is presented in the Methods and Materials section. Future work should focus on the development and implementation of field, experimental and modelling strategies to assess and account for all the relevant parameters (e.g., soil moisture, soil compaction and cohesive forces, soil bulk density, precipitation rate thresholds, surface roughness).

Materials and methods
Chronology of the 2011-2012 CC eruption. After 20,49,58 . The eruption produced plumes with heights up to 12 km above the vent and generated about 1 km 3 of rhyolitic tephra, mostly during the first 15 days of eruption, that was dispersed towards the Argentinian Andes and the Patagonian steppe by prevailing westerly winds. Long-lasting eruptive dynamics combined with strong and shifting winds produced a complex primary tephra-fallout sequence of Scientific RepoRtS | (2020) 10:14456 | https://doi.org/10.1038/s41598-020-71022-w www.nature.com/scientificreports/ deposits over a wide area [50][51][52] . The main primary tephra-fallout deposit consists of 4 main units subdivided into 13 layers 52 . Unit I, that corresponds to the climactic phase (4 and 5 June), was characterised by plumes 9-12 km high above the vent, that dispersed to the SE and generated a lapilli-bearing bed set with a total volume of ca. 0.75 km 352 . Unit II corresponds to a variable wind direction period when plumes of 6-10 km above the vent, dispersed abruptly to the N and SE, producing a lapilli-bearing deposit of a total volume ca. 0.21 km 3 . Unit III was primarily produced and sedimented mainly from the 6 to 15 June and was characterised by fluctuating plume heights that dispersed to the E-ESE, and generated a fine-grained deposit consisting of 5 layers (K1 to K5) 50,52 . The volume of the K2 layer was estimated to be ca. 0.05 km 352 . Unit IV was deposited after the 15th June, and is composed of a millimetres-thick fine ash deposit and is found only in very proximal areas (< 20 km from the vent). On 20 June 2011, an effusive phase started with the emission of viscous lava that lapsed until mid-2012 although its emplacement continued until January 2013 59 . Similar to the case of the 1991 Hudson eruption (Chile), one of the most important long-term consequences of this eruption has been the ash wind-remobilisation 14,20,22 , particularly in the arid and semi-arid regions of the Patagonian steppe. Given the associated stratigraphic position, the fine-grained material and the widespread sedimentation, Unit III represents the Unit most affected by remobilisation in the Argentinian region 23 .
Sampling strategy and data collection. In this study we used dedicated samplers (BSNE type 21 ) installed in 2011 by the National Agricultural Technology Institute of Argentina (INTA) as part of the National Soil Research Program, to systematically collect horizontal airborne material syn-deposition 60 . These collectors, located in 4 strategic erosion sites (Fig. 2), consist of a wind-oriented upright with a weather vane and 3 samplers at fixed heights, 0.15 m, 0.5 m and 1.5 m (Fig. 3a). Based on the same dataset, Panebianco et al. 21 thoroughly analysed the effects of topography, landscape geomorphology and vegetation cover on the total mass flux; in addition, Dominguez et al. 23 studied the physical features of particles (i.e., grainsize and shape) to infer information on the transport mechanisms involved. Here, we investigate the timescales of mass flux through a comprehensive analysis of the mass flux decay timescales as a function of particle size and height of sampling collection of the site S4, accounting for a total of 42 samples (14 collection periods at three different heights). Since the studied area is difficult to access, the periods between sample collections are not constant, and therefore the temporal resolution of our dataset varies from about 1-20 months (See details in the ESM- Table S3). Meteorological and soil parameters (i.e., direction and velocity of wind, total precipitation, volumetric soil moisture and roughness length) have been obtained from the ERA-Interim reanalysis data at surface (i.e., 10 m height) of the ECMWF dataset (European Center for Medium-range Weather Forecasts).
Grainsize analysis. Granulometry analyses were performed with a laser diffraction particle size analyser (CILAS 1180) at the University of Geneva. It is important to note here that these are disturbed samples due to the transport from the collector to the laboratory. We apply 60 seconds of ultrasound before measurement in order to disaggregate fine particles, but no chemical processes have been used. Grainsize distributions are expressed in half-phi intervals. The equivalence between phi and metric scales, as well as the representative size value for each interval, corresponding to the mid half-phi value, is presented in the ESM- Table S2.
Wind friction velocity. The ultimate force driving wind erosion depends on the transfer of momentum from the atmosphere to the surface, which is described by the momentum flux or the shear stress of the wind 38 . This is related to the wind friction velocity by, where τ s is the surface shear stress and ρ a is the air density (taken to be 1.112 kg m −3 ). Roughness elements on natural surfaces exert a drag on the flow due to the difference in pressure between the windward and the leeward sides of the element. Consequently, the transfer of momentum from the flow to the surface involves a complex process depending on the viscous and pressure forces on the roughness elements 41 . Near the surface, horizontal wind speed is commonly described by a log profile, under neutral stability conditions. It is accepted though, that above some height z, the mean wind profile U(z) is given by 61 , where κ = 0.4 is the von Karman constant, z 0 is the roughness length and depends on the geometric features of roughness elements (e.g., height, width, etc.), and z is the height above the surface. Following Stull 61 , this can be re-arranged into: We can obtain the forecast roughness length, z 0 , from the ERA-Interim dataset, taking z = 10 m.
Threshold friction velocity. Based on the models of Owen 39 and Iversen et al. 43 , Shao and Lu 44 derived a simple expression for the threshold friction velocity as a function of particle size, d, as follows, where f (Re * t ) is a function of particle Reynolds number at the threshold friction velocity (0.0123 from fitted experiments 44 ), ρ p is the particle density (in kg m −3 ), ρ a is the air density (1.112 kg m −3 ), g is the gravity constant (9.81 m s −2 ) and γ is a constant accounting for inter-particle cohesion. Shao and Lu 44 treat γ as a free parameter and fit the Eq. (9) to data obtained from wind-tunnel experiments 43 and found acceptable agreement for 1.65 × 10 −4 < γ < 5 × 10 −4 kg s −2 . However, there is not yet reliable experimental data to estimate γ for d < 50 µm 41 . According to the Eq. (9), u * t is dependent on the particle density, but the density of vesicular particles produced by volcanic eruptions can vary considerably as a function of grainsize 62,63 . Consequently, highly vesiculated large pumices have lower densities compared to small ash particles with lower vesicularity. Density studies suggest that for particles finer than 6 φ , the density is approximately constant and fragments can be considered as dense as the skeletal density, without vesicles (i.e., dense rock equivalent, DRE) 62,63 . In order to account for this density variation, here we linearly interpolate the density from 1 φ (i.e., 1270 kg m −3 for K2 fragments) to 6 φ (i.e., 2690 kg m −3 , DRE) using measurements from Pistolesi et al. 52 (See ESM- Table S4). The threshold friction velocity is strongly dependent on the soil moisture and vegetation cover. In this study, following Folch et al. 17 , we apply a correction factor accounting for soil moisture, f w (w) , according to the parameterisation of Fecan et al. 45 , as follows, where where w g is the gravimetric soil moisture, and w ′ is the maximum amount of water that can be adsorbed. The gravimetric soil moisture is related to the volumetric soil moisture w through w g = w ρ w /ρ b , with w being obtainable from the ERA-Interim dataset; and ρ w and ρ b being the water and soil bulk densities (here we use 998 kg m −3 and 900 kg m −351 , respectively). We used a parameterisation for w ′ depending on the clay content Cs (i.e., particles < 2µm), according to the expression, w ′ = 0.0014 Cs 2 + 0.17 Cs 45 . Since average clay content of the primary distal deposit is 4.13% 23 , w ′ is 0.7%.
Streamwise mass flux of particles. The streamwise (horizontal) mass flux Q is defined as the vertical integral of the saltation streamwise mass flux density q at each given height z [38][39][40][41][42] and can be expressed as, For field measurements, the average mass flux has been calculated according to, where q z is the time-averaged mass flux at a given height z (units of kg m −2 day −1 ), m z is the mass measurement at each height; A = 1 × 10 −3 m 2 is the surface area of the collector opening and dt is the number of days per period of collection. We fit a linear relationship of the form mz + n, the measured streamwise Q msd is then given by The theoretical streamwise mass flux, Q th , has been calculated according to the parameterisation of Shao et al. 40 based on the model of Bagnold-Owen 39 (Eq. 15). Additionally, we also imposed two meteorological conditions; first, the total precipitation rate p must not exceed a threshold p t = 0.01 mm h −1 , following Leadbetter et al. 30 , and second, the wind friction velocity u * must exceed the threshold friction velocity u * t , such that where C o is the dimensionless Owen Coefficient, typically close to 1 40 . Q th is expressed in kg m −1 s −1 .
if p < p t and u * ≥ u * t (d) www.nature.com/scientificreports/ This theoretical framework allows us to calculate Q th for uniformly sized particles; however, volcanic primary deposits have dispersed size distributions. Since the effects of both particle interactions and sorting of natural soils are very complex, we follow Shao et al. 46 and assume that Q th is not significantly altered by the presence of a range of particle sizes. Consequently, we calculate Q th for each particle size fraction d i , from the mass fraction, f i , of that particular size class in the parent soil (i.e., the distal tephra-fallout grainsize distribution). The total streamwise mass flux is given then by the sum of all size fractions, Fitting procedure for the double exponential decay. Mass flux decay has been analysed under the hypothesis of a double exponential dependence in time (Eq. 2). The uncertainty associated with the four fitted parameters (i.e., τ 1 , τ 2 , b, c) is determined by the analysis of the non-linear least squares method, imposing a confidence interval of 68% (i.e., 1 standard deviation) in the Matlab fitting toolbox. A comprehensive analysis for each size class has been done including the independent evaluation of each of the two terms of the double exponential combined with the error bars derived from this procedure. We highlight two particular cases: (i) particles larger than 180 µ m are practically not present in the primary deposit (see ESM- Table S4) but they are found in the samplers, either because they are coarser volcanic particles coming from the prevailing westerly winds with time (CC primary proximal deposit is coarser than distal one 23 ), or because they are old soil sediments. As a consequence, the mass flux decay of these size-classes is not consistent within this study and they are not shown in the Fig. 4. (ii) particle sizes of 26 to 53 µ m show strongly pronounced peak and a rapid decay that τ 1 is not able to describe. As a result, error bars associated with these classes are considerably larger than other classes, and, in the case of 1.50 m the fit does not converge for τ 1 . It is important to notice that these classes are the most abundant in volume in the primary deposit (ESM -Table S4), hence it is not surprising that the first mass flux peaks were quite strong. Uncertainty analysis. Aeolian erosion of volcanic particles is a very dynamic process. Although meteorological seasonality has been considered (based on the ECMWF ERA-Interim Re-Analysis dataset), several uncertainties, particularly associated with the volumetric soil moisture estimation, w, arose from the reanalysis model. On the other hand, due to lack of data, we have taken soil parameters to be static (i.e., initial values of primary deposit: soil bulk density, ρ b , and average clay content, C s ), neglecting changes due to seasonal variations, stochastic wind fluctuations and compaction effects. One of the major caveats in this study is the irregular and large time window of sample collection. The Patagonia steppe is a remote area of difficult access and it is not feasible to establish a more systematic sampling. For future studies, we strongly recommend both a temporal resolution of the sampling suitable for the possible mass flux decay and a regular monitoring of physical conditions of the parent soil surrounding the dust-collectors (e.g., soil density, moisture, thickness, grainsize). Table 1 summarises the main caveats associated with the analysed parameters. Given the complexity of aeolian dynamics due to the turbulence of gusty winds and the effect of roughness elements (e.g., vegetation), these factors have not been considered in this study.  Table 1. Caveats concerning the estimation of the parameters required for the calculation of the threshold friction velocity u * t and the streamwise mass flux Q th .

Parameter Abbreviation Units Caveats
Cohesive forces 44 γ kg s −2 This term considers cohesive forces linearly proportional to particle diameter, however, van der Waals, electrostatic, liquid and chemical forces, which are sensitive to soil properties (e.g., particle shape, surface texture, mineralogy, deposit sorting, bonding agents, such as moisture and soluble salts) are not considered 44 . All these parameters are very important in order to account for compaction of primary volcanic deposits.
Soil bulk density ρ b kg m −3 This is probably the most important parameter changing over time. It is strongly related to the infiltration of water and thus soil moisture. It could indirectly account for compaction of the deposit. In this study, it is considered constant and corresponds to the primary deposit (Unit III) density at the beginning of the eruption.
Volumetric soil moisture w % It is obtained from soil models that are not updated once the eruption occurs. From the ECMWF dataset, "volumetric soil moisture layer 1" corresponds to 1-7 cm depth. However, it is only the most superficial layer (few mm) and generally much drier than the one to be remobilised.
Max. amount of water that can be adsorbed w' % Only based on the clay content 45 , it does not account for infiltration or permeability of volcanic particles.
Precipitation rate threshold p t mm h −1 The precipitation threshold inhibiting remobilisation processes is not well constrained. Here we use 0.01 mm h −1 , following the UK Met Office modelling strategy 30 .
Particle density ρ p kg m −3 Particle density is strongly dependent on the grainsize and vesicularity, which is particularly important in volcanic particles. It is difficult to measure density for particles smaller than 1φ . There are linear 62 and sigmoidal models 63 to fit density distributions. Here we use a linear model from 1 φ particles and DRE measurements.
Roughness length z 0 m The abrupt change in surface roughness is not considered by the atmospheric forecast models. However, it is strongly dependent on the large amount of loose particles on the surface due to fresh volcanic deposits. It has important implications on the wind friction velocity. Here we use the forecast surface roughness of the ERA-Interim Re-Analysis dataset.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creat iveco mmons .org/licen ses/by/4.0/.