Modeling volcanic ash resuspension – application to the 14–18 October 2011 outbreak episode in central Patagonia, Argentina

Volcanic fallout deposits from the June 2011 Cordón Caulle eruption on central Patagonia were remobilized in several occasions months after their emplacement. In particular, during 14–18 October 2011, an intense outbreak episode generated widespread volcanic clouds that were dispersed across Argentina, causing multiple impacts in the environment, affecting the air quality and disrupting airports. Fine ash particles in volcanic fallout deposits can be resuspended under favorable meteorological conditions, particularly during strong wind episodes in arid environments with low soil moisture and poor vegetation coverage. As opposed to eruption-formed ash clouds, modeling of resuspensionformed ash clouds has received little attention. In consequence, there are no emission schemes specially developed and calibrated for resuspended volcanic ash, and few operational products exists to model and forecast the formation and dispersal of resuspension ash clouds. Here we implement three dust emission schemes of increasing complexity in the FALL3D tephra dispersal model and use the 14– 18 October 2011 outbreak episode as a model test case. We calibrate the emission schemes and validate the results of the coupled WRF–ARW (Weather Research and Forecasting – Advanced Research WRF)/FALL3D modeling system using satellite imagery and measurements of visibility (a quantity related to total suspended particle concentration at the surface) and particulate matter (PM 10) concentration at several meteorological and air quality stations located at Argentina and Uruguay. Our final goal is to test the capability of the modeling system to become, in the near future, an operational forecast product for volcanic ash resuspension events.


Introduction
Resuspension and dispersal of volcanic ash by wind is of concern to human health and the environment (Baxter, 1999;Wilson et al., 2012). Micron-size ash particles suspended at low atmospheric levels deteriorate the quality of air and cause severe impacts on the local population and animals, causing irritation of mucosae, chronic respiratory symptoms resulting from inhalation of ash, and dispersion of toxic chemicals (see Baxter, 1999). Other hazards include disruption to ground transportation systems and airports (see Guffanti et al., 2009). All these impacts can occur far away from the original deposit region because, like dust clouds, resuspended ash clouds can be dispersed over large distances. The lift of ash and subsequent formation of clouds is enhanced under particular meteorological (such as strong winds and surface friction velocities) and ambient (such as low soil moisture, no vegetation) conditions, and can occur in both fresh and relic ash fallout deposits. For example, favorable meteorological conditions occurred during 20-21 September 2003 resulted in continuous resuspension of relic volcanic ash from the Katmai volcano and the formation of ash clouds that were transported up to 230 km into the Gulf of Alaska, affecting operations at the Kodiak Airport (Hadley et al., 2004). Other examples include resuspension of ash from  (Hobbs et al., 1983), Soufriere Hills in Montserrat (Hincks et al., 2006), Mount Hudson in the southern Argentinean Patagonia , Eyjafjallajökull in Iceland (Thorsteinsson et al., 2012;Leadbetter et al., 2012), and Cordón Caulle in the central Argentinean Patagonia, which is addressed in this paper.
Modeling of atmospheric dispersal and sedimentation of ash clouds from volcanic eruptions has been an active topic of research during the last two decades, and gained heightened interest in the aftermath of the major civil aviation disruptions following the 2010 Eyjafjallajökull and the 2011 Cordón Caulle eruptions. Several Volcanic Ash Transport and Dispersion Models (VATDM) (e.g.,  run operationally at the Volcanic Ash Advisory Centers (VAACs) and other institutions for forecasting purposes, and the underlaying modeling strategies are being substantially reviewed and improved (Bonadonna et al., 2012). In contrast, modeling of resuspended ash clouds has received little attention. This is surprising because, even if potentially less hazardous, resuspension-formed clouds can also trigger substantial impacts (see section 4). In fact, modeling of resuspended ash has been recognized as a research priority during a recent joint WMO (World Meteorological Organization)-VAAC modeling workshop (NOAA, Washington DC, 5-9 November 2012). The examples in literature regarding modeling of resuspended ash clouds are very scarce. Barsotti et al. (2010) proposed a simple fit based on data from Hincks et al. (2006) to estimate resuspended PM 10 surface concentrations at various populated locations around the Etna volcano. More recently, Leadbetter et al. (2012) modeled resuspension events form the 2010 Eyjafjallajökull deposits using the UK Met Office Lagrangian dispersion model NAME-III (Numerical Atmospheric dispersion Modeling Environment). Given the similarities between volcanic ash and mineral dust, Leadbetter et al. (2012) considered a simple dust emission scheme to compute the mass flux of resuspended ash depending on the wind friction velocity and the occurrence of precipitation. This has been used by the London VAAC/UK Met Office to provide an operational warning forecast service for resuspended ash to IMO (Icelandic Met Office).
The objective of this paper is to test different dust emission schemes of increasing complexity for implementation in volcanic ash transport models. To this purpose, we implement three different schemes in the FALL3D model (Costa et al., 2006;Folch et al., 2009) and consider the Cordón Caulle 14-18 October 2011 outbreak episode as a test case. We use a meteorological station located in Buenos Aires (1380 km from the volcano) to calibrate the emission of ash. This is necessary because uncertainties exist in both source strength parameters and formulation, including that the emission schemes have been originally developed for mineral dust rather than for volcanic ash, and poor constraints on the properties and grain size distribution of ash particles along the deposit. After calibrating, we compare the model results with satellite retrievals and particulate matter surface concentrations inferred from observations at different meteorological and air quality stations across Argentina. The final goal is to test the modeling strategy before its implementation as an operational forecast product at the Argentinean National Meteorological Service (SMN) and Buenos Aires VAAC.

Previous considerations on modeling volcanic ash resuspension
Modeling of volcanic ash resuspension requires the implementation of emission schemes in VATDMs. Emission schemes give the mass flux of resuspended (ash) particles depending on meteorological conditions, soil moisture, terrain roughness, and characteristics of the fallout deposit (size and density of particles, grain-size distribution, deposit thickness). Typically, soil moisture is obtained from the Numerical Weather Prediction (NWP) model driving the VATDM (WRF-ARW in our case). This can introduce an unknown uncertainty as NWP models do not include updated information about recent fallout deposits, which can substantially alter moisture and roughness of the original soil. On the other hand, fallout deposits are heterogeneous, introducing a second source of uncertainty regarding granulometric properties. Heterogeneities in grain-size distribution and particle properties exist because the variation of settling velocity with particle size and density (causing strong granulometric variations along the dispersal axis), occurrence of different eruptive phases (e.g., variable eruption rates, column heights, degree of magma fragmentation, etc.), or transport under unsteady heterogeneous wind fields. Moreover, even if a fallout deposit is well characterized through a dedicated field campaign, information is rarely available in a regular grid, but rather across transects or along the most accessible deposit parts (e.g., along roads). It follows that, from a modeling perspective and in practical terms, an option to obtain the granulometric characteristics at all deposit grid points is to run a preliminary simulation for the fallout using the total grain size distribution, typically reconstructed from field data.

Dust emission schemes
Saltation impact represents the most effective mechanism for resuspension of smaller-size particles in soils (Shao et al., 1993). When the intensity of wind blowing across a granular soil exceeds a certain threshold, grain particles begin to saltate. Experiments with sand-sized particles show that the impact of saltating mid-size grains (larger than about 50 µm) when falling back to ground breaks the cohesive forces of smaller particles, enhancing their suspension. For this reason, the emission rate (vertical flux of particles), defined as the mass emitted per unit of area and time, strongly depends on the horizontal (saltation) flux of larger particles. In recent years, various emission schemes for mineral dust have been proposed and implemented in atmospheric transport models to simulate long-range transport and deposition of mineral dust. This section summarizes different dust emission schemes that we will test for volcanic ash. We follow an approach similar to that of Darmenova et al. (2009) and Kang et al. (2012), who tested different dust emission schemes in the WRF and WRF-Chem models respectively.

Threshold friction velocity
The friction velocity u * is a reference velocity used in surface boundary layer theory to scale the shear stress of a fluid, assumed proportional to the square of the mean velocity. The threshold friction velocity, defined as the wind friction velocity at which soil erosion initiates (Greeley and Iversen , 1985), depends on physical properties of soil particles (size and density) and on surface conditions such as soil moisture and roughness. Soil moisture reinforces soil particle cohesion and therefore inhibits erosion and resuspension. The presence of rough elements on the ground (e.g., non erodible rocks or vegetation) also decreases the emission rate of particles because the elements act as wind shelters and absorb part of the momentum of the wind, leading to a decrease of wind shear stress acting on the erodible surface. A precise quantification of the friction velocity is challenging and requires data on soil properties and meteorological conditions, which may be unavailable on a local scale. Simple dust emission schemes assume a constant threshold friction velocity over the region of interest regardless of the particle size. For example, Leadbetter et al. (2012) modeled the resuspension of ash after the 2010 Eyjafjallajökull eruption in Iceland assuming a threshold friction velocity of 0.4 m s −1 for all particle diameters in the range 1-10 µm combined with a precipitation rate cutoff of 0.01 mm hr −1 . The cutoff accounts for a critical soil moisture above which the emission of particles is assumed to be inhibited. However, the constant friction velocity approach can be inadequate over large areas because particle properties and soil characteristics are likely to vary substantially in space. For example, in the particular case of volcanic deposits, a decrease in the mean particle size is expected at more distal locations. In order to account for heterogeneous particle and soil characteristics, the threshold friction velocity u * t can be expressed at each point as (Shao, 2001) where u * ts (d) is the threshold friction velocity on a bare dry surface for particles of size d (for irregular particles d is assumed to be the equivalent particle diameter), and f w (w) and f λ (λ) are correction functions for soil moisture w and surface roughness λ, respectively (f w (w) ≥ 1 and f λ (λ) ≥ 1). A number of experimental and theoretical studies have determined parameterizations for threshold friction velocities on bare dry surfaces. Marticorena and Bergametti (1995) fitted experimental data obtained in wind tunnels by Iversen and White (1982) for various particle densities (from 210 to 1135 kg m −3 ) and diameters (from 12 to 1290 µm), and derived the following expression for u * ts depending on particle size and density: u * ts = 0.129K (1.928Re 0.092 −1) 0.5 0.03 < Re ≤ 10 0.129K(1 − 0.0858e −0.0617(Re−10) ) Re > 10 (2) with K = ρ p gd ρ a 1 + 0.006 ρ p gd 2.5 and Re = 1331 × d 1.56 (the lower bound of the fit corresponds to particles of ≈ 10 µm in size). In the expressions above, ρ p and ρ a are particle and air densities (expressed in g cm −3 ), g is gravity (in cm s −2 ), d is the particle size (in cm), Re is the Reynolds number parameterized as a function of the particle size, and u * ts is given in cm s −1 . In turn, Shao and Lu (2000) derived an expression for u * ts considering spherical particles with a cohesion force proportional to particle size: where γ is a parameter ranging between 1.65 × 10 −4 and 5 × 10 −4 kg s −2 (a value of 3×10 −4 kg s −2 is assumed here). Using experimental data from literature, Fecan et al. (1999) parameterized the increase in the threshold friction velocity due to soil moisture in arid and semi-arid regions and derived the following expression: where w is the maximum amount of water that can be absorbed (depending on soil texture) and w g is the gravimetric soil moisture (water to soil bulk mass ratio), w g = wρ w /ρ b with ρ w and ρ b being the water and soil bulk densities respectively, and w is the volumetric soil moisture in % (water to soil bulk volume ratio). The values of w can be obtained from a NWP model. When the gravimetric soil moisture content w g is close to but smaller than w (dry soil), interparticle capillary forces are not strong enough to increase the erosion threshold significantly. The value of w ranges from 0 % for sand to ≈ 30 % for pure clay (Fecan et al., 1999). It is difficult to give a value of w for volcanic fallout deposits because the amount of water that can be absorbed strongly depends on porosity which, in turn, depends on variable factors controlling tephra formation (e.g., magma rehology, degree of fragmentation, amount of volatiles in magma, etc.). In this study we assume that w = 10 %. Figure 1 shows the dependency of the threshold friction velocity u * t on particle size with and without soil moisture correction. Note that particles in the range of 30-200 µm are more likely suspended. On the other hand, several parameterizations have been proposed for the so-called drag partition coefficient, the inverse of f λ (λ) (e.g., Marticorena and Bergametti, 1995), depending on the roughness length over a smooth surface and

Westphal Marticorena Shao
w =00% w =25% 1: Dependency of the threshold friction velocity u * t on particle size (in µm) according to Shao and Lu (2000) (red line, ) and Marticorena and Bergametti (1995) (blue line, eq. 2). Solid lines are without moisture correction (w = 0). Dotted show the effect of moisture using eq. 4 with w = 25% and w = 10%. For comparison, the solid black line shows the tant threshold friction velocity u * ts = 0.25 m/s used in WE scheme and the dotted black line its correction when a ture of w = 25% is considered. Fig. 1. Dependency of the threshold friction velocity u * t on particle size (in µm) according to Shao and Lu (2000) (red line, Eq. 3) and Marticorena and Bergametti (1995) (blue line, Eq. 2). Solid lines are without moisture correction (w = 0). Dotted lines show the effect of moisture using Eq. (4) with w = 25 % and w = 10 %. For comparison, the solid black line shows the constant threshold friction velocity, u * ts = 0.25 m s −1 , used in WE scheme and the dotted black line its correction when a moisture of w = 25 % is considered. the aeolian roughness length. However, it is very difficult to assess these parameters in the case of ash fallout deposits. An obvious difficulty is that the deposition of ash modifies the underlying surface (to an extent that depends on the thickness of deposit) affecting any information embedded in the land modules of NWP models (e.g., roughness length). For this reason the influence of terrain roughness on the threshold friction velocity will not be considered in this study.

Horizontal (saltation) flux
The horizontal flux of saltating particles, i.e., the stream-wise flux of saltating particles integrated along the vertical, measures the "intensity" of saltation, which strongly affects the emission rate of smaller-size particles. Following the theory of saltation and experimental results from Owen (1964), Shao et al. (1993) proposed the following parameterization: where F H is the horizontal (saltation) flux (units of kg m −1 s −1 ) of saltating particles of size d s , and c o is an empirical dimensionless constant close to 1. A similar expression was proposed by Marticorena and Bergametti (1995) and Marticorena et al. (1997) after the seminal work of White (1979): with c w = 2.61 according to the original wind tunnel experiments (White, 1979) and c w ≈ 1 according to successive corrections (Marticorena et al., 1997).

Vertical flux (emission rate)
The simplest dust emission parameterizations depend only on meteorological conditions (typically on a power of the friction velocity). For example, Westphal et al. (1987) measured the vertical flux of aerosol dust particles (<10 µm radius) from sandy, loamy, and clay soils depending on the friction velocity and obtained the following least squares fit: where F V is the vertical flux (in kg m −2 s −1 ), occurring only above a (constant) threshold friction velocity of u * t = 0.3 m s −1 in the particular experiment of Westphal et al. (1987). An important limitation of Eq. (7) is that the vertical flux does not depend on particle size or soil moisture. Although very simplistic, this parameterization can be useful when information on soil characteristics (e.g., particle sizes and densities, moisture, roughness, etc.) is poorly constrained or unavailable. A slightly more sophisticated approach consists of modeling the emission rate as a function of the difference between the friction velocity and the threshold friction velocity. This embeds all the information on soil properties within the threshold friction velocity and allows one to compute fluxes depending on particle size. For example, from the parameterization in Marticorena et al. (1997b): where K is a soil texture coefficient equal to K = 5.4 × 10 −4 m −1 from the experiments of Gillette et al. (1997). This parameterization has been used, among others, by Draxler et al. (2001) to simulate PM 10 concentrations from dust storms using the HYSPLIT model. Finally, more sophisticated emission schemes consider particle-particle interaction by saltation bombardment and, in some cases, by aggregate disintegration (e.g., Shao, 2001). Shao et al. (1993) considered that the vertical flux (emission rate) of particles of size d caused by the saltation bombardment of particles of size d s (d s ≥ d) is proportional to the horizontal flux of saltating particles: where α (units of m s −2 ) is a coefficient of blasting efficiency determined experimentally (Shao and Leslie, 1997;Shao, 2001):  with d and d s expressed in mm. Given the particle grain size distribution p(d) at each point of the deposit and assuming a discretization in a series of bins ( p(d bin ) = 1), the emission rate of particles of size d can be obtained by summing up the contribution from saltating particles of all sizes equal or larger than d: where d max is the maximum considered size (typically few hundreds of µm). Note that, when u * u * t , the vertical flux F V (d) becomes proportional to u 3 * in both parameterizations Eqs. (8) and (11) (with F H given by either Eqs. 5 or 6).

Tested emission schemes
In this study, we test three different emission schemes (of increasing sophistication) for the case of volcanic ash.
-Emission scheme 1, hereafter referred as the WE scheme. Consists of computing the emission rate using Eq. (7) and a threshold friction velocity u * t independent of particle size (Westphal et al., 1987).
-Emission scheme 2, hereafter referred as the MB scheme. Consists of computing the emission rate using Eq. (8) and the threshold friction velocity u * t (d) using Eq.
Given that these schemes have been developed theoretically and calibrated experimentally for mineral dust rather than for volcanic ash, we multiply the emission rates by a correction factor φ that will be determined by comparing model values with observations, in our case form the Cordón Caulle 14-18 October 2011 outbreak episode.

The 14-18 October 2011 resuspension event
The June 2011 eruption from Puyehue-Cordón Caulle Volcanic Complex (CCVC) in Chile blanketed with volcanic ash a vast area of the Argentinean Central Patagonia (Collini et al., 2013). Recurrent ash resuspension events occurred during the following months due to the strong (≈ 100 km/h) westerly and southwesterly winds that typically blow over Patagonia during the austral spring (Lässig et al., 1999;Peri et al., 2002). The small villages spread sparsely across the Patagonian steppe (e.g., Ing. Jacobacci in the Río Negro province, see Fig. 2) were heavily impacted by this phenomenon, which hindered ordinary activities and forced inhabitants to remain indoors during the strongest wind episodes (Wilson et al., 2012). By mid-October 2011, a particularly strong resuspension episode occurred during the passage of a southwestern frontal system crossing northern Patagonia with surface wind speeds of 65-85 km h −1 and maximum gusts of about 95 km h −1 measured by Bariloche and Neuquén meteorological stations. As a result, a widespread ash cloud reaching the 850 hPa atmospheric level (1.5 km elevation roughly) was formed and dispersed rapidly east-northeast across Argentina (Damiani, 2011). Impacts occurred at a national level. Main routes in northern Patagonia (e.g., the famous Route 40 linking Bariloche city with the Neuquén province, see Fig. 2) were closed during 15-16 October due to very low visibility and after the occurrence of accidents. This disruption affected the transportation of persons and goods in the southern part of the country. During and after the afternoon of 15 October, the ash cloud reached the provinces of Río Negro, La Pampa, and west of Buenos Aires. In the early morning of 16 October (around 10:00 UTC), a dense ash cloud was clearly visible over the sky of the metropolitan area of Buenos Aires (Ciudad Autónoma de Buenos Aires, CABA), at nearly 1375 km from CCVC, and a perceptible ash layer was deposited at ground level, blanketing the city, as reported by the meteorological station METARs. The air quality stations of the Government of CABA (GCBA) registered the ash cloud from its arrival and measured a dailyaveraged PM 10 level of 252 µg m −3 . This largely exceeds the US EPA NAAQS (Environmental Protection Agency National Ambient Air Quality Standard) as well as the national legislation standard (law 1356, Government of Buenos Aires) 124 Folch et al.: Modeling volcanic ash resuspension limit of 150 µg m −3 PM 10 concentration for 24 h exposure. For the first time since the June 2011 CCVC eruption, the national authorities issued warnings regarding transportation in the Buenos Aires province due to the low visibility owing to the presence of suspended ash at low atmospheric levels. Around midday, flights from the Aeroparque Jorge Newbery and Ezeiza international airports, located in CABA, were suspended. Up to 146 flights were canceled on Sunday 16 October alone. The morning after, some scheduled flights began to depart after the cleaning of platforms and runaways but, nevertheless, the main commercial airlines did not resume their operations until Monday 17 October in the afternoon. During this day, the ash cloud covered the south of Uruguay disrupting the Montevideo international airport (40 canceled flights). Flight disruptions affected also the Argentinean airports located in the cities of Córdoba, Mendoza and San Luis until 18 October afternoon because a branch of the ash cloud moved towards these inner provinces making the sky almost invisible.

Modeling strategy
We modeled the 14-18 October 2011 Cordón Caulle resuspension event using the FALL3D dispersal model (Costa et al., 2006;Folch et al., 2009) coupled offline with the Weather Research and Forecasting (WRF-ARW) meteorological model (Skamarock et al., 2008). FALL3D uses 4-D meteorological fields generated offline and volcanological inputs to produce time-dependent variables like airborne ash concentration, ash cloud column load or ground deposit thickness. The new version of the code, FALL3D-7.0 , includes new capabilities to resolve a "continuum" of sources over a region, as opposed to a single point source or a set of vertical point sources typically considered for eruption columns. As already mentioned in Sect. 2, our modeling strategy builds upon a preliminary eruption simulation run to characterize the deposit (extent, thickness and granulometry at each point) followed by the simulations of resuspension.

Deposit characterization
As a starting point, we ran the WRF-ARW/FALL3D modeling system for the period 4-20 June 2011, during which most deposition of tephra from Cordón Caulle eruption occurred. We considered a Gaussian total grain size distribution (TGSD) discretized in 10 bins ranging from −1 (2 mm) to 8 (4 µm) and a linear dependency of particle density with diameter (end-member density values of 1000 and 2200 kg m −3 ). Note that the largest particles cannot be transported substantially by resuspension but, nonetheless, we considered them because they play a key role in the SH emission scheme (see Eq. 11). FALL3D-7.0 can handle ash aggregation phenomena but not particle disaggregation, Complex (CCVC) is indicated by a triangle.  which typically occurs when (fragile) aggregates impact the ground. For this reason, and because our primary interest is to model subsequent resuspension, particle aggregation was not considered for the deposit simulation. This introduces some uncertainty because the real amount of fine particles in the proximal deposit (originally fallen as aggregates and then disaggregated) can be larger than that predicted by the model. Eruption column heights in the model oscillate daily from around 10 km a.s.l. (above sea level) down to less than 3 km a.s.l. by 20 June 2011 (Collini et al., 2013). A comparative study between modeled deposit and field observations can be found in Osores et al. (2012). Fig. 3 shows the modeled deposit that we consider as the potential source area for resuspension. Note that, as observed in the field (see Collini et al., 2013), the main deposition lobe is directed southeast.

Meteorological driver and emission schemes
We used the modeled deposit to initialize the resuspension simulations using a WRF-ARW run from 14-18 October 2011 as the meteorological driver. This WRF-ARW run uses 0.5 • Global Forecast System (GFS) analysis and forecasts supplied by the National Centers for Environmental Prediction (NCEP) of the National Oceanic and Atmospheric Administration (NOAA) as initial and boundary conditions. The WRF-ARW model was configured with a horizontal resolution of 12 km, 38 vertical levels, the Ferrier scheme for microphysics, the Rapid Radiative Transfer Model (RRTM) for   long-wave radiation, the Dudhia model for shortwave radiation, the NOAH land surface model, the Betts-Miller-Janjic scheme for convection, and the Mellor-Yamada-Janjic for the planetary boundary layer. The computational domain on the inner WRF-ARW nest spanned between 21 • -48 • S; 30 • -90 • W. For illustrative purposes, Fig. 4 compares the WRF-ARW predicted wind velocities (first model layer, 10 m above ground) with measurements at three different stations located at Bariloche (M1, 108 km from CCVC), Neuquén (M2, 386 km from CCVC) and Buenos Aires (Aeroparque airport, M13, 1379 km from CCVC). Note how two different velocity peaks, responsible for two resuspension events, are clearly visible at M1 and M2. The predictions of WRF-ARW at station M13 on CABA show the largest discrepancy during 16 October, where differences of up to 4 m s −1 exist, concurrently with the plume arrival at CABA. In our opinion, this discrepancy might be due to the Río de la Plata breeze effects not well resolved by the NWP model. Figure 5 shows the emission rate integrated over the deposit as predicted by the three emission schemes over the 14-18 October 2011 period. Two major resuspension events are evident, the first starting on 14 October at around 18:00 UTC and the second, more intense, starting 24 h later. Note in Fig. 5a how these events coincide with values of u * 126 Folch et al.: Modeling volcanic ash resuspension Folch et al.: Modeling volcanic ash resuspension 6: Results for total suspended particle concentration (in mg/m 3 ) at the M13 meteorological station after applying the ction factor φ to the different emission schemes. The resulting values for φ are given in Table 2. 1: Location and altitude information of the stations used for the ground observation comparisons. Meteorological stations M1 to M17) provide visibility and wind velocity measurements. Air quality stations (from C1 to C3) provide direct urements of P M 10 concentration at CABA. All stations belong to Argentina, except M16 that is located in Uruguay.  Fig. 6. Results for total suspended particle concentration (in mg m −3 ) at the M13 meteorological station after applying the correction factor φ to the different emission schemes. The resulting values for φ are given in Table 2.
(spatially averaged u * ) above 0.4 m s −1 approximately. In the MB and SH schemes, u * t results from emission model parameterizations, whereas in the WE scheme this quantity is an input to be defined (see Eq. 7). Based on this previous analysis we set u * ts = 0.25 for the WE scheme that, once corrected by moisture, results in values of u * t ≈ 0.4 m s −1 . This agrees well with the value of 0.4 m s −1 used by Leadbetter et al. (2012) for Eyjafjallajökull in Iceland.

Ash transport modeling
For the resuspension runs, we set a FALL3D-7.0 domain spanning between 30 • -46 • S and 49 • -73 • W with a horizontal model resolution of 0.1 • and 14 vertical levels ranging from 10 to 4000 m above terrain (model layer thickness increases gradually in order to have finer resolution within the atmospheric boundary layer, where the ash cloud was mostly confined). The ash emission schemes in FALL3D-7.0 assume that the resuspended ash is distributed uniformly along the vertical up to a maximum height fixed by the user (250 m a.g.l. (above ground level) in our case). Particle granulometry changes at each deposit grid point, as determined by interpolating results of the preliminary eruption simulation. For resuspension, we consider the same 10 particle bins but fix the maximum size of particles that can be resuspended to 250 µm. This allows 7 particle classes to be resuspended and transported. However, larger particles have to be considered because they influence the emission rate in the SH scheme, although we verified that these cannot travel significant distances. Finally, in order to perform the evaluation of the kinematic turbulent fluxes, the diagonal components of the eddy diffusivity tensor in FALL3D have been parameterized using the similarity theory option for the ver-tical component and the CMAQ (Community Multiscale Air Quality) modeling system for the horizontal diffusion (Byun and Schere, 2006).

Model calibration
Our first quantitative comparisons indicated that, in most cases, the simulations overestimated observed surface particle concentrations. There might be multiple reasons for this, including (1) a poor characterization of the deposit, (2) the non-consideration of the surface roughness in the computation of the threshold friction velocity, (3) overestimation of the friction velocity by WRF-ARW or, (4) the inaccuracy of dust emission schemes when applied to volcanic ash. In order to calibrate the emission rates we considered observations made at the M13 meteorological station (see Table 1 for location information). We determined the correction factor φ for each emission scheme fitting simulation results to observations at this particular station. Results are shown in Fig. 6 and Table 2. Note that φ > 1 for the SH scheme and φ < 1 for WE and MB, meaning that the original SH formulation underestimated source strength whereas WE and MB overestimated. Although the factor φ was specifically determined from measurements at CABA, the value found was consistent with observations in the other meteorological stations, except for those in the deposit region (see Sect. 6.2).

Comparison with space-based measurements
In order to have a first qualitative verification of model results we used images from the Moderate Resolution Imaging Spectroradiometer (MODIS) sensor on-board the AQUA/TERRA satellites using the brightness temperature difference (BTD). The BTD algorithm, also known as the split window technique (Prata, 1989), computes the difference of brightness temperatures (derived from the inverse Plank's function) between the 10.7 and 12 µm frequency bands. Given the reverse absorption of volcanic ash particles (negative BTD), the technique allows for discrimination between meteorological and volcanic ash clouds. The efficiency of detection increases when the ash cloud is not opaque, has low water/ice contents, and contains small particles (Prata et al., 2001). Additionally, we considered volcanic ash advisories (VAAs) issued by the Buenos Aires VAAC. These are text messages that identify the observed and forecasted positions of ash clouds (ICAO, 2011). In particular, we used volcanic ash graphics (VAGs), a graphical depiction of the VAAs showing the edges of the polygons encompassing the ash cloud location. VAAs and VAGs from the Buenos Aires VAAC were based on observations from GOES-12 (Geostationary Operational Environmental Administration) and NOAA-19 satellites (METAR and airline pilot reports were considered in some cases).  Table 2. Characteristic quantities for each emission scheme. The correction factor φ is adjusted to match the maximum measured concentration of C TSP = 3.91 mg m −3 at the M13 station. The total emitted mass, mean ( e i ) and root mean square ( e 2 i ) of residuals at all stations (see Sect. 6.2) are also reported. These can be considered a measure of model to data bias and error respectively.  Figure 7 compares remote sensing observations and simulation results. The VAG polygons have been superimposed to the FALL3D predicted column mass at the appropriate times for both MB and SH schemes. The model predicts the formation and evolution of the cloud matching the satellite data, and shows how the ash located in northern Patagonia is transported towards the east and northeast reaching the coastal areas few hours after the passage of the frontal system. A large portion of central Argentina, the Río de La Plata region and southern Uruguay are affected. Note in Fig. 7a how the emission of ash starts shortly after the simulation begins, as expected from our previous analysis (see Fig. 5). The larger area of the VAG polygon in Fig. 7a is explained by the occurrence of a previous minor resuspension event, which was detected and reported by the VAAC. This event was partially modeled because we considered a model spin-up of 6 h.
In general, we found a good qualitative agreement between model results and BTD images. The results of the different schemes are similar, although the SH scheme gives higher concentrations close to the source region. It is interesting to note that, when the simulated cloud passes over the northeast of the Buenos Aires province (between 16 October at 19:00 UTC and 17 October at 01:00 UTC), the cloud is shifted ahead of the observed polygon (see Fig. 7d and e). These differences in cloud arrival times may be explained if one considers that WRF-ARW overpredicts wind velocities during this period (see Fig. 4c). Finally, it is worth noting that the little VAG triangle near the CCVC vent in Fig. 7d and e does not correspond to any resuspension phenomenon but reflects minor eruptive events from the CCVC remaining at that time.

Comparison with ground-based measurements
In this section we compare ground measurements of total suspended particulate matter concentration and PM 10 (denoted by C TSP and C PM 10 , respectively) with the FALL3D  values at the first model layer (10 m above ground level). To this purpose, we used information of visibility, wind speed and direction, and present and past weather phenomena reported in the meteorological message SYNOP (surface synoptic observations) and METAR (Meteorological Aerodrome Report) issued by the stations of the Argentinean National Meteorological Service network and the National Direction of Meteorology of Uruguay (see Fig. 3 and Table 1). The SYNOP is a numerical code (called FM-12 by WMO) used for reporting weather observations made by manned and automated weather stations. We used the visibility observations to yield an estimation of the TSP (total suspended particle) concentration using the following empirical relationship (Shao et al., 2003): where C TSP is the TSP concentration (in µg/m 3 ) and D is the visibility in km. Besides this data, we also considered measurements of respirable suspended particulate matter (PM 10 ) from the EPA Air Quality Monitoring Stations Network (GCBA, 2013). This network consist of three stations located at the northeast (C1), downtown (C2), and southwest (C3) of CABA. The instruments installed at these stations are the Thermo Model FH62 C14, which continuously measure the mass concentration of particulate utilizing a beta rays attenuation technique. The instruments meet US and International Particulate Monitoring Regulations and are US EPA certified to agree with the international air quality regulations. In the measurements, we subtracted a background value of 32 µg m −3 , corresponding to the averaged anthro-pogenic contribution recorded in October 2010, in order to estimate the contribution owing to the presence of ash. Additionally, we tried also other sources of data like LIDAR from CEILAP at Villa Martelli, province of Buenos Aires, and raw and processed data from the AERONET network but, because of the presence of clouds, we did not have goodquality data during the considered period. Figure 8 shows a comparison between simulated C TSP using the three different schemes and observations at 15 different stations. As observed from the figure, best quantitative agreement is obtained for stations registering the cloud arrival on 16 October (e.g., M4, M6 or M7). In contrast, for far-field stations (e.g., M10 or M17) where the plume arrived on 17 October, the model correctly predicts the arrival time but underestimates the observed values by a factor of 3-5. Conversely, for the stations close to the emission points (e.g., M1 at Bariloche or M2 at Neuquén, not shown) all the emission schemes overestimate substantially.
In order to evaluate if the cloud arrival times were correctly predicted by the model, we defined a characteristic event time as the time when C TSP reaches half of its maximum value. The comparison between simulated and observed characteristic event times shows good agreement, as observed in Fig. 9 for the WE scheme. Because the arrival time is mainly controlled by the meteorological model, characteristic event times given by the other two schemes were similar. Figure 10 compares modeled and observed maxima of C TSP (associated with the plume arrival at each station). Although some scattering exists, differences are in general within a factor of 2 (dashed line in Fig. 10). In order to quantify the performance of the different schemes, we define the residual of the ith observation as e i = (C mi − C oi )/C oi , where C mi and C oi are the modeled and observed values of C TSP , respectively. The mean ( e i ) and root mean square ( e 2 i ) of residuals are reported in Table 2 (note that e i = 0 is an indicator of model bias and e 2 i quantifies the error). As shown in Fig. 10, C TSP at points where the plume arrived on 17 October (triangles) are clearly underestimated by the model, contributing with negative residuals. From values in Table 2 we conclude that, for this particular episode, the WE and MB schemes give similar results, with WE showing the lower bias and error. In contrast, the SH scheme tends to overestimate C TSP (higher mean residuals) and shows larger differences with observations. Additionally, and regarding PM 10 concentration, we observed that MB and SH widely underestimate measurements in CABA (stations C1, C2, and C3) whereas WE gives the right order of magnitude. These low values can be attributed to several factors but, in our opinion, are explained by an incorrect prediction of the threshold friction velocity, which in the MB and SH schemes depends on soil moisture provided by WRF-ARW. In fact we also analyzed results using the MB scheme without soil  moisture correction, and found this gives much better results for PM 10 concentration in CABA (see Fig. 11).

Conclusions
We have implemented and tested three different dust emission schemes in the FALL3D-7.0 tephra dispersal model in order to simulate the 14-18 October 2011 ash resuspension events in central Patagonia. The modeling strategy combines a preliminary simulation to characterize the fallout deposit followed by the simulations of resuspension using FALL3D driven offline by WRF-ARW. Even if not specifically developed for volcanic ash, all the emission schemes give very promising results when comparing the simulated clouds with satellite images and the predicted concentrations of PM 10 and TSP at ground with measurements at different stations.
The three schemes present a similar qualitative behavior, but differ in the amount of emitted ash. Total resuspended mass is 8.4 × 10 10 kg for SH; 7.7 × 10 10 kg for MB and 3.3 × 10 10 kg for WE. For comparison, Collini et al. (2013) estimated a total amount of erupted ash mass of 1 − 5 × 10 12 kg during the period 4-19 June 2011. This amount of mass can give an idea of the hazards and disturbances caused by resuspended ash transported during the period considered in this study.
Remarkably, we find better agreement with observations using the simplest emission scheme (WE), although MB and SH showed also a good agreement. This result highlights the role of the sensitivity of more complex emission schemes to input parameters. In fact, relevant magnitudes for the source strength in the MB and SH schemes, such as granulometry or soil moisture, were obtained from modeling and are subject to large uncertainties. On the other hand, the simplest WE scheme seems more attractive from an operational point of view given its versatility. For example, an operator could easily set u * ts according to a specific situation in order to match partial observations with model results. Our next research steps are to explore how to better constrain inputs for more sophisticated schemes and to investigate how to improve model accuracy in the near-field areas.