Volcanic ash modeling with the online NMMB-MONARCH-ASH v 1 . 0 model : model description , case simulation , and evaluation

Traditionally, tephra transport and dispersal models have evolved decoupled (offline) from numerical weather prediction models. There is a concern that inconsistencies and shortcomings associated with this coupling strategy might lead to errors in the ash cloud forecast. Despite this concern and the significant progress in improving the accuracy of tephra dispersal models in the aftermath of the 2010 Eyjafjallajökull and 2011 Cordón Caulle eruptions, to date, no operational online dispersal model is available to forecast volcanic ash. Here, we describe and evaluate NMMB-MONARCH-ASH, a new online multi-scale meteorological and transport model that attempts to pioneer the forecast of volcanic aerosols at operational level. The model forecasts volcanic ash cloud trajectories, concentration of ash at relevant flight levels, and the expected deposit thickness for both regional and global configurations. Its online coupling approach improves the current state-of-the-art tephra dispersal models, especially in situations where meteorological conditions are changing rapidly in time, two-way feedbacks are significant, or distal ash cloud dispersal simulations are required. This work presents the model application for the first phases of the 2011 Cordón Caulle and 2001 Mount Etna eruptions. The computational efficiency of NMMB-MONARCH-ASH and its application results compare favorably with other long-range tephra dispersal models, supporting its operational implementation.


Introduction
Explosive volcanic eruptions can eject large quantities of particulate matter (tephra) that, along with other aerosol droplets and trace gases, are carried upwards into the atmosphere by the buoyant eruption column and are then dispersed by winds aloft (e.g., Sparks et al., 1997).Tephra particles smaller than 2 mm in diameter, technically defined as volcanic ash (Schmid, 1981), can spread over large distances away from the source, forming ash clouds that jeopardize air traffic (Casadevall, 1993), airports (Guffanti et al., 2009), and, for very large eruptions, alter both atmospheric composition and chemistry (Myhre et al., 2013;Self, 2006).Tephra transport and dispersal models (TTDMs; e.g., Folch, 2012) are used to simulate the atmospheric transport, dispersion, and ground deposition of tephra and to generate operational short-term forecasts to support civil aviation and emergency management.The recent eruptions of Eyjafjallajökull (Iceland) in 2010 and Cordón Caulle (Chile) in 2011 have reinforced the importance of tephra dispersal in the context of global aviation safety.In addition to short-term forecast, other model applications include the reconstruction of past events, studying the impact of volcano eruptions on climate, probabilistic tephra hazard assessments, or simulation of recent eruptions for model evaluation purposes.For any of those cases, TTDMs require a driving numerical weather prediction model (NWPM) or a meteorological reanalysis dataset for the description of the atmospheric conditions and an emission or source model for the characterization of the eruption column (Fig. 1).
Traditionally, TTDMs have evolved decoupled (offline) from NWPMs.In the offline strategy, the meteorological driver runs a priori and independently of the TTDM to produce the required meteorological fields at regular time intervals.Meteorological data are then furnished to the TTDM, which commonly assumes constant values for the meteoro- logical fields during each time interval or, at most, performs a linear interpolation in time.Although the offline approach is operationally advantageous, there is concern that it can lead to a number of accuracy issues (e.g., inaccurate handling of atmospheric processes) and limitations (e.g., neglect of feedback effects) that can be corrected by online approaches (Grell et al., 2004).These inconsistencies are especially important when meteorological conditions change rapidly in time or for long-range transport.However, uncertainties arising from offline systems have received little attention, even if the experience from other communities (e.g., air quality) highlights the importance of coupling online dispersal and meteorological models (e.g., Baklanov et al., 2014;Grell and Baklanov, 2011).To date, only the Weather Research and Forecasting model coupled to Chemistry (WRF-Chem; Grell et al., 2005) includes a coupled functionality that allows simulation of emission, transport, dispersion, transformation, and sedimentation of pollutants released during volcanic activities (Stuefer et al., 2013).
In this paper we describe and evaluate NMMB-MONARCH-ASH, a new online meteorological and atmospheric transport model for simulating the emission, transport, and deposition of ash (tephra) particles released from volcanic eruptions.The model predicts ash cloud trajectories, concentration of ash at relevant flight levels, and the expected deposit thickness for both regional and global domains.The novel online coupling in NMMB-MONARCH-ASH allows for solving of both the meteorological and aerosol transport concurrently and interactively at every time step.This coupling strategy aims at improving the current state-of-the-art tephra dispersal models, especially in situations where meteorological conditions are changing rapidly in time, twoway feedbacks are significant, or distal ash cloud dispersal simulations are required.The model builds on the Multiscale Online Nonhydrostatic AtmospheRe CHemistry model (NMMB-MONARCH; formerly known as NMMB/BSC-CTM) (Badia et al., 2017;Jorba et al., 2012;Pérez et al., 2011) to represent the transport of volcanic particles.Its meteorological core, the Non-hydrostatic Multiscale Model on a B grid (NMMB; Janjic and Black, 2007;Janjic and Gall, 2012;Janjic, 2005;Janjic et al., 2011), allows for nested global-regional atmospheric simulations by using consistent physics and dynamics formulations.The final objective in developing NMMB-MONARCH-ASH is twofold.On the one hand, at a research level, we aim at studying the differences between the online-offline modeling strategies.Moreover, a second version of the model is projected to quantify the feedback effects of dense volcanic ash clouds from large explosive eruptions on the radiative budget and local meteorology.On the other hand, at an operational level, the low computational cost of the NMMB dynamic core presented in this work suggests that NMMB-MONARCH-ASH could be applied for more accurate online operational forecasting of volcanic ash clouds.Consequently, the focus on developing an online volcanic ash model is timely.
The remainder of the paper is arranged as follows: Sect. 2 summarizes the modeling background and the standard physical schemes employed in NMMB-MONARCH-ASH; Sect. 3 provides a comprehensive description of the ashrelated modules, including details about the emission, transport, and deposition of volcanic particles; Sect. 4 validates the regional and global configurations of the model for the 2001 Mount Etna and 2011 Cordón Caulle long-lasting eruptions; Sect. 5 discusses the implementation and performance of the model for its operational use; and finally, Sect.6 provides a summary and conclusion of this work.
2 Modeling background NMMB-MONARCH-ASH is a novel online multi-scale meteorological and atmospheric transport model developed at the Barcelona Supercomputing Center (BSC).The model attempts to pioneer the forecast of volcanic aerosols by embedding a series of new modules in the BSC's operational system for short to midterm chemical weather forecasts (NMMB-MONARCH) developed at the BSC in collaboration with the US National Centers for Environmental Prediction (NCEP) and the NASA Goddard Institute for Space Studies.The development of the volcanic ash module follows the implementation of the mineral dust (Pérez et al., 2011) and sea-salt (Spada et al., 2013) modules in NMMB-MONARCH and allows for a range of different physical parameterizations for research and operational use.The system allows for feedback processes among gases, aerosol particles, and radiation and includes a gas-phase module to simulate tropospheric gasphase chemistry (Badia et al., 2017;Jorba et al., 2012).
Its meteorological core, the NMMB, is a fully compressible meteorological model with a non-hydrostatic option that allows for nested global-regional atmospheric simulations by using consistent physics and dynamics formulations.The standard physical and numerical schemes employed in NMMB are summarized in Table 1.The non-hydrostatic dynamics were designed to avoid overspecification.The cost of the extra non-hydrostatic dynamics is about 20 % of the cost of the hydrostatic part, both in terms of computer time and memory (Janjic, 2001(Janjic, , 2003)).The numerical schemes for the hydrostatic and nonhydrostatic options available in the NMMB dynamic solver were designed following the principles found in Janjic (1977), were developed and modified thereafter (Janjic, 1979(Janjic, , 1984(Janjic, , 2003)), and are summarized in Janjic and Gall (2012).The Arakawa B-grid horizontal staggering is applied in the horizontal coordinate, employing a rotated latitude-longitude coordinate for regional domains and latitude-longitude coordinate (Janjic, 2003) with polar filtering for global domains.Rotated latitude-longitude grids are employed for regional simulations in order to obtain more uniform grid distances.In this particular case, the Equator of the rotated system runs through the middle of the integration domain, reducing the longitudinal grid size as the southern and the northern boundaries of the integration domain are approached (Janjic and Gall, 2012).In the vertical, the Lorenz-staggering vertical grid is used with a hybrid sigma-pressure coordinate.The general time integration philosophy in NMMB uses explicit schemes when possible for accuracy, computational efficiency, and coding transparency (e.g., horizontal advection), and it used implicit schemes for very fast processes that would otherwise require a restrictively short time step for numerical stability with explicit differencing (e.g., vertical advection and diffusion, vertically propagating sound waves).The NMMB model became the North American Mesoscale (NAM) operational meteorological model in October of 2011, and it has been computa-tionally robust, efficient, and reliable in operational applications and pre-operational tests since then.In high-resolution numerical weather prediction (NWP) applications, the efficiency of the model significantly exceeds those of several established state-of-the-art non-hydrostatic models (e.g., Janjic and Gall, 2012).

The volcanic ash module
The ash module is embedded within the NMMB meteorological model and solves the mass balance equation for volcanic ash, taking into account (i) the characterization of the source term (emissions), (ii) the transport of volcanic particles (advection-diffusion), and (iii) the particle removal mechanisms (sedimentation-deposition).The coupling strategy of this module can be turned on or off, depending on the solution required (online vs. offline).The online version of the model solves both the meteorological and aerosol transport concurrently and consistently (online coupling).This strategy allows the particle transport to be automatically tied to the model resolution time and space scales, resulting in a more realistic representation of the meteorological conditions.In contrast, the offline approach uses an "effective wind field" in which meteorological conditions (e.g., wind velocity, mid-layer pressure) are set to constant and are only updated at specific coupling intervals (i.e., time for which meteorological fluctuations are not explicitly resolved).This strategy replicates the offline coupling effect of traditional dispersal models used at operational levels (e.g., coupling intervals of 1 or 6 h).The conservativeness of the model is evaluated to ensure that the ash transport scheme is consistent with the mass conservation equation.

Source term
Explosive volcanic eruptions release large amounts of particles into the atmosphere.These particles, commonly known as tephra, mix with ambient air to form an eruption column or volcanic plume.To forecast the ash cloud movement and provide actual ashfall concentrations, tephra dispersal models require a complete characterization of the parameters describing the source term.These parameters are generally referred to as eruption source parameters (ESPs) and include the eruption start and duration, column height, mass eruption rate (MER), vertical distribution of mass, and the particle grain size distribution (GSD).ESPs vary not only from one eruption to another but also during the different eruptive phases of a single event.
Typically, the eruption starting time, duration, and column height are inferred or constrained from visual or satellite observations.However, other parameters such as GSD, MER, or the vertical distribution of mass in the column are not available in real time and must be inferred from previous events of similar characteristics (e.g., Mastin et al., 2009 Monin and Obukhov (1954); Obukhov similarity theory with Zilitinkevich thermal Zilitinkevich (1965); Janjic (1994Janjic ( , 1996) ) roughness length Land surface, heat, and moisture LISS model Vukovic et al. (2010) surface flux Planetary boundary layer-free Mellor-Yamada-Janjic scheme Mellor and Yamada (1982); Janjic atmosphere (1996,2001)

Convective adjustments
Betts-Miller-Janjic scheme Betts and Miller (1986); Janjic (1994Janjic ( , 2000) ) Uncertainties in source parameter values are a key factor limiting the accuracy of ash-cloud model forecasts (Bonadonna et al., 2015a).The characterization of each ESP in NMMB-MONARCH-ASH is described in the following subsections.

Mass eruption rate
The MER gives the mass released by unit of time and defines the eruption intensity.Its characterization in NMMB-MONARCH-ASH is achieved by employing a series of empirical correlations between (observed) column height and eruption rate, which, according to plume similarity theory, scale roughly as the fourth power of the column height.Be-cause of this strong dependence, uncertainties within 20 % in the determination of column height can translate into uncertainties up to 70 % for the MER (e.g., Biass and Bonadonna, 2011).Averaged column heights of eruptions that have not been directly observed are typically derived from characteristics of tephra deposits (e.g., Bonadonna and Costa, 2013;Carey and Sparks, 1986;Pyle, 1989) or are derived from model inversion (e.g., Connor and Connor, 2006;Marti et al., 2016;Pfeiffer et al., 2005).

Vertical distribution of mass
The vertical distribution of mass in the column at the vent location is key when representing the plume, especially if wind shear exists with elevation at the volcano (Lin, 2012).To determine this distribution of mass, NMMB-MONARCH-ASH allows for the following geometrical distributions: (i) point source, where mass is released as a single source point at a certain height above the vent, H plume ; (ii) top hat, where mass is released along an umbrella-type slab of user-defined thickness; and (iii) the so-called Suzuki distribution (Suzuki, 1983;Pfeiffer et al., 2005), which assumes a more complex vertical distribution of mass release along the eruption column; where S is the mass per unit of time (kg s −1 ) released at a given height z above the vent, MER is the total mass eruption rate, H plume is the column height above the vent, and A and λ are the so-called Suzuki parameters.The parameter A dictates the height of the maximum particle release (concentra-tion), whereas λ controls how closely mass distributes around this maximum.Any of the three options above can be combined independently with the different options for MER estimation.In NMMB-MONARCH-ASH, the terrain following hybrid sigma-pressure vertical levels of the model must be converted to elevations for each model integration time step in order to interpolate MER from the discrete source points into the nodes of the model grid.

Grain size distribution
The impact of explosive volcanic eruptions on climate and air traffic strongly depends on the concentration and GSD of pyroclastic fragments injected into the atmosphere (e.g., Girault et al., 2014).Grain size distribution is normally reconstructed by volcanologists from grain size data at individual outcrops, ranging from basic unweighted average of the GSD at individual sparse outcrops to various integration methods of grain size data (e.g., Rose and Durant, 2009).The particle grain size distribution in NMMB-MONARCH-ASH is specified through an input file, which defines the particle bin properties (bin mass fraction, diameter, density, and shape factor).In volcanology, grain size distributions are given in terms of the , defined as d = 2 − , where d is the particle diameter in mm.The granulometry file in the model can be furnished by the user (typically derived from field data) or generated by an external utility program that produces Gaussian and bi-Gaussian distributions in (log-normal in diameter d) (Costa et al., 2016;Folch et al., 2009).

Particle aggregation
The total grain size distribution (TGSD) erupted at the vent can be altered in case of particle aggregation, which dramatically impacts particle transport dynamics, thereby reducing the atmospheric residence time of aggregating particles and promoting the premature fallout of fine ash.For computational purposes, particle aggregation in NMMB-MONARCH-ASH is assumed to take place mainly in the eruption column, where particle concentration and water contents are higher (the subsequent formation of aggregates downstream in the ash cloud under the appropriate atmospheric conditions is not contemplated by the model).The model considers aggregates as another particle class (bin), introduced as a standard source term by solving either (i) a series of simple analytical expressions based on field observations or (ii) a more sophisticated wet aggregation model originally proposed by Costa et al. (2010).
The analytical expressions available in the model modify the user-given particle grain size distribution by assuming that a certain mass fraction of each granulometric class forms a new aggregate class added to the TGSD.Despite the obvious limitations (obviates the physics of aggregation processes), these field-based simplistic approaches are advantageous in that only the source term has to be modified in order to account for aggregation.Table 3 provides an overview of these options.In addition to these empirical aggregation schemes, NMMB-MONARCH-ASH also includes a wet aggregation model originally proposed by Costa et al. (2010).This option allows for wet aggregation in the column, providing an intermediate solution between the unaffordable allsize class approach and the empirical solutions presented before.The model is based on a solution of the classical Smoluchowski equation, obtained by introducing a similarity variable and a fractal relationship for the number of primary particles in an aggregate.It also considers three different mechanisms for particle collision: Brownian motion, ambient fluid shear, and differential sedimentation.Table 4 provides an overview of the governing equations of this wet aggregation model.

FPlume model
A more sophisticated approach to obtain MER and the mass distribution in the column from the conditions at the vent consists of solving a 1-D radially averaged buoyant plume theory (BPT) model for mass, momentum, and energy.These 1-D plume models are more useful in operational roles and broad exploratory investigations (Costa et al., 2015;Devenish et al., 2012).For that reason, NMMB-MONARCH-ASH is coupled with the 1-D FPlume model (Folch et al., 2015), a 1-D cross-section-averaged plume model that accounts for plume bent over, entrainment of ambient moisture, effects of water phase changes on the energy budget, particle fallout and re-entrainment by turbulent eddies, and variable entrainment coefficients fitted from experiments.The model also accounts for particle aggregation in the presence of liquid water or ice that depends on column dynamics, particle properties, and amount of liquid water and ice existing in the column (Folch et al., 2010).This allows the plume model to predict an effective grain size distribution depleted in fines with respect to that erupted at the vent.For a complete definition of the governing equations of FPlume, refer to Folch et al. (2015).FPlume has two solving strategies where the model (i) solves directly for column height for a given MER or (ii) solves iteratively for MER for a given height.For any case, the following inputs need to be provided in the ash input file in NMMB-MONARCH-ASH: eruption start and duration, vent coordinates and elevation, conditions at the vent (exit velocity, temperature, magmatic water mass fraction, and total grain size distribution), and total column height or mass eruption rate.

Particle advection-diffusion
Transport of volcanic ash by advection and turbulent diffusion is analogous to that of atmospheric tracer (e.g., moisture) transport (Janjic et al., 2009) in NMMB.Tracer advection is Eulerian, positive-definite, and monotonic.The Adams-Bashforth scheme is used for horizontal advection and the Crank-Nicolson scheme is used for vertical advection.For the horizontal diffusion, the model uses a secondorder scheme with two types of parameterized dissipative processes: explicit lateral diffusion (often called horizontal diffusion, a second-order nonlinear Smagorinsky-type approach; Janjic, 1990) and horizontal divergence damping (Janjic and Gall, 2012).Plumes from high-intensity eruptions can be injected high into the stratosphere, reaching a maximum column height and intruding laterally at the neutral buoyancy level (NBL) as a gravity current (Sparks et al., 1997).This current can spread at velocities exceeding those of the surrounding winds, affecting tephra transport and deposition near the source.As larger particles are removed by deposition and air is entrained, the plume density decreases and momentum reduces such that, at a certain distance, atmospheric turbulence and wind advection become the dominant atmospheric transport mechanisms (Baines and Sparks, 2005).Neglecting the gravitational spreading of the umbrella cloud in tephra dispersal simulations could misrepresent the interaction of the volcanic ash cloud and the atmospheric wind field for highintensity eruptions and for proximal deposition of tephra (Mastin et al., 2014).To account for the gravity-driven transport, NMMB-MONARCH-ASH is coupled with the model of Costa et al. (2013), describing cloud spreading as a gravity current.This parameterization calculates an effective ra-dial velocity of the umbrella spreading as a function of time or cloud radius.The effective radial velocity of the umbrella spreading is then combined with the wind field velocity centered above the vent in the umbrella region to calculate the contribution of the gravitational spreading to the total cloud spreading.To estimate the radial distance at which the critical transition between gravity-driven and passive transport occurs, the umbrella front velocity is compared with the mean wind velocity at the NBL, estimating the Richardson number.Table 5 provides an overview of the governing equations of the gravity current model embedded in NMMB-MONARCH-ASH.

Particle sedimentation and dry deposition
Particle sedimentation in NMMB-MONARCH-ASH is governed by the terminal velocity of settling particles.This fall velocity is sensitive to particle size and atmospheric conditions, determining the residence time of ash particles in the atmosphere.NMMB-MONARCH assumes that the settling velocities of aerosols (mineral dust, sea salt, etc.) follow the Stokes law for spherical particles corrected by the Cunningham slip factor.The Stokes law applies to the creeping or Stokes flow regime, in which the drag force is proportional to particle velocity and holds only for Reynolds numbers Re ≤ 0.1.This regime is justified for small particles  (Suzuki and Koyaguchi, 2009) spreading N = Brunt-Väisälä frequency (atm.ambient stratification) q = volumetric flow rate into the umbrella region (10) M = efficiency of air entrainment rate into the 32 Table 5. Governing equations for NMMB-MONARCH-ASH gravity current model. !"!!

Gravity current scheme Eq. Parameters
Effective radial velocity of the umbrella spreading  Ri > 1: gravity-driven regime passive transport) Ri < 0.25: passive transport regime and aerosols (< 20 µm).However, calculating fallout times based on settling according to Stokes law is less adequate for coarse ash (> 64 µm), which sediments much faster.In addition, ash particles are not spherical, which complicates and further slows fallout.In order to properly simulate a wider spectrum of particle sizes, NMMB-MONARCH-ASH adds a new sedimentation module that covers the turbulent regime (Re ≥ 1000) in which the drag force is proportional to the square of the particle velocity.In this case, the gravitational particle settling velocity, v s (in m s −1 ), can be expressed as where ρ a and ρ p denote air and particle density, respectively, d is the particle equivalent diameter, and C d is the drag coefficient (depending on the Reynolds number).Strictly, the expression above is valid for spherical particles in the turbulent regime, but it is often generalized to the whole range of Re numbers and particle shapes by defining the drag coefficient properly.Table 6 provides an overview of the different settling velocity models available in NMMB-MONARCH-ASH, each relying on different empirical evaluations of drag coefficient.Dry deposition, acting at the bottom layer of the model, is a complex process depending on physical and chemical properties of the particle, the underlying surface characteristics, and micro-meteorological conditions.Dry deposition in the model is based on that originally proposed by Zhang et al. (2001).This parameterization has been updated to account for the different settling velocities available for volcanic particles -Eq.( 13).The dry deposition velocity in the model, v d (in m s −1 ), is given by where R a is the aerodynamic resistance of the particle, and R s is the surface resistance (both in s m −1 ).These terms take into account all the effects of the lowermost layer of the atmosphere, such as turbulence (R a ), Brownian diffusion, impaction, and interception (R s ).It is worth mentioning that for most of its residence time, airborne volcanic ash lies above the near-surface atmospheric layers where gravitation dominates, implying that, in most cases, dry deposition has little influence on model results.

Mass conservation
Mass conservation is a critical requirement for any atmospheric transport algorithm.Nonconservative schemes can significantly underestimate or overestimate concentrations, especially for long-time integrations in which it is critical that the tracer advection scheme is consistent with the mass continuity equation (Jöckel et al., 2001).Most mesoscale meteorological models use observation and/or analyzed fields or global model results as initial conditions, and therefore they are not very sensitive to slowly accumulated mass inconsistencies as reinitializations remove accumulations.However, For spherical particles only GANSER (Ganser, 1993) For spherical and non-spherical particles = shape factors = average between the min and max.axis = sphere volume = particle sphericity (= 1 for spheres) WILSON (Pfeiffer et al., 2005;Wilson and Huang, 1979) , ,  = particle aspect ratio  =  +  2 !! = particle semi-axes DELLINO (Dellino et al., 2005) For larger particles only = Archimedes number = gravity acceleration = particle shape factor = dynamic viscosity = particle equivalent diameter, = particle density = air density (14) For spherical particles only (Arastoopour et al., 1982) GANSER (Ganser, ReK 1 K 2 (15) For spherical and nonspherical particles 1993) For spherical particles only GANSER (Ganser, 1993) For spherical and non-spherical particles = shape factors = average between the min and max.axis = sphere volume = particle sphericity (= 1 for spheres) WILSON (Pfeiffer et al., 2005;Wilson and Huang, 1979) , ,  = particle aspect ratio  =  +  2 !! = particle semi-axes DELLINO (Dellino et al., 2005) For larger particles only = Archimedes number = gravity acceleration = particle shape factor = dynamic viscosity = particle equivalent diameter, = particle density = air density dispersal models are usually very sensitive to mass inconsistencies set in previous simulations or spin-up fields as initial conditions, thereby accumulating mass inconsistencies.In addition to mass conservation, monotonicity and prevention of nonphysical under-and overshoots in the solution are also highly desirable characteristics in transport schemes (Rood, 1987).For these reasons, the model includes a conservative, positive-definite (i.e., tracer is a positive scalar), and monotone (i.e., entirely increasing) Eulerian scheme for advection.
The positive definiteness in the model is guaranteed by advecting the square root of the tracer using a modified Adams-Bashforth scheme for the horizontal direction and a Crank-Nicolson scheme for the vertical direction.The conservation of the tracer is achieved as a result of the conservation of quadratic quantities by the advection scheme.Monotonization is applied a posteriori to eliminate new extrema (Janjic et al., 2009).The conservative nature of NMMB-MONARCH-ASH is evaluated by calculating the mass flux at the boundaries (for regional domains) of the computational domain, the airborne mass, and the mass deposited on the ground to verify mass conservation at each time step (e.g., < 0.5 % mass creation for a 30-day simulation).

Numerical performance
The high computational efficiency of the NMMB meteorological driver allows for the application of nonhydrostatic dynamics at a global scale (Janjic et al., 2009) and supports the use of NMMB-MONARCH-ASH in an operational forecast of volcanic ash clouds.Model parallelization is based on the well-established Message Passing Interface (MPI) library.The computational domain is decomposed into subdomains of nearly equal size in order to balance the computational load, where each processor is in charge of solving the model equations in one subdomain.The Eulerian schemes in the model require relatively narrow and constant width halos (i.e., data points from the computational domain of neighboring subdomains that are replicated locally for computational convenience), which simplify and reduce communications.
To measure the time to solution required, we compute the parallel speed-up (computation speed) of the model, that is, the performance gains of parallel processing in comparison to serial processing: where S is the computed speed-up value and t is the simulation run time, employing P processors instead of running it serially (P = 1).
To evaluate the efficiency of the model while using the computational resources, the parallel efficiency of the model is computed by looking at the ratio between the parallel speed-up over P : Parallel efficiency is used as a metric to determine how far the model's speed-up is from the ideal.If the speed-up is ideal, the efficiency is 1, regardless of how many cores the program is running on.If the speed-up is less than ideal, the efficiency is less than 1.

Simulations and validation
The forecast skills of NMMB-MONARCH-ASH have been tested for several well-characterized eruptions, including the Pinatubo 1991 (Philippines), Etna 2001 (Italy), Chaitén 2008 (Chile), and Cordón Caulle 2011 (Chile) eruptions (e.g., Marti et al., 2013Marti et al., , 2014)).Here, we present two applications of the model for the ash dispersal forecast of weak and strong long-lasting eruptions.Section 4.1 summarizes the results of the regional and global simulations for the first days of the 2011 Cordón Caulle eruption.This event represents a suitable case study of strong long-lasting eruptions with changing winds, which is useful to evaluate the advantages of the online approach for operational forecast.In a parallel effort, Sect.4.2 summarizes the results from the regional configuration of the model for the 2001 Etna eruption.This eruption is a good example of a weak, long-lasting eruption, useful when evaluating the sedimentation mechanisms of the model against well-characterized tephra deposits.

The 2011 Cordón Caulle eruption
The 2011 Cordón Caulle eruption was a typical midlatitude Central and South Andean eruption, where dominating winds carried ash clouds over the Andes, causing abundant ash fallout across the Argentine Patagonia.In addition to the significant regional impacts on agriculture, livestock, and water distribution systems, this eruption stranded thousands of passengers due to air traffic disruptions in the Southern Hemisphere, thereby causing important economic losses to airlines and society (e.g., Raga et al., 2013;Wilson et al., 2013).This event is evidence of the global nature of the volcanic ash dispersion phenomena and highlights the need for accurate realtime forecasts of ash clouds.The Cordón Caulle volcanic complex (Chile, 40.5 • S, 72.2 • W, vent height 1420 m a.s.l.) reawakened on 4 June 2011 around 18:30 UTC after decades of quiescence.The initial explosive phase spanned more than 2 weeks, generating ash clouds that dispersed over the Andes.The climactic phase (∼ 27 h) (Jay et al., 2014) was associated with a ∼ 9 km (a.s.l.) high column (Osores et al., 2014).For the period between 4 -14 June, numerous flights were disrupted in Paraguay, Uruguay, Chile, southern Argentina, and Brazil.The two major airports serving Buenos Aires and the international airport in Montevideo, Uruguay, were closed for several days, along with airports in Patagonia (Wilson et al., 2013).A detailed chronology of the eruption can be found in Collini et al. (2013) and Elissondo et al. (2016), the stratigraphy and characteristics of the resulting fallout deposit are described in Pistolesi et al. (2015) and Bonadonna et al. (2015b), and a summary of the environmental impacts of the eruption is discussed in Raga et al. ( 2013) and Wilson et al. (2013).
Here, we describe the synoptic meteorological situation during the first 2 weeks of eruptive activity (Fig. 2), and we give a brief chronology of the events in order to compare them with the predictions of the model.The eruption developed as a long-lasting rhyolitic activity with plume heights above the vent between 9 and 10 km high a.s.l.(4-6 June), 4 and 9 km during the following week (7-14 June), and < 6 km after 14 June (Global Volcanism Program, GVP, http://www.volcano.si.edu;Siebert et al., 2010).The first major episode, on 4 June (18:45 UTC), resulted in an ash cloud (9-10 km) that reached the Chile-Argentina border within 1 h of the eruption.On 5 June, E-SE winds drove the plume to the Atlantic Ocean (1800 km away from the source), leaving a large area of Argentina territory affected by ashfall.On 6 June, the plume changed its direction abruptly toward the N-NE, reaching the northern regions of the Argentine Patagonia, and then shifted direction again towards the SE, threatening the Buenos Aires air space.On 7 June, a second episode resulted in a plume (4-9 km) dispersing ash further to the north of Argentina, leading to a more recognizable shift of winds to the E-SE.On 8 June, the volcanic cloud (9-10 km a.s.l.) dispersed towards the NE with a bend toward the SE 400 km from the source.On 9 June, the plume had a NE direction, reaching the city of Buenos Aires and the northern boundary of Paraguay following a frontal zone passage through Patagonia.This resulted in major air traffic disruption at the two international airports that service the city: Aeroparque (AEP) and Ezeiza (EZE), which remained intermittently closed during the following 15 days.Later during the day, the wind turned SE, dispersing ash over Uruguay, Brazil, and Paraguay.Ash cloud continued to change in direction over the next 6 days, with clouds following the ridge structure to the NE and SE.Plots show the direction (vector) and velocity (contours m s −1 ) of the wind at 9100 m a.g.l.(above ground level) (300 hPa circa).Meteorological data obtained from the NMMB meteorological forecast driven with ERA-Interim reanalysis at 0.75 • horizontal resolution.

Model setup
The model domain for the regional run is presented in Table 7 and consists of 268 × 268 grid points covering the northern regions of Chile and Argentina using a rotated latitudelongitude grid with a horizontal resolution of 0.15 • × 0.15 • and 60 vertical layers.The top pressure of the model was set to 21 hPa (∼ 34 km) with a mesh refinement near the top (to capture the dispersion of ash) and the ground (to capture the characteristics of the atmospheric boundary layer).The computational domain spans in longitude from 41 to 81 • W and in latitude from 18 to 58 • S. Runs were performed with the online version of NMMB-MONARCH-ASH from 3 June 2011 at 00:00 UTC to 21 June 2011 at 00:00 UTC.The integration time step for the meteorological core and aerosol transport was set to 30 s.The dynamic time steps for the long-and shortwave radiations were computed every 120 time steps.Feedback effects of ash particles on meteorology and radiation were not included in this run.The meteorological driver was initialized with wind fields from the ERA-Interim reanalysis at 0.75 • × 0.75 • resolution as initial and 6 h boundary conditions.In order to reduce the errors in meteorological conditions, they were reinitialized every 24 h with a spin-up of 12 h.Daily eruption source parameters (ESP) were obtained from Osores et al. (2014), who estimated column heights for each eruptive pulse using the imager sensor data from the GOES-13 satellite, applying the cloud-top infrared (IR) image technique (Kidder and Vonder-Haar, 1995).Mass flow rate released along the column was derived from column heights based on Mastin et al. (2009), assuming a Suzuki vertical distribution of mass typical of explosive Plinian eruptions (A = 4; λ = 5).Grain size distribu- tion was obtained from Collini et al. (2013) and discretized in 10 bins ranging from −1 (2 mm) to 8 (4 µm), with a linear dependency of particle density on diameter ranging from 1000 to 2200 kg m −3 .Particle sphericity was set to a constant standard value of 0.9 for all bins.The percentage aggregation model was used to update the TGSD with a new bin for aggregates, resulting in a total of 11 bins.

Validation of results against satellite imagery
Model results for the airborne mass concentration of ash were validated using qualitative and quantitative comparisons with data obtained using two different techniques.We performed a qualitative comparison between the simulated column mass (g m −2 ) from the model and the NOAA-AVHRR satellite imagery provided by the high-resolution picture transmission (HRPT) division of the Argentinian National Meteorological Service. Figure 3 shows how the model predictions for cloud trajectory and arrival times are in agreement with observations, capturing the three major dispersion episodes.It should be noted that these types of images are not directly comparable because the MODIS ash detection threshold and the reflectivity coefficients of volcanic ash are not well constrained.However, the figure illustrates the capability of the model to predict the variation of the cloud position with time.
Column mass simulations were also validated against ash mass loadings presented by Osores et al. (2015), who retrieved ash-contaminated pixels detected on the basis of the concept of reverse absorption (Prata, 1989a, b), i.e., those pixels with brightness temperature differences between 11 and 12 µm (BTD 11-12 µm) that are lower than 0.0 K. To minimize the presence of false positives, pixels with a BTD 11-12 µm > −0.6 K and clear sky pixels were removed.Mass loadings were mapped up to 15 g m −2 based on an approach that combines the satellite data with lookup tables of brightness temperatures obtained with a radiative transfer model and optical properties of andesite volcanic rocks (Prata, 2011).Figure 4 shows a good quantita-  tive agreement between the model results and the airborne ash mass loadings described above.

Validation of results against fallout deposit
Tephra was mostly deposited eastward from the source during the first 72 h of the event within an elongated area between 40-42 • S and 64-72 • W. Results from the NMMB-MONARCH-ASH ash deposition forecast were validated against (i) a detailed characterization of the proximal deposit for the first 72 h of the eruption and (ii) an isopach map derived from measurements taken for the period beginning on 4 until 30 June (Collini et al., 2013).
To evaluate the simulated computed thicknesses (centimeters) by the model near the vent during the first 72 h of the event, model results were compared against a comprehensive classification of the proximal deposit presented by Pistolesi et al. (2015), who constrained the stratigraphic sequence of the deposit in different units (phases).Here, we constrain the deposit to the first three units of their work, corresponding to the first 72 h of the eruptive even and including (i) Unit I, which contains coarser-grained layers A-B, representing the very first stage of the eruption within the first 50 km from the vent, and layers A-F associated with the first 24-30 h of the eruption (afternoon of 4 to morning of 5 June); (ii) Unit II, which contains layer H, a fine pumice lapilli layer that was emplaced starting on the night of 6 June; (iii) Unit III, which encloses layer K2, the easiest to identify from several coarser (fine-lapilli) grain size layers, and being associated with the morning of 7 June.Figure 5 shows that NMMB-MONARCH-ASH can reproduce the deposit presented by Pistolesi et al. (2015) both in time and space.Key sections located along the dispersal area (e.g., San Carlos de Bariloche -SCB, 90 km from the vent; Ingeniero Jacobacci -IJ, 240 km east of the vent) were used as geographic references.
To evaluate the model performance at the end of our simulation, model results were also validated against an isopach map derived from measurements taken from the 4 to 30 June presented by Collini et al. (2013).Deposit load variations produced by remobilization were not considered in this analysis.Figure 6 shows good agreement between the modeled deposit load (kg m −2 ) at the end of the simulation and the measured ground deposit isopachs (kg m −2 ) on 30 June from Collini et al. (2013).
The model resulted in a cumulative mass of ∼ 4.2 × 10 11 kg.This value is in agreement with previous works, where total mass was either modeled (Collini et al., 2013) or estimated by empirical fits (Bonadonna et al., 2015b).Ashfall forecast with the NMMB-MONARCH-ASH model represented the overall deposit load for the 2011 Caulle eruption well.

Global simulation
For this simulation, the global domain was configured using a regular latitude-longitude grid with a horizontal resolution of 0.75 • × 1 • and 60 vertical layers.The ash distribution is simulated between 3 and 21 June 2011 using the ERA-Interim reanalysis at 0.75 • × 0.75 • resolution as initial and 6 h boundary conditions.Meteorological conditions for the global runs were also reinitialized every 24 h.The atmospheric model's fundamental time step was set to 180 s, while the rest of the model variables and grain size distribution remained the same as in the regional simulation.Figure 7 shows the global dispersal of ash for the 2011 Cordón Caulle eruption at different times of the simulation.As it can be inferred from this figure, by 10 June, the plume entered the Australian and New Zealand airspace (Fig. 7b), covering more than half of the Southern Hemisphere.At that point, the Civil Aviation Authority of New Zealand warned pilots that the ash cloud was between 20 000 and 35 000 ft (6 to 11 km), the average cruising level for many aircraft (Sommer, 2011).Before the end of our simulation, on 13 June the ash cloud had completed its first circle around the globe.This is in agreement with satellite images reported by the Darwin Volcanic Ash Advisory Centre (Darwin VAAC, 2011).Finally, results from the global simulation are also in agreement with those from our regional run.

Forecasting impacts on civil aviation
NMMB-MONARCH-ASH can furnish values of airborne concentration at relevant flight levels (FL), defined as the vertical altitude (expressed in hundreds of feet) at standard pressure at which the ash concentration is measured.This information is particularly important for air traffic management and can be used to decide alternative routes to avoid an encounter with a volcanic cloud.Airborne concentration at FL050 (5000 ft on nominal pressure ∼ 850 hPa in pressure altitude) is relevant for the determination of flight cancelations and airport closures, while concentrations at FL300 (30 000 ft) are critical to assist flight dispatchers while planning flight paths and designing alternative routes in the presence of a volcanic eruption.The model runs as if responding to an eruptive event, i.e., we only used the semiquantitative data available at that time as volcanological inputs.
Figure 8 shows the airspace contamination forecasted by the model during 6-7 June at flight levels FL050 and FL300, within a latitude band between 20 and 55 • S. Model results show the volcanic cloud twisting in different directions during that period of time, achieving critical concentration values within a wide area east of the Andes range.On 6 June, simulation results show the volcanic cloud at high atmospheric pressure (∼ 30 000 ft or 300 hPa) moving northwards and the cloud at lower atmospheric pressure threatening the main international airports that service the region of Buenos Aires (Fig. 8a).On the morning of 7 June, the ash cloud present at lower atmospheric pressure changed its direction towards the SW, ultimately affecting part of Patagonia and Chile (Fig. 8b), while higher ash clouds started their course around the globe (Fig. 8c).These results suggest that the cancellation of multiple flights in several Argentinean airports during this time was justified.It is important to point out that, for this work, our objective is not to perform a detailed study of the Caulle eruption but to use it as a blind test to confront short-term model predictions and semiquantitative syn-eruptive observations.

The 2001 Mount Etna eruption
Mount Etna is the most active volcano in Europe and constitutes a continuous hazard for eastern Sicily.Since 1980, Mount Etna has injected large volumes of pyroclasts into the atmosphere (between 10 4 and 10 7 m 3 per event) over more than 160 eruptive episodes (Scollo et al., 2012).The explosive activity of Mount Etna reached its climax in 2001 and 2002-03 when two major flank eruptions occurred, both characterized by long-lasting explosive activity (Branca and Del Carlo, 2005).The 2001 event represents a good case to evaluate the deposition mechanisms of the model against the well-characterized tephra deposit reported in Scollo et al. (2007).The explosive activity at the 2570 m vent had three main phases characterized by phreatomagmatic, magmatic, and vulcanian explosions.The eruption started with a series of phreatomagmatic explosions during the first days of the eruption.These explosions were followed by a second eruptive phase characterized by strombolian and Hawaiian style explosions during 19-24 July.The explosive activity continued until 6 August with a series of vulcanian explosions.Tephra fallout associated with the explosive activity during 21-24 July represented a major source of hazard for eastern Sicily.Flight operations were canceled at the Catania and Reggio Calabria airports on 22 and 23 July.A detailed chronology of the eruption can be found in Scollo et al. (2007).Volcanic plumes were captured by the Multi-angle Imaging SpectroRadiometer (MISR) on board NASA's Terra spacecraft and were analyzed with stereo matching techniques to evaluate the height of the volcanic aerosol with a precision of a few hundred meters (Scollo et al., 2012).
Here, we validate NMMB-MONARCH-ASH against the tephra deposit produced from the 2570 m vent for that period of time, and we compare the model performance against simulation results from the FALL3D model (Costa et al., 2006;Folch et al., 2009) for the same event.FALL3D is a Eulerian model for transport and deposition of volcanic ash particles that solves a set of advection-diffusion-sedimentation equations (one equation for each particle class) on a structured terrain-following grid using a second-order finite differences explicit scheme.The FALL3D model is used at the Buenos Aires and Darwin Volcanic Ash Advisory Centers (VAAC) in operational forecasts.

Regional simulation
Model setup Two regional domains were used to simulate the first phase of the 2001 eruption of Mount Etna (Table 8).The first domain (Regional 1), used to reconstruct the tephra deposit, consists of 101 × 101 grid points covering the SE flank using a rotated latitude-longitude grid with a horizontal resolution of 0.05 • × 0.05 • and 60 vertical layers.Similar to the Cordón Caulle simulations, the top pressure of the model was set to 21 hPa (∼ 34 km) with a mesh refinement near the top and ground.The computational domain spans in longitude from 12.5 to 17.5 • E and in latitude from 35.25 to 40.25 • N. Simulation runs were performed with the online version of NMMB-MONARCH-ASH from 21 July 2001 at 00:00 UTC to 25 July 2001 at 00:00 UTC.The integration time step for the meteorological core was set to 10 s.The meteorological driver was initialized with ERA-Interim reanalysis meteorological data at 0.75 • × 0.75 • resolution as initial and 6 h boundary conditions.A spin-up of 12 h was used to prepare the meteorological conditions for the run.Each daily model run was reinitialized with the corresponding reanalysis of the model tracers' output from the previous day and the associated eruption source parameters.Meteorological conditions were reinitialized every 24 h.The grain size distribution and eruption source parameters were obtained from Scollo et al. (2007), who assumed a Suzuki vertical mass dis- tribution located at the middle of the eruption column (A = 2; λ = 1) and employed the Mastin et al. (2009) empirical relationship to characterize the MER and the Voronoi tessellation method to obtain the grain size distribution.Finally, sensitivity analyses were performed against the different aggregation schemes available in the model.In all cases, the TGSD was updated with a new bin for aggregates, resulting in a total of 8 bins.
A second regional domain (Regional 2) was used to evaluate tephra dispersal between 21 and 25 July.In this case, the domain consisted of 201 × 201 grid points covering a computational domain spanning in longitude from 41 to 81 • E and in latitude from 18 to 58 • S.This domain used a coarser horizontal resolution of 0.1 • × 0.1 • and 60 vertical layers.The integration time step for the meteorological core was set to 30 s.The rest of the model setup was kept the same as in the first regional domain (Regional 1).

Validation of results against fallout deposit
At the end of the second explosive phase, a continuous tephra layer covered Etna's flanks between Giarre and Catania (from E to S). Ash deposition results from the model were validated against 47 samples collected between 25 and 26 July from measured areas on flat open spaces where the deposit did not show any reworking.The computed tephra dispersal and deposition from NMMB-MONARCH-ASH was able to reproduce the bilobate shape of the real deposit with the two axes oriented toward Acireale and Aci Castello towns.Figure 9 compares the simulated deposit load (kg m −2 ) at the end of the run against the isopachs map derived from measurements taken from 21 to 24 July (Scollo et al., 2007).The model resulted in a cumulative mass of ∼ 1.18 × 10 9 kg.This value is in agreement with the results obtained from Scollo et al. (2007).

Model intercomparison: NMMB-MONARCH-ASH vs. FALL3D
To validate the model performance of NMMB-MONARCH-ASH for its operational implementation, we compare the tephra deposition results of the model against those of the operational FALL3D model for the reconstruction of the 2001 Mount Etna eruption.For this comparison we ran both models using the same meteorological and volcanological initial conditions (Table 8). Figure 10 shows the simulated thicknesses (vertical axis) for both transport models against the observations (horizontal axis) presented in Scollo et al. (2007).The model improved the tephra distribution results from FALL3D simulations for the same event (R 2 ; 0.80/0.62),reducing the RMSE (0.014/0.24) and bias (0.02/0.6) and the computational time by 1 order of magnitude.In particular, all values simulated with NMMB-MONARCH-ASH plot inside the region between 5 and 1/5 (dashed orange line) times the observed mass at each station.The greatest differences perceived against the observations for both models belong to those points located at distances less than 15 km from the vent associated with the uncertainty in the ESPs.The mean value of the relative error between the computed values and observed data is 64 %, which improves those from FALL3d (91 %), and is comparable with that of Scollo et al. (2007), who obtained a 57 % by deposit best fitting using the HAZMAP dispersion model.

Operational forecast with NMMB-MONARCH-ASH
The Barcelona Supercomputing Center is currently working on a modeling integrated system to provide operational forecast of volcanic ash with NMMB-MONARCH-ASH.The system includes a preprocessing tool (prepares the model for real-data simulations), an executable file to run the model, and a user-based postprocessing utility tool.Figure 11 shows a simple schematic representation of the operational implementation of the model.The outcomes of this model- ing system are currently being evaluated against two operational models: (i) the NOAA/ARL Hybrid Single Particle Lagrangian Integrated Trajectory Model (HYSPLIT; Draxler and Hess, 1998) -used at the Washington VACC -and FALL3D (Costa et al., 2006;Folch et al., 2009), used at the Buenos Aires and Darwin VAACs.This section introduces the structure of the operational NMMB-MONARCH-ASH system.Preliminary results for the model intercomparison against FALL3D are described in Sect.4.2.2.

The preprocessing system
The preprocessing utility system consists of a set of programs whose collective role is to prepare the model for real-data simulations.Programs are grouped to preprocess geographical, meteorological, and climatological inputs and interpolate those to the model grid(s).The preprocessing system uses three main programs: runfix, degrib, and runvariable.
-Runfix defines the model domain(s) and interpolates static geographical data to the model grid(s).In addition to computing the latitude and longitude of the rotated grid points, this program interpolates soil categories, land use types, terrain height, annual mean deep soil temperature, monthly albedo, maximum snow albedo, and slope category.
-Degrib extracts the necessary meteorological fields from files formatted with gridded binary (GRIB), used as an initial condition for global simulations and as initial and boundary conditions for single regional domains (i.e., not nested with a global domain).GRIB files contain time-varying meteorological fields obtained from another regional or global NWPM.In addition to the available NCEP's North American Model (NAM) or Global Forecast System (GFS) model, the program has been updated to include European Centre for Medium-Range Weather Forecasts (ECMWF) ERA-Interim reanalysis data (Dee et al., 2011) as forcing.
-Runvariable interpolates the meteorological fields extracted by debgrib to the model grid(s) defined by run-fix and prepares the climatological schemes.This program generates the initial and boundary conditions that are ingested by NMMB using the NOAA Environmental Modeling System (NEMS; Janjic, 2005;Janjic and Black, 2007), a high-performance software superstructure and infrastructure based on the Earth System Modeling Framework (ESMF) for use in operational prediction models at NCEP.

The ash module I/O files
The model takes three run-specific input files, including the model input file (nmmb.inp),which defines the computational and physical schemes needed by the meteorological core, the atmospheric model's fundamental time step, and the parameterization for chemical processes and radiative schemes for aerosol tracers (including ash), amongst other properties of the model.For long-lasting eruptions, the model performs restart runs initializing the tracers from the previous day's history file.
the ash input file (ash.inp), which defines those parameters employed in the ash module.The user-defined parameters include (i) the characterization of the source term: eruption source type, column height and determination of the mass eruption rate, eruption duration, aggregation processes, and particle settling velocity model.In the event of various eruptive phases, the respective ESPs for each phase can be defined; (ii) the settings to turn the gravity current model on or off, altering the particle transport in the umbrella cloud; and (iii) the definition of the coupling strategy (on vs. offline) employed by the model.
the granulometry input file (ash.tgsd), which specifies the diameter, density, sphericity, and relative mass fraction of each particle bin.This information is typically obtained from field data or created by external utility programs for idealized grain size distributions.If aggregation is active, a new bin class for aggregates is added to the granulometry input file.
Once a simulation is concluded, NMMB-MONARCH-ASH writes the following output files: -A log file (ash.log) is written containing information about the run, including a summary of the computed volcanic ash source and mass balance statistics for each time step, and errors and warnings if any.
-A forecast results file (problemname.nc) is written in NetCDF format containing, amongst other variables, the total column mass concentration (g m −2 ) and ground deposition (kg m −2 ) for all bins, the concentration at different flight levels (g m −3 ) and the aerosol optical depth.This information can be processed using several open-source programs to generate plots and animations.Alternatively, the post-process utility program NMMB2GMT has been developed to generate basic General Mapping Tools (GMT) scripts automatically.
-A restart file (nmmb.hst)that is used to initiate a new run using the ash concentrations from a previous simulation is written.

The postprocess system
The postprocess utility tools are designed to interpolate outputs from the NMMB-MONARCH-ASH native grid(s) to National Weather Service (NWS) standard levels (pressure, height, etc.) and standard output grids (Lambert Conformal, polar-stereographic, etc.) in NetCDF format.The system also includes the NMMB2GMT program, which uses the GMT software (Wessel and Smith, 1991) to produce similar plots to the Volcanic Ash Graphics (VAG) used by Volcanic Ash Advisory Centers in operational forecasts.

Scalability analysis
To optimize a future operational implementation of the model, we aim to minimize the time to solution, avoiding communication overhead.In this context, we evaluate the model scalability (scaling efficiency) for its regional and global configurations by performing a strong scalability test, in which the problem size of our simulation (e.g., model domain and resolution) remains fixed while increasing the number of processing cores.Figure 12 shows the parallel speedup (S, Eq. 19) and efficiency (E, Eq. 19) of the modeling system for a global simulation of the climactic phase for the 2011 Cordón Caulle (Table 7).On the MareNostrum-III supercomputer, maximum efficiency for the global simulation described in Table 7 is reached between 32 and 40 nodes (16 CPUs each), with a parallel efficiency of 0.6.The scalability analysis was performed on all the available source term and sedimentation schemes in the model.The relative computational cost associated with the main processes in NMMB-MONARCH-ASH is presented in Fig. 13.Processes represented include meteorological prediction, volcanic ash transport and sedimentation forecast, aggregation of particles, gravity current effects, and the restart phase.The restart phase represents the CPU time employed to rerun the preprocess system every 24 h of simulation.This figure suggests that the computational increase (CPU time) associated with the ash module can vary from 5 to 55 %, depending on the number of computational nodes employed.It is important to note that, depending on the settling velocity model employed, up to 60 % of the time allocated to the ash module is assigned to the sedimentation term.
Results from the scalability analysis show that the model performance (in terms of speed-up) depends on the problem size as well as on the domain partitioning topology.In that context, the relative computational cost of the model's me- teorological core (NMMB) is evaluated as a function of its domain decomposition (e.g., distribution of processing units for the horizontal domains -nodes i and j ).For this analysis the bin-performance dependency of the model is considered, therefore evaluating only the cost of one bin of ash.Results from this analysis suggest that, for an optimal simulation using 32 nodes, the computational cost of the meteorological core decreases over 10 % when the weight of the decomposition is focused on the j nodes (e.g., more computational resources assigned for the fast Fourier transformation algorithm).The best domain decomposition resulted in 6(i) × 84(j ) + 8(w), where i and j are the number of processors employed in the horizontal and vertical domains respectively, and w is the number of writing processors.
For operational purposes, the computational time employed to provide ash dispersal forecast using NMMB-MONARCH-ASH is evaluated for the global simulation with 1 bin of ash.The maximum time required by the model to perform a 24 h forecast, running all the available processes (e.g., advection, diffusion, sedimentation) every time step (180 s) is less than 3 min when using the best domain decomposition presented before (6 × 84 + 8).This time can be further optimized for operational purposes, i.e., calling the model physics less frequently in order to save computational time.As a general rule of thumb, the adjustment time step in seconds for the meteorological core can be taken as 2.25 times the grid spacing in kilometers.For higherresolution model runs made without parameterized convection, a time step in seconds of about 1.9 to 2.0 times the grid spacing may be more appropriate (Janjic and Gall, 2012).

Cost-benefit analysis
Employing online models for operational dispersal forecast requires larger computational resources and is not always feasible at all operational institutes.Nevertheless, due to the increase in computing power of modern systems, one can argue that such gradual migration towards stronger online coupling of NWPMs with TDMs poses a challenging but attrac-tive perspective from the scientific point of view for the sake of both high-quality meteorological and volcanic ash forecasting.
The focus on volcanic aerosols integrated systems in operational forecast is timely.Experiences from other communities (e.g., air quality) have shown the benefits from twoway online meteorology-chemistry modeling.For example, the importance of the different feedback mechanisms for meteorological and atmospheric composition processes have been previously discussed for models developed in the USA (Zhang, 2008) and Europe (Baklanov et al., 2014).These benefits have been recently stressed by several studies covering the analysis of the aerosol-transport and aerosol-radiation feedbacks onto meteorology from the air quality model evaluation international initiative (AQMEII) in its phase 2 (Alapaty et al., 2012;Galmarini et al., 2015) and the EuMetChem COST Action ES1004 (EuMetChem, http://eumetchem.info)Demonstrating these benefits, however, requires running the online model with and without feedbacks over extended periods of time.For the particular case of volcanic aerosols, further research is still required to quantify the benefits posed by online couple models over traditional offline TTDM on both atmospheric transport and the radiative budget.The Barcelona Supercomputing Center is currently working to quantify these benefits with the NMMB-MONARCH-ASH model and assess how the magnitude of the model forecast errors implicit in the offline approach compares with other better-constrained sources of forecast error, e.g., uncertainties in eruption source parameters.Preliminary results from this study indicate that meteorology-transport inconsistencies from offline models can be, in some cases, on the same order of magnitude as those associated with the eruption source parameters.In terms of computational cost, the computational efficiency of the model's meteorological core allows for online integrated operational forecasts that employ an equivalent computational time as FALL3D for the same computational domain and number of processing cores.
We present NMMB-MONARCH-ASH, a new online multiscale meteorological and transport model developed at the Barcelona Supercomputing Center (BSC) to forecast the dispersal and deposition of volcanic aerosols.The objective of this model is to improve the current state-of-the-art tephra dispersal models, especially in situations where meteorological conditions fluctuate rapidly in time, two-way feedbacks are significant, or long-range ash cloud dispersal predictions are necessary.The model predicts ash cloud trajectories, concentration of ash at relevant flight levels, and the expected deposit thickness for both regional and global domains.NMMB-MONARCH-ASH solves the mass balance equation for volcanic ash by means of a new ash module embedded in the BSC's operational system for short-midterm chemical weather forecast (NMMB-MONARCH).In addition to volcanic ash, the system is also capable of forecasting the dispersion of other atmospheric aerosols (e.g., dust, sea salt, black carbon, organic aerosol, sulfates).Its multi-scale capability allows for nested global-regional atmospheric transport simulations, taking into account the characterization of the source term (emissions), the transport of volcanic particles (advection-diffusion), and the particle removal mechanisms (sedimentation-deposition).The model has been shown to be robust and scalable to arbitrary domain sizes (regional to global) and numbers of processors.
The forecast skills of NMMB-MONARCH-ASH have been validated against several well-characterized eruptions, including the Etna 2001 (Italy), Chaitén 2008 (Chile), Cordón Caulle 2011 (Chile), and Pinatubo 1991 (Philippines) eruptions (e.g., Marti et al., 2013Marti et al., , 2014)).To evaluate the online coupling strategy and the multi-scale capability of the model, this paper summarizes the regional and global configurations of the model to forecast the dispersal of ash for the first days of the 2011 Cordón Caulle eruption (strong long-lasting eruption with rapid wind changes).In addition, to evaluate the sedimentation mechanisms of the model, this work also includes the results from the regional configuration of the model for the first phase of the 2001 Etna eruption, a good case study of a weak long-lasting eruption with well-characterized tephra deposits.Simulation results demonstrate that NMMB-MONARCH-ASH is capable of reproducing the spatial and temporal dispersal variability of the ash cloud and tephra deposits.

Figure 1 .
Figure 1.Schematic representation of the main components of an atmospheric transport model.Red text shows model specifications for the transport of volcanic ash.
7) α m = mean sticky efficiency decay per unit +A DS θ 4/D f n n tot = number of particles available to aggregate Number of For Brownian motion: A B = − 4k b T 3µ 0 k b = is the Boltzmann constant 1.38 × 10 −23 m 2 kg s −2 K particles available T = absolute temperature to aggregate µ 0 = dynamic viscosity of air Kernels Ambient fluid shear: A S = 2 S γ 4 3 S = fluid shear Differential sedimentation: γ = relationship of particle diameter to volume fractal radial velocity as a function of time () or cloud radius (R) = empirical constant (≈0.2) (Suzuki and Koyaguchi, 2009) = Brunt-Väisälä frequency (atm.ambient stratification) = volumetric flow rate into the umbrella

Figure 2 .
Figure2.Meteorological synoptic situation during the first 2 weeks (4-14 June) of the 2011 Caulle (white star) activity over South America.Plots show the direction (vector) and velocity (contours m s −1 ) of the wind at 9100 m a.g.l.(above ground level) (300 hPa circa).Meteorological data obtained from the NMMB meteorological forecast driven with ERA-Interim reanalysis at 0.75 • horizontal resolution.

Figure 3 .
Figure 3. Composite image of NMMB-MONARCH-ASH results for dispersion of ash for the 2011 Caulle eruption at different time slices.Simulation results are compared against split window algorithm NOAA-AVHRR satellite images (bands 11-12 microns).Contours indicate ash column load (g m −2 ) resulting from integrating the mass of the ash cloud along the atmospheric vertical levels.

Figure 5 .
Figure 5. Left panels: isopach maps in centimeter of layers A-B, A-F, H, and K2 from the 2011 Caulle eruption.Dashed lines infer the limit of the deposits presented in Pistolesi et al. (2015).Right panels: corresponding NMMB-MONARCH-ASH computed thicknesses (centimeters).Key locations in blue include San Carlos de Bariloche (SCB) and Ingeniero Jacobacci (IJ), 90 and 240 km east of the vent.

Figure 6 .
Figure 6.Left panel: measured ground deposit isopachs (kg m −2 ) for the 2011 Caulle eruption between 4 and 30 June 2011.Dashed lines infer the limit of the deposits (modified from Collini et al., 2013).Right panel: predicted deposit load (kg m −2 ) with NMMB-MONARCH-ASH at the end of the simulation.Key locations in blue include San Carlos de Bariloche (SCB; 90 km from the vent), Ingeniero Jacobacci (IJ; 240 km east of the vent), and Trelew and Viedma (∼ 600 km SE and NE of the vent, respectively).

Figure 9 .
Figure 9. Left panel: isomass map of the tephra deposit for the 2001 Mount Etna eruption formed between 21 and 24 July.Curves are given in kg m − 2. Coordinates are given in UTM-Datum ED50 (Scollo et al., 2007).Right panel: modeled deposit load (kg m −2 ) with NMMB-MONARCH-ASH at the end of the event.

Figure 10 .
Figure 10.Simulated versus observed thicknesses for the reconstruction of the 2001 Etna eruption with NMMB-MONARCH-ASH (circles) and FALL3D (crosses).The solid bold line represents a perfect agreement, while the dashed and solid thin orange lines mark the region that is different from observed thicknesses by a factor of 5 (1/5) and 10 (1/10), respectively.

Figure 11 .
Figure 11.Schematic representation of the operational implementation of NMMB-MONARCH-ASH.

Figure 12 .
Figure 12.NMMB-MONARCH-ASH scalability results.Top panel: parallel speed-up (S, computational speed) for meteorology only (blue) and for meteorology and dispersal combined (red).Bottom panel: parallel efficiency (E) vs. number of computation nodes employed.

Table 1 .
Main characteristics of the NMMB-MONARCH-ASH meteorological solver.

Table 2 .
Options implemented in NMMB-MONARCH-ASH to estimate the mass eruption rate from column height.Unless otherwise noted, the units for all parameters are in SI.

Table 3 .
Ash aggregation options in NMMB-MONARCH-ASH from analytical solutions based on field observations.Default aggregate properties can be modified by the user.n/a = not applicable.
k = number of particles of diameter k in an aggregate Number of

Table 7 .
Model configuration for the 2011 Cordón Caulle regional and global runs.The regional run used a horizontal resolution of 0.15 • × 0.15 • with a 30 s dynamic time step, while the global domain used a horizontal resolution of 1 • × 0.75 • with a 180 s dynamic time step.

Table 8 .
Model configuration for the 2001 Mount Etna regional simulations.Regional Run1 used a horizontal resolution of 0.1 • × 0.1 • with a 30 s dynamic time step, while Run2 used a finer horizontal resolution of 0.05 • × 0.05 • with a 10 s dynamic time step.