Simulations and parameterisation of shallow volcanic plumes of Piton de la Fournaise , La Réunion Island , using Méso-NH version 4-9-3

Abstract. In mesoscale models (resolution ~ 1 km) used for regional dispersion of pollution plumes the volcanic heat sources and emissions of gases and aerosols, as well as the induced atmospheric convective motions, are all sub-grid-scale processes (mostly true for weak effusive eruptions) which need to be parameterised. We propose a modified formulation of the EDMF scheme (eddy diffusivity/mass flux) proposed by Pergaud et al. (2009) which is based on a single sub-grid updraft model. It is used to represent volcano induced updrafts tested for a case study of the January 2010 summit eruption of Piton de la Fournaise (PdF) volcano. The validation of this modified formulation using a reference large eddy simulation (LES) focuses on the ability of the model to transport tracer concentrations up to 1–2 km above the ground in the lower troposphere as is the case of majority of PdF eruptions. The modelled volcanic plume agrees reasonably with the profiles of SO2 (sulfur dioxide) tracer concentrations and specific humidity found from the reference LES. Sensitivity tests performed for the modified formulation of the EDMF scheme emphasise the sensitivity of the parameterisation to ambient fresh air entrainment at the plume base.


Introduction
A critical factor in successfully monitoring and forecasting volcanic ash and gases dispersion is the height reached by eruption clouds, which is mainly controlled by the eruptive mass flux (e.g.Kaminski et al., 2011) but is also af-fected by environmental factors, such as wind shear and atmospheric vertical stability (Glaze and Baloga, 1996;Graf et al., 1999;Bursik, 2001;Tupper et al., 2009).The term "volcanic plume" refers to both the vertical buoyant column of gas/ash above the eruptive vent and the following horizontal transport of pollutants at the regional to hemispheric scales by the wind flow.Therefore, there is a need of numerical prediction systems coupling volcanic plume dynamics and atmospheric circulation models.An attempt of such a system was proposed by Kaminski et al. (2011) for the deep tropospheric 2010 eruption of the Eyjafjallajökull volcano in Iceland.
The convective scale of a volcanic plume corresponds to the unstable region where intense but localised sensible and latent heat fluxes released by pyroclasts, gases and lava near eruptive vents generate convection which transports energy and pollutants to higher altitudes through buoyant plumes.Throughout the course of this convection, mixing of the plume with the atmosphere takes place at different levels of altitude through entrainment and detrainment.This process allows for the distribution of pollutants over a certain vertical range.
Piton de la Fournaise (PdF) is one of the world's most active volcanoes (Lenat and Bachelery, 1987) with an average of one eruption every 8 months in the last 50 years (Peltier et al., 2009).Most of the studies undertaken for deep volcanic injection have been applied to stratospheric injections, which are mostly performed by large explosive volcanoes (comprehensive review by Robock, 2000).However, much less is known about the environmental and atmospheric impacts and fates of weak volcanic plumes injected into the troposphere (Mather et al., 2003;Delmelle et al., 2002).PdF can create a major source of tropospheric air pollution as was the case during the eruption of April 2007 (Tulet and Villeneuve, 2011).Details on the island areas affected by pollution during this eruption can be found in Viane et al. (2009) and Bhugwant et al. (2009).The air-quality standard for ecosystem and human health protection was exceeded for sulfur dioxide (SO 2 ) at several inhabited locations of the island (Bhugwant et al., 2009).Suzuki et al. (2005) developed a three-dimensional numerical fluid-dynamics model to explicitly simulate volcanic plumes and explore different dynamical regimes as function of the ejection velocity and the mass discharge rate of the volcanic material.However, the spatial resolution is very fine in their model, with a horizontal grid spacing well below the ones used in pollution dispersion models at regional scale (at best 1 km).Such a model with presumably high numerical cost is thus not applicable for air-quality prediction purposes.
Simulations of atmospheric plumes from intense heat source points have been performed using the Méso-NH (Lafore et al., 1998) model to represent the impact of forest fires on the dynamics and chemistry of the atmosphere.A study by Strada et al. (2012) simulated forest fire plumes at 1 km resolution which showed good agreement with observations where high sensitivity to the atmospheric stability was observed.
Simulations of buoyant eruptive columns, chemistry dispersal in the proximal environment and the volcanic cloud tracking at regional scale can be based on similar numerical and conceptual approaches as the ones used for the study of forest fire plumes.However, volcanic eruptive vents usually cover small areas and in (at best) kilometric-resolution models used for air-quality purposes (simulation or forecasts), the localised heat source is diluted in the model grid; hence, no convection is explicitly generated.
Several types of atmospheric movements are sub-grid processes, and they are incorporated into atmospheric models through appropriate parameterisation schemes.In order to determine the evolution of volcanic plumes in the atmosphere, numerical models need to consider two different scales: 1. an implicit/convective scale corresponding to the convective plume above the erupting volcano, whose processes are sub-grid even at fine resolutions (> 500 m), and 2. an explicit/dispersion scale that corresponds to the dispersion of the volcanic plume in the atmosphere.
In mesoscale models used for regional dispersion of pollution plumes (target resolution ∼ 1 km), the volcanic heat sources and emissions of gases and aerosols as well as the induced atmospheric convective motions are all sub-grid-scale processes which need to be parameterised.In this article we first briefly describe an existing sub-grid shallow convection scheme by Pergaud et al. (2009) used in the atmospheric model Méso-NH for conventional weather simulations.This scheme is based on a single sub-grid convective updraft approach, whereby the updraft vertical development is calculated step by step from one vertical model level to the level above.Therefore, the updraft needs to be initialised at the ground level, and this relies on local atmospheric turbulence in the scheme formulation as per Pergaud et al. (2009).In Sects.2.2.2 and 2.2.3, we propose a specific adaptation of this scheme whereby the size and intensity of the volcanic heat source serve as alternative initial conditions at ground level for the modelled updraft.
Due to the computational efficiency of a one-dimensional (1-D) model and the ability to isolate a column of atmosphere for study, 1-D modelling is an ideal configuration to develop and test parameterisations (Randall et al., 1996).Due to specific constraints (see Sect. 2.3.4), the new parameterisation is tested for an eruption observed at PdF in January 2010 not really in a 1-D model but actually in the central column of a 3 × 3 column model.This central column can be seen as a quasi-1-D model with open lateral boundary conditions and thereafter referred to as SCM (single column model).The choice of 1 km as horizontal resolution for SCM simulations is because it is the target resolution of future forecast models running over Réunion Island.
Simultaneously, a three-dimensional (3-D) large eddy simulation (LES) is performed (10 m resolution) using the same size and intensity of the eruption as prescribed for the SCM simulations with adapted convection scheme.Observations relating to the case study are used to evaluate the LES simulation which is further used as reference to validate the SCM results.As outlined by Pergaud et al. (2009), both Siebesma et al. (2003) and Brown et al. (2002) have shown that LES are robust for representing shallow cumulus convection.This methodology has been used in the past by the Global Energy and Water Cycle Experiment Cloud System Study (GCSS) (Browning, 1993) and also applied to test the scheme used to parameterise shallow convection (Pergaud et al., 2009).

January 2010 summit eruption of Piton de la Fournaise
An eruption took place on 2 January 2010 around 10:20 UTC at the summit of PdF located at 2632 m a.s.l. as detected by the monitoring networks of the Piton de la Fournaise Volcanological Observatory (OVPF/IPGP) (Roult et al., 2012).At 10:27 UTC a small and diluted gas plume was first visible and a vertical plume rapidly formed above the crater at 10:57 UTC.Up to seven lava fountains erupted together from the same number of vents along a fracture on the west Dolomieu crater wall (length of fracture about 60 m) and the highest lava fountain of about 30 m was emitted by the largest vent located in the middle of the fracture.Lava flows were fed by magma flowing from the vents (Fig. 1) but also from hot fountain products that were remobilised after falling to the ground.According to Roult et al. (2012) the eruption emitted 1.2 × 10 6 m 3 of lava in about 9.6 days; the mass flow rate decreased exponentially after the beginning of the eruption and the fountaining and gas plumes described here only occurred till 4 January 2010, whereafter mostly effusive lava flows were observed.The vertical plume above the crater (Fig. 1; right) was relatively steady implying low winds and a level of neutral buoyancy was reached at approximately 1300 m above the fountain top (from observation).For the development of a parameterisation this case study is the least complex one compared to other eruptions of PdF since 2000.There was a welldeveloped vertical gas column which was weakly affected by horizontal winds and the topography of the area.

Description of the volcanic plume parameterisation
It is well understood that a volcanic eruption plume enters into an atmosphere that has a pre-existing stratification in terms of temperature, moisture content and wind (Bjornsson et al., 2011;Petersen et al., 2012).There are three dynamically distinct regions related to volcanic plumes (Sparks, 1986): 1. the gas thrust region, where the dynamics is dominated by the exit velocity at the vent and the flow near the vent is driven upward by its initial kinetic momentum; 2. buoyancy-driven convective region which covers most of the height of the plume; and 3. umbrella cloud region, where vertical motion is small and the plume disperses horizontally due to wind impacts.
For the purpose of modelling volcanic clouds using Méso-NH, we are predominately interested in the convective region of the volcanic cloud.For the kind of effusive eruption under consideration in this study, the gas thrust region extends only over few metres (see, e.g.Fig. 2 for an eruption comparable to January 2010).For simplicity, it will be assumed that the thrust region is very short compared to the total vertical extension of the plume and that the plume is primarily driven by buoyancy.Thus, the plume will be assumed to be convective from the ground level to its top.
The current updraft model used in Méso-NH defined by Pergaud et al. (2009), Sect.2.2.1, is not adapted to volcanic plumes.We propose an adaptation of the updraft scheme (Sects.2.2.2 and 2.2.3) which is applied to volcanic plumes and consists mainly in a modification of the updraft initialisation at ground level (z grd ) using values inspired from terrain observations.

Sub-grid cloud parameterisation as per Pergaud et al. (2009)
The basic idea of the EDMF (eddy diffusivity/mass flux) approach is to represent vertical transport of matter and energy that occurs at the sub-grid scale in numerical simulations of the convective boundary layer (CBL) with resolutions of written in the form of −K φ ∂φ ∂z , where −K φ is a diffusion coefficient and φ the average of any model variable φ (temperature, tracer mixing ratio, etc.) over the local grid cell (Holton, 2004).
A grid box can contain multiple convective updrafts.For simplicity a single updraft is considered carrying the properties of the ensemble of updrafts.This is known as the massflux (MF) approach.The fraction of the total area of a grid box that is covered by the updraft is known as the fractional updraft area (a u ).The corresponding net vertical flux for φ over the grid cell takes the form of M u ρ (φ u − φ), where M u is the updraft mass flux, φ is the mean value and φ u is the updraft value of the variable φ.
Both ED and MF approaches have been combined in a single EDMF parameterisation such that nonlocal sub-grid transport due to strong updrafts is taken into account by MF, while the remaining transport is taken into account by ED (Siebesma and Teixeira, 2000;Hourdin et al., 2002;Soares et al., 2004;Siebesma et al., 2007;Pergaud et al., 2009;Witek et al., 2011).In our approach for volcano-induced convection we only modify the MF scheme (Sect.2.2.2).
The two key parameters determining the mass-flux profile are entrainment (ε) and detrainment (δ), expressed as fractions of the updraft mass flux (M u ) per unit height.This simply leads to the following steady state mass-flux continuity equation: The mass-flux evolves along the vertical at a rate given by the difference between the ε and δ rates.The definition of entrainment/detrainment rates is the crucial point in EDMF parameterisation as it is at this level that the physical coupling between turbulent mixing and mass flux is performed.
In Pergaud et al. (2009) the mass-flux profile depends on the vertical velocity of the updraft (w u ), whose vertical evolution is affected in turn by a buoyancy force (B u ), and a drag term where the entrainment of environmental air, namely lateral mixing, is accounted for: The updraft buoyancy acceleration is evaluated in relation to the difference of virtual potential temperature (θ v ) between the updraft and its environment: B u = g(θ u, v −θ v )/θ v ; coefficients c 1 and c 2 are usually set to 1 (Simpson and Wiggert, 1969).Independent solutions of Eqs. ( 1) and (2) permit calculating the vertical variation of the updraft fractional area, that is used to diagnose the cloud fraction, hence defining the sub-grid condensation scheme in the EDMF framework.

Modified EDMF -updraft initialisation
Firstly, in the current EDMF parameterisation w u is initialised at the ground level (z grd ) using turbulent kinetic energy (TKE) (e) as w 2 u (z grd ) = 2 3 e(z grd ), which is bound to local meteorology.However, this computation is not applicable to volcanic plumes as vertical velocity in this case does not depend on the atmosphere through the TKE.During volcanic eruptions, a mixture of gases, magma fragments, crystals and eroded rocks is injected into the atmosphere at high velocity, pressure and temperature.The diverse and unpredictable variability of eruptive styles depends mostly on the complex rheology of magma and the nonlinear processes leading to the fragmentation of the viscous melt into a mixture of gases and particles (Gonnermann and Manga, 2007).Nonetheless, the explosive character of a magmatic eruption like that of January 2010 is associated with the rapid decompression and the consequent abrupt expansion of gases in the magma (Parfitt and Wilson, 2008).In order to simplify, we consider the vertical velocity of the updraft w u (z grd ) as the vertical velocity of the lava fountain (a variable that is mostly known from observation).The input data mentioned in this section (used for updraft initialisation) and the following sections are listed in Table 1.
Secondly, the updraft fraction area is simply initialised as the ratio of the fissure surface (S Fis,SCM ) by the model cell surface (S MNH ).
Now, as w u (z grd ) and a u (z grd ) are both known and are independent of one another, using a similar principle as in Pergaud et al. (2009), the mass flux at the ground can be calculated such that M u (z grd ) = ρ mix (z grd ) a u (z grd ) w u (z grd ). (5) The ground level density of the updraft, ρ mix (z grd ), is approximated by a mixture of the two main gases at PdF (H 2 O and SO 2 ) considered as perfect gases, such that ρ mix (z grd ) = P (z grd ) T u (z grd )R mix , where P (z grd ) is the ambient pressure at ground level, T u (z grd ) is the temperature of the updraft at ground level and R mix represents the specific gas constant of the mixture, composed mostly of water vapour and SO 2 : ), where R is the universal gas constant; and M H 2 O and M SO 2 , are the mixing mass ratios and molar mass of water vapour and SO 2 , respectively.Equation (5) uses ρ mix rather than using density of dry ambient air (as in the standard formulation from Pergaud et al., 2009, our Eq. 3).Indeed, in magmas like those erupted in 2010 the gas melange is dominated by water vapour, i.e. about 80 % of the melange mass (Di Muro et al., 2014) and the remaining 20 % is that of SO 2 (i.e.[H 2 O] = 0.8 and [SO 2 ] = 0.2 kg kg −1 ).This gives a H 2 O / SO 2 ratio of  4, which is the ratio expected by simple closed system degassing of PdF shallow magmas.This value is at the lower end of the range actually measured by the OVPF geochemical network (Allard et al., 2011).Note also that ρ mix (z grd ) is formulated as a perfect gas mixture, which implicitly assumes that no solid fraction is present in the atmospheric plume.This is a reasonable assumption for such an eruption, whereby volcanic ash represents a small fraction of the volcanic plume.

Modified EDMF -basal lateral mass exchange
Entrainment of ambient air through turbulent mixing plays a central role in the dynamics of eruption plumes, primarily because the plume density is controlled by the mixing ratio between ejected gas/material and ambient air (Woods, 1988).Furthermore, the amount of air entrained controls the heights of eruption columns (Suzuki and Koyaguchi, 2010).In the current EDMF (Sect.2.2.1); the mass-flux entrainment of the updraft ε at the ground level has a constant value of 0.02 m −1 whereas δ is zero.
In this sub-section we present the modifications to the input method of ε and δ such that, for some height z above the ground, a desired mass of ambient air may be entrained into the updraft and conversely a desired mass of the updraft may be expelled.Above this height ε and δ are both calculated as defined by Pergaud et al. (2009) and the coexistence of entrainment/detrainment both continue to feed the vertical evolution of M u .
The importance of adjusting the ground level ε and δ will become more apparent in Sect.3.However, the modifications are presented below.Figure 3 assembles all modifications made to the EDMF model along with the input variables (marked in red) used at ground level.
Let M env represent the mass flux of environmental air that enters the updraft between levels z grd and z grd + z.Hence updraft mass flux at (z grd + z) is simply defined as Let α = rearranging Eq. ( 6), If ε and δ are constant between z grd and (z grd + z), then, by integrating Eq. ( 1) between z grd and (z grd + z), Eq. ( 7) can be rewritten as Finally, using Eqs.( 7) and ( 8), For a desired fraction α of ambient air entrained in the volcanic gas column at z grd + z, an infinity of entrainment and detrainment rate combinations can be prescribed such that Eq. ( 9) is respected.

Simulation set-up and configuration
For our chosen case study, three sets of simulations were run as depicted in Fig. 4.
1. Section 2.3.2 describes the 3-D spin-up simulation which is used to generate background atmospheric profiles used for both the LES and SCM simulations described below.
2. Section 2.3.3 details the LES considered as the reference.

Common features to all simulations
The Méso-NH model (version MNH-4-9-3) is used in this study; it is a mesoscale non-hydrostatic atmospheric model able to simulate convective motion and flow over sharp topography.This model has been jointly developed by Laboratoire d'Aérologie (UMR5560 UPS/CNRS) and Centre National de Recherches Météorologiques -Groupe d'études de l'Atmosphère Météorologique, CNRM-GAME (UMR3589 CNRS/Météo-France), and is designed to simulate atmospheric circulations from small-scale (type -LES) to synoptic-scale phenomena (Lafore et al., 1998).All Méso-NH related documentation and articles along with various model versions are available at http://mesonh.aero.obs-mip.fr.Different sets of parameterisations have been introduced for cloud microphysics (Cohard and Pinty, 2000), turbulence (Bougeault and Lacarrere, 1989) and convection (Bechtold et al., 2001).The shallow convection in Méso-NH is parameterised according to Pergaud et al. (2009) while for the purposes of this study no deep convection parameterisation was activated.The ISBA (Interactions Soil-Biosphere-Atmosphere scheme) (Noilhan and Mafhaf, 1996) is the scheme used for land surfaces in order to parameterise exchanges between the atmosphere and the ground providing surface matter and energy fluxes to the atmosphere.The turbulent scheme implemented in Méso-NH is a full 3-D scheme that has been developed by Cuxart et al. (2000) with regards to both LES and mesoscale simulations.Kessler's warm microphysical scheme (Kessler, 1969) was activated during the simulation.Méso-NH can be used for idealised as well as real case studies and for the purpose of this article we focus on idealised case studies.For all simulations performed, a vertical grid composed of 72 levels in the Gal-Chen and Sommerville (1975) coordinates is used, with a vertical mesh stretched from 40 m at the ground to 600 m at the model top.

3-D spin-up simulation to generate background profiles
A 3-D spin-up simulation is performed to generate the background profiles which are used for SCM and LES.Two, twoway grid-nested domains with horizontal mesh sizes of 4 and 1 km are used (Fig. 4a).Both domains have 100 points in x and y.The initial state for the simulation, as well as the boundary conditions updated every 6 h for the outermost model, is provided by analyses from the French operations forecasting system for the Indian Ocean, ALADIN-Réunion (9.6 km resolution; Montroty at al., 2008).The simulation starts 1 January 2010 at 00:00 UTC and ends 2 January 2010 at 18:00 UTC using a time step of 1 and 0.25 s for the 4 and 1 km resolution models, respectively.
Figure 5 shows the vertical profiles of temperature ( • C), potential temperature (K) and water vapour mixing ratio (g kg −1 ) as simulated at 10:50 UTC on 2 January 2010 by the spin-up simulation for the local area of interest (location of the PdF volcano).The ambient atmosphere is dry with water vapour concentration just under 8 g kg −1 at the ground and decreasing with altitude.The tropopause is found at about 16 km above ground level (a.g.l.hereafter, where the ground level corresponds to about 2.6 km above the sea) which corresponds well to tropical climates.The 0 • C isotherm is located at 2.7 km a.g.l.Those profiles are then used as initial and steady background conditions for our LES and SCM simulations.
The vertical structure of trade winds over Réunion Island was investigated by Lesouef (2010) and Lesouef et al. (2011).The trade wind inversion located at about 4 km a.s.l.(Taupin et al., 1999) is described as a consequence of the descending branch of the Hadley cell circulation (Lesouef et al., 2011) where easterly winds prevail in the lower levels while westerly winds prevail in upper levels.It coincides with a temperature inversion or at least a layer of enhanced vertical static stability.This is found in Fig. 5 (middle) at about 2 km a.g.l.(4.6 km a.s.l.) as an increased gradient of potential temperature.This stable layer can behave as a bar-rier for development of clouds (Hastenrath, 1991) but also for plumes generated through our simulations.
It should be noted that the wind profiles as obtained from the spin-up simulation appeared to be unrealistic since the wind near the ground (7-8 m s −1 ) was clearly overestimated; as a consequence, a strong tilt in the volcanic plume above the crater was simulated in a first tentative LES using these wind profiles (not shown) but clearly not observed in reality (Fig. 1).In addition, no strong wind above the caldera rim of Bellecombe has been reported by Météo France for the simulation period (average 10 m wind of 2.5 ± 0.9 (1σ ) m s −1 , with a maximum hourly mean of 4.9 m s −1 ).For this reason, instead of using wind profiles from the spin-up simulation, we prescribed a vertically uniform and very slow wind field (u(z) = 0.1 m s −1 and v(z) = 0 m s −1 ) as background wind in the LES and for consistency also in the SCM.

LES simulations
An LES model has such a high resolution that it can resolve not only convective motions but also the largest eddies (responsible for the major part of the turbulent transport).This section describes the set-up of the LES simulation considered as reference used to validate the EDMF parameterisation for volcano-induced convection.
Table 2 summarises the configuration of the LES model.Its horizontal physical domain is chosen to extend over 1 km × 1 km, which exactly corresponds to the size of the central column in our quasi-1-D model (see Sect. 2.3.4).Thus, horizontal averages of the LES fields will be taken as references to validate the quasi-1-D model output profiles.The LES horizontal resolution is 10 m × 10 m, such that convection can be explicitly resolved.Due to the short simulation duration, radiative processes are neglected.As the model domain is quite small, and also to avoid complex topographic effects, the local topography of the volcano is not taken into account (except the fact that the model ground is at the correct altitude, such that the ground pressure is 78 695 Pa), depicting a flat domain for simplifying the model (as also done for the SCM model detailed in Sect.2.3.4).).Altitude displayed is above ground level (a.g.l.), which is at 2600 m a.s.l.The surface mass and heat fluxes representing the volcanic mass and energy source in the LES are prescribed for one single surface cell (i.e. S Fis,LES = 100 m 2 ; Fig. 4c) with a correction factor Corr = 1.2 such that the input fluxes are consistent with that of the SCM model, where the volcanic fissure covers an area S Fis,SCM = 120 m 2 .
Let F H 2 O be the vapour mass flux (kg m −2 s −1 ) prescribed for this particular surface cell at every model time step, then and similarly for the SO 2 mass flux (kg m −2 s −1 ): The variables here (all at the ground level) have the same definitions and values as in Sect.2.2.2 (see also Table 1).Finally, let F s represent the sensible heat flux (W m −2 ), then where C p, mix is the specific heat capacity of the mixture at constant pressure such that C p, mix = 4R mix (H 2 O and SO 2 being both triatomic gases); T u and T are the temperatures of the updraft and of ambient air outside the updraft, respectively.An appendix at the end of the text presents detailed derivations of Eqs. ( 10)-( 12).Steady surface fluxes are used as volcanic input in LES runs.Their values are summarised in Table 1.

SCM simulations
Table 3 shows the configuration of the SCM model.The volcanic updraft is simulated only in a single central grid column of size 1 km × 1 km; however, the total number of grid columns used is 3 × 3 (Fig. 4b).This is simply to allow for the use of open lateral boundary conditions and hence avoid matter and energy to accumulate in the model.The adapted EDMF model in Sect.2.2.2 is used to run SCM simulations.The variables used as volcanic input in SCM runs were presented in Sect.2.2.2 and are summarised in Table 1 along with their values.As mentioned earlier, since the gas melange in the eruption column consists of 80 % of H 2 O and 20 % of SO 2 , the SCM model is simply initialised with [H 2 O] = 0.8 kg kg −1 and [SO 2 ] = 0.2 kg kg −1 in the updraft at ground level.
Common to both LES (Sect.2.3.3) and SCM runs, -the wind profiles obtained from the spin-up simulation are not used as background conditions, instead a prescribed uniform wind field is used (u = 0.1 and v = 0 m s −1 ; see Sect.2.3.2); -radiative processes are neglected.
Finally, the comparison method of LES and SCM simulation results is as follows: fields horizontally averaged over the full LES simulation domain (1 km × 1 km) provide vertical profiles, which are compared to vertical profiles from the central SCM grid column of 1 km × 1 km.This is sketched in Fig. 4b and c.

Results and analysis
In this section, results obtained from the 1-D SCM and 3-D LES of the case study are presented and analysed.

Demonstration of the need of specific heat source to generate deep plumes
A first most obvious question is whether we need to parameterise volcanic updraft.Figure 6 shows results from four simulations: Fig. 6a and b shows simulation results from LES without and with volcanic heat sources, respectively, whereas Fig. 6c and d show results from the SCM model without and with volcanic heat source, respectively.Results for Fig. 6b follow the initialisation of the volcanic heat source as outlined in Sect.2.3.3 above and results from Fig. 6d follow the initialisation of the volcanic heat source as outlined in Sect.2.2.2.All four simulations have been initialised with a passive SO 2 tracer as outlined in Table 1 and used as a tracer pollutant injected into the atmosphere.
In simulations with no volcanic heat source, SO 2 tracer is simply diffused to a few hundreds of metres above the ground and the majority of the tracer remains at low altitude (Fig. 6a, c).Results from the reference LES simulation with volcanic source (Fig. 6b) shows an uplift of tracer to higher altitudes, with maximum concentration around 1.0 km above the ground after 90 min, and almost no tracer above 2 km.The SCM simulation with modified EDMF (M.EDMF) results (Fig. 6d) also show tracer lifted to much higher altitudes with the majority of the concentration levelled off at around 7.25 km.The overall tracer concentrations are vertically distributed between 4 and 11 km above the ground.It is clear that without modifications to EDMF and without initialising the LES simulation with specific volcanic heat sources, the two models are not capable to transport tracer concentrations to higher altitudes.
Although both Fig. 6b and d show transport of tracer to higher altitudes, it is evident that, in terms of maximum detrainment height of the tracer (1.4 and 7.25 km respectively) and its vertical profile, at this stage the M.EDMF results are poorly comparable to that of the LES (the reference simulation) -the plume generated by M.EDMF being much too deep.
Hereafter, the height at which there is a maximum detrainment of the tracer will be referred to as the "maximum injec- tion height".The sensitivity of the latter against entrainment and detrainment at the base of the updraft will be investigated, with the aim at obtaining better agreement between the reference LES and M.EDMF simulations.

Influence of entrainment/detrainment at the base of the updraft
It is well known that both entrainment and detrainment have an impact on the updraft development because they affect buoyancy at all updraft levels (Woods, 1988;Glaze et al., 1997;Graf et al., 1999;Kaminski et al., 2005;Carazzo et al., 2008).
Figure 7 shows the updraft temperature profile for the plume generated in Fig. 6d.The temperature of the plume taken through infrared (IR) imagery for a similar PdF eruption in October 2010 (as no IR imagery is available for January 2010) is shown in Fig. 2. In the buoyant region of the plume, the IR imagery shows a very rapid decrease from several hundred degrees Celcius to mostly ambient temperature within the first tens of metres above the eruptive vent.The temperature decrease in the modelled updraft (Fig. 7) is much slower (temperatures below 200 • C only encountered well above 1 km above the ground).The comparison is very crude and qualitative but at least it shows that the updraft temperature at the base of our simulated plume is not in a correct range of temperatures and consequently the plume is too buoyant and too deep, as not enough fresh ambient air is entrained in the plume base.To correct this discrepancy there is a need to modify the entrainment and detrainment rates at the base of M.EDMF model (as described in Sect.2.2.3).
The question of fresh air entrainment at the base of highly buoyant plumes is actually relevant for all types of hightemperature surface sources inducing convection in the atmosphere, i.e volcanoes but also combustions and in particular biomass fires (Rio et al., 2010).Volcanic or combustion hot gases are extremely buoyant and without entrainment of a large part of fresh air at the base of the buoyant updraft, the latter would accelerate dramatically and, by need of vertical mass conservation, its section would become much thinner than the area of the ground heat source at some altitude above the ground.This is clearly not what is observed in reality, neither in volcanic nor fire plumes.To account for actual observed plumes, the concept of a basal feeding layer in which strong entrainment of fresh air occurs must be introduced.The main questions are how deep is this feeding layer, and how to model entrainment in this layer.Rio et al. (2010) proposed the idea that the entrainment in the feeding layer exactly compensates the narrowing of the plume coverage due to acceleration.They apply this constraint over the full depth of the atmospheric well-mixed boundary layer.
In the present work we keep it as simple as possible.We started from the simple observation that a dominant part of fresh air has been already entrained into the plume within few tens of metres above the ground (Fig. 7).The simplest solution was thus to prescribe a desired fraction α of fresh air at the top of the first model layer (here 40 m above the ground).The relationship between α and the entrainment/detrainment coefficients in the first model layer has been established earlier (Eq.9, Sect.2.2.3).
To compare the Rio et al. (2010) approach with ours, we estimated the fraction α of entrained fresh air at 40 m above the ground, using their assumption (constant updraft section between the ground and 40 m).Transposed to our notations, their Eq.( 15) reads and this yields Assuming no detrainment in the layer, Eq.9 yields With our numerical values (Table 1) the result is α = 68 %, which supports the idea that a dominant fraction of fresh air is entrained into the updraft within a few tens of metres.
The sensitivity of our model to a range of prescribed values of ε and δ in the first model level (i.e.within z above the ground) is here discussed in Fig. 8. Firstly, assuming no detrainment (δ = 0), the plume maximum injection height is found to decrease from about 10 km to almost 0 for ε in the range 0-2 z −1 (Fig. 8a).Beyond 2 z −1 , the volcanic plume does not take off from the ground.A correct value with respect to the LES reference simulation (Fig. 6b) is achieved for ε ≈ 1.8 z −1 .
Figure 8b reveals that the vertical plume development is mostly independent of detrainment.This is not unexpected, since the altitude reached by the plume is in great part driven ).An ensemble of 13 × 9 SCM simulations was considered to build these graphs, with dimensionless values of ε z ranging from 0 to 3, and for δ z from 0 to 2 (by 0.25 steps in both cases).(a) δ = 0 case: altitude where the maximum SO 2 mass detrainment rate is found as function of ε z.The horizontal dashed line is the maximum injection height (1.0 km) found in the reference LES.(b) Contour plot showing the altitude where the maximum SO 2 mass detrainment rate is found as function of both ε z and δ z.(c) Contour plot of SO 2 mass fraction in ambient air in the grid cell at this altitude.In all panels, the black dot marks the location in the ε-δ space of the best-fitted SCM simulation with respect to the reference LES simulation (see text for details).
by its initial buoyancy, the latter being affected only by entrainment but not by detrainment (which does not change the updraft intensive properties such as temperature and water vapour mass fraction).Considering the SO 2 tracer mass fraction at the altitude of maximum detrainment (Fig. 8c), it is found to decrease with increasing detrainment in the first model layer (since less SO 2 mass is left available in the updraft).
Figure 9a shows the SO 2 mass fraction vertical profiles resulting from both the reference LES and the best M.EDMF simulation.The better adjustment in terms of maximum injection height is found for δ = 0 and ε = 1.82 z −1 , corresponding to α = 0.838.The peak SO 2 concentration is found lower (by about 40 %) and vertically more distributed than in the LES.It is clear in Fig. 8c that adding detrainment in the first model level would not improve the result, since detrainment tends to dilute the peak concentration.Despite quantitative imperfection, the M.EDMF model is able to inject a volcanic tracer at the right altitude and at the right order of magnitude in terms of concentration, providing appropriate tuning of the basal entrainment parameter.This simulation is thereafter referred to as the best-fitted M.EDMF simulation.
Up to here, SO 2 mass fractions have served to adjust the ε and δ parameters.Figure 9b also shows the anomaly profiles of the water vapour mixing ratio ([H 2 O]) for the reference LES and the best-fitted M.EDMF simulation.(The anomaly is here defined with respect to the initial water vapour mixing ratio profile, i.e. as [H 2 O](z, t 90 ) − [H 2 O](z, t 0 ), where t 0 and t 90 are simulation times at 0 and 90 min, respectively.)At near ground level, the M.EDMF shows a lower water vapour mixing ratio than the LES model (this is owing to the modification in entrainment at z = z, which imposes a strong increase of the updraft mass flux in the first model level and, in turn, strong divergence of the water vapour mass flux which results in a negative source term for [H 2 O] at the grid scale).At higher altitudes (≥ 0.5 km) the M.EDMF simulation shows comparable agreement and differences with respect to the LES reference as for SO 2 .Figure 9c and d shows the maximum detrainment observed at about 1 km a.g.l. which coincides with the maximum of SO 2 tracer concentration (Fig. 9a) and water vapour anomaly (Fig. 9b).The entrainment and detrainment both reach nearly 0 at around 4 km a.g.l.indicating the maximum height of the updraft and vanishing SO 2 mass fraction at this height (Fig. 9a).
In our simulation, ad hoc fresh air entrainment is prescribed only in the first model layer.The question that arises is whether the fraction α required to achieve the correct injection altitude is dependent of z.To address this issue, a sensitivity experiment was performed whereby a M.EMDF simulation was run with doubled vertical grid spacing ( z = 80 m near the ground) and the same value for ε z = 1.82.This means that α is still equal to 0.838 but now this dilution factor is valid at 80 m above the ground.The resulting SO 2 and H 2 O profiles are shown in Fig. 9a and b, respectively (magenta curves).In terms of peak altitude, intensity and vertical distribution, the profiles are very similar as in the M.EDMF simulation with z = 40 m.This result sug-gests that the required α value is not (or weakly) resolution dependent and (again) that the plume final height is primarily sensitive to its initial buoyancy at the top of the feeding layer, but much less to the depth of the latter.Clearly, further work is needed on the question of the influence of entrainment and detrainment -not only in the first model layer but also at higher levels -on the plume characteristics, but this is the subject of future improvements for our model.

Conclusions
In order to represent deep convective injections of volcanic emissions into the low to mid troposphere in case of effusive eruptions, the EDMF parameterisation by Pergaud et al. (2009) has been adapted.The adapted EDMF scheme takes into account the intense and localised input of sensible and latent heat near eruptive vents and induces a sub-grid convective plume.
We have shown the need to input the specific heat source in order to generate deep plumes using the Méso-NH model by adapting the EDMF scheme.LES simulations were also initialised using water vapour mass flux, sensible heat flux and SO 2 mass flux for the same area and intensities as for the M.EDMF model.In absence of appropriate terrain observations, the LES simulation (considered as a reference) was used to validate the EDMF parameterisation for volcano induced convection (i.e.M.EDMF model).The LES and M.EDMF models have both been successful in generating deep plumes and hence transporting SO 2 tracer to higher altitudes.We have further demonstrated the need to modify the existing lateral mass exchanges a few tens of metres above the localised heat source in the SCM model as without this modification the plumes generated are too deep because of overestimated temperatures a few tens of metres above the ground.The sensitivity of our model to lateral mass exchanges at 40 m above the ground (first model level above the ground) have been presented while further aiding us to tune our model such that SCM results (for SO 2 tracer concentrations) are coherent with the results obtained from LES.
Entrainment of ambient air in a volcanic plume is largely known to be one of the key parameters affecting its buoyancy.Since the first experiments by Morton et al. (1956), extensive research (modelling studies or laboratory experiments) has been deployed to constrain this sensitive parame-ter (e.g.Wright, 1984;Hunt and Kaye, 2001;Kaminski et al., 2005;Carazzo et al., 2008).Although great advances have been made by differentiating between the different regimes (volcanic jets, strong plumes and collapsing columns), it is clear from the comprehensive review found in Tate (2002) and Matulka et al. (2014) that this is still an area of open research.For our case, SO 2 concentrations have served to adjust the parameterisation parameters (prescribed ε and δ within the first model level).The best fit compared to the LES SO 2 profile was obtained with no detrainment and a large fraction of fresh air incorporated into the plume (δ = 0, and ε such that α = 83.8%).The resulting humidity profiles in the LES and SCM show a good agreement as well.
As this parameterisation has been used in an idealised and controlled set-up for one particular case study (January 2010 summit eruption), further work needs to be undertaken whereby the parameterisation is tested for different configurations (i.e.changes in volcanic heat sources; idealised and real case simulations).Furthermore, further investigation is needed on how entrainment and detrainment should be formulated, not only at the base but also at all levels of the updraft.Ideally, a formulation valid at all levels and for a large variety of eruption cases should be sought.

Figure 1 .
Figure 1.January 2010 summit eruption of Piton de la Fournaise: the 60 m long fissure on the inner cliff of Dolomieu summit crater emits lava flows towards the bottom of the caldera.The < 30 m high fountains (left) are the source of the ca. 1 km high vertical plume (right) of gas and vapour.Transport and sedimentation of solid particles are mostly confined to the lowest portion (< 100 m) of the plume.Pictures provided by the Piton de la Fournaise Volcanological Observatory (OVPF/IPGP).

Figure 2 .
Figure 2. Temperature ( • C) of the October 2010 eruption of PdF through infrared imagery provided by OVPF.The temperature scale (right colour bar) ranges between 20.8 and 538 • C. As approximate spatial scale in this image, the crater diameter is about 25 m.(The temperature indication in the upper right corner corresponds to the central pixel marked as a cross.)

Figure 3 .
Figure 3. Figure displaying the input data, the mass flux at ground level (z grd ), and the mass flux at level z grd + z after the incorporation of environmental air mass.The input variables of the model are highlighted in red.

Figure 4 .
Figure 4.The interconnection in terms of the simulation domain between the three sets of simulations performed: Spin-up, SCM and LES.The single cell corresponding to the fissure is tagged for LES.

Figure 5 .
Figure 5. Meteorological profiles from the 3-D spin-up model (1 km resolution) at the location of the January 2010 summit eruption on 2 January 2010 at 10:50 UTC.Vertical profiles of temperature ( • C), potential temperature (K) and water vapour mixing ratio at the grid scale (g kg −1).Altitude displayed is above ground level (a.g.l.), which is at 2600 m a.s.l.

Figure 6 .
Figure 6.(a) and (b): vertical profiles of SO 2 tracer mass fractions (g kg −1 ) horizontally averaged over the 1 km × 1 km domain of LES simulations.(c) and (d): vertical profiles of SO 2 mass fractions in the central grid cell (1 km × 1 km) of the SCM simulations.All simulations were inputted with the same SO 2 surface mass flux.Colour code for all panels: red -30 min, blue -60 min, and green -90 min after model initialisation.

Figure 8 .
Figure 8. Sensitivity of plume characteristics to various entrainment (ε) and detrainment (δ) values prescribed for the first model level (within z = 40 m a.g.l.).An ensemble of 13 × 9 SCM simulations was considered to build these graphs, with dimensionless values of ε z ranging from 0 to 3, and for δ z from 0 to 2 (by 0.25 steps in both cases).(a) δ = 0 case: altitude where the maximum SO 2 mass detrainment rate is found as function of ε z.The horizontal dashed line is the maximum injection height (1.0 km) found in the reference LES.(b) Contour plot showing the altitude where the maximum SO 2 mass detrainment rate is found as function of both ε z and δ z.(c) Contour plot of SO 2 mass fraction in ambient air in the grid cell at this altitude.In all panels, the black dot marks the location in the ε-δ space of the best-fitted SCM simulation with respect to the reference LES simulation (see text for details).

Figure 9 .
Figure 9. Results from the best-fitted SCM simulation with respect to the reference LES simulation (Fig. 6b).The best adjustment was found for δ = 0 in the first model layer and a value of entrainment ε such that the fraction of fresh ambient air entrained at the top of the first model layer is α = 0.838.(a) Compared SO 2 profiles in the SCM and LES.(b) Compared profiles of water vapour anomaly (defined as departure from the initial profile).(c) Profiles of dry air mass entrainment and detrainment rates (ε M u and δ M u , respectively).Where detrainment dominates, the difference is shaded in magenta.(d) Profiles of water vapour and SO 2 mass detrainment (δ [H 2 O] M u and δ [SO 2 ] M u , respectively).

Table 1 .
Variables and values used for LES and SCM models.