Numerical model of gas dispersion emitted from volcanic sources

An Eulerian model for passive gas dispersion based on the K-theory for turbulent diffusion, coupled with a mass consistent wind model is presented. The procedure can be used to forecast gas concentration over large and complex terrains. The input to the model includes the topography, wind measurements from meteorological stations, atmospheric stability information and gas flow rate from the ground sources. Here, this model is applied to study the distribution of the CO2 discharged from the hot sources of the Solfatara Volcano, Naples, Italy, where the input data were measured during a 15 day campaign in June 2001 carried out to test an Eddy Covariance (EC) station by Osservatorio Vesuviano-INGV, Naples. Mailing address: Dr. Antonio Costa, Dipartimento di Scienze della Terra e Geologico-Ambientali, Università degli Studi di Bologna, Italy; e-mail: costa@ov.ingv.it () Present affiliation: Centre for Environmental and Geophysical Flows, Department of Earth Sciences, University of Bristol, Wills Memorial Building, Queen’s Road, Bristol BS8 1RJ, U.K.


Introduction
Many volcanic and non-volcanic areas in Italy emit a huge amount of gas into the atmosphere.One of the most frequent gases discharged from both volcanic (e.g., Solfatara volcano) and non-volcanic sources (e.g., Central Italy vents), is carbon dioxide (CO 2 ) which has a molecular weight greater then air.
For this reason, under stable atmospheric conditions and/or in the presence of topographic depressions, its concentration can reach high values resulting in lethal effects to humans or animals.In fact, several episodes of this phenomenon were recorded in Central Italy (Rogie et al., 2000) and in different areas.One of the most tragic examples was the degassing of Lake Nyos, Cameroon in 1986 when a dense cloud of carbon dioxide hugging the ground suffocated more than 1700 people in one night (Clarke, 2001).
The cloud dispersion of gas denser than air released from a given source is governed by the gravity and by the effects of lateral eddies which increase the mixing with air around the edges of the plume, decreasing the density.In the initial phase the negative buoyancy controls the gas dispersion and the cloud follows the ground (gravitational phase).In this phase, the dispersion of heavy gas is markedly different from a passive or a positively buoyant gas dispersion.When the density contrast is not important, gas dispersion is basically governed by the wind and atmospheric turbulence (passive dispersion).
Although from a theoretical point of view, gas dispersion can be fully studied by solving the complete equations system for mass, momentum and energy transport, in actual practice, different simplified models able to describe only specific phases or aspects are commonly used.These models range from the simplest analytical Gaussian models to the more complex Computational Fluid Dynamics (CFD) models.
Gaussian plume models assume that the crosswind concentration has a Gaussian profile with the standard deviation of the distribution function of the downwind distance from the source and of the atmospheric conditions.A Gaussian plume profile is also assumed in the Puff model but the release is divided into a number of different puffs and each puff is independently modeled.The final concentration at any point is found by a superposition of all puffs.Although these models are frequently used, several studies have shown that their validity is limited (see Willis and Deardoff, 1976;Pasquill and Smith, 1983;Briggs, 1985) and, for example, other analytical non-Gaussian solutions are in better agreement with the observed data (see e.g., Hinrichsen,1986;Lin and Hildemann, 1996).
Another common approach is given by the Box (or Similarity) models which describe the integral properties of plume.A set of differential equations for averaged mass, momentum and energy balance is solved along the plume using different simplifying similarity assumptions (e.g., Blackmore et al., 1982).SLAB (Ermak, 1990), HEGADAS (Witlox, 1994) and DE-GADIS (Spicer and Havens, 1989) are popular examples of these similarity models.
The most complete and computationally most expensive models are given by the threedimensional CFD models based on the transport theory of mass, momentum, energy and species.This approach is able to simulate dispersion of an heavy gas accounting for obstacles and topographic effects, variation of atmospheric conditions and wind direction, etc.
A compromise between the complexity of CFD models and the simpler integral models is given by the shallow layer approach which uses depth-averaged variables to describe the flow behavior (Hankin and Britter, 1999;Venetsanos et al., 2003).These models are used to describe gravity driven flows of dense gas over complex topography.
Finally, another approach to describe gas dispersion in different way is given by the Lagrangian particle models that simulate the random walk of many distinct mathematical particles in the mean flow (see Pareschi et al., 2001) for an application in a volcanic environment).
For a less dense gas which is passively driven by wind advection and atmospheric turbulence, simpler models based on the advectiondiffusion equation can be used, like the model we are presenting in this paper.The model is based on an efficient fully explicit solver of the advection-diffusion equation (Sankaranarayanan et al., 1998) coupled with a mass consistent wind model (Douglas and Kessler, 1990).
In particular the model is proposed as an operational procedure for forecasting gas concentration over a complex terrain such as the area surrounding the Solfatara crater in the Campi Flegrei.In fact, Solfatara volcano, releases 1500 td −1 of hydrothermal CO 2 through soil diffuse degassing from a relatively small area (0.5 km 2 ) (Chiodini et al., 2001;Cardellini et al., 2003).
The study of the dispersion of «cold» CO 2 discharged from many vents in Central Italy (Rogie et al., 2000) needs a model able to capture the different behavior of dense gravity driven gas flows such as those based on the shallow layer approach.Instead, in the Solfatara case, temperature of the emitted gases is relatively high (fumaroles temperature is between 96°C and 162°C (Chiodini et al., 2001) and, flux weighted temperature of diffusing soil is 66°C).This implies a density decrease that almost balances the increase due to the greater molecular weight (M CO2 /M air = 44/29 while T CO2 /T air ∼ 400/300).Moreover low wind conditions are rarely recorded at Solfatara (during the last six years daily average wind intensity w>1 m/s were 75% of the total and the cases with w<0.1 m/s 0.05% of the total) and in any case the passive plume assumption is legitimate since the condition ∆ρ=(1−M CO2/Mair)C=15/44C<0.001 kgm −3 (Mohan et al., 1995) is already observed at few tens of centimeters level (C is the CO2 concentration).These reasons show that the simple model proposed here, based on the advection-diffusion equation, is able to describe the relevant physical features of the gas dispersion in the case of hot sources of gases, as in Solfatara example.

The model
A full description of gas dispersion into the atmosphere requires the solution of mass, mo-mentum, energy and species transport equations of Macedonio and Costa (2002).This implies a relatively high computational effort and sometimes it is impossibile to study the gas dispersion in complex and large domains, especially in the surface layer where a large range of length scales are involved.
In order to reduce the computational time we introduce some assumptions.Our main objective is to solve the advection-diffusion equation for the gas concentration c=〈c〉+c′ where V=(u x +u x ′, u x +u y ′, u z +u z ′) is the wind field and 〈Q〉 is a source term.Terms inside the brackets 〈〉 represent the average part, terms with the symbol prime «′» their turbulent fluctuations.
As, for example, in Prabha and Mursch-Radlgruber (1999a,b), we do not solve the complete set of the coupled equations for mass, momentum, energy and concentration, but we use the wind field as given by a diagnostic wind model which produces a zero three-dimensional divergence velocity field consistent with the measured values, avoiding artificial generation or loss of gas.Finally, turbulent terms are parameterized according to the K-theory.

The Diagnostic Wind Model (DWM)
To evaluate the wind field, we use the Diagnostic Wind Model (DWM) developed by the United States Environmental Protection Agency (EPA).This model generates the wind components (U, V, W) in a terrain following the coordinates system.
The DWM needs topographic data, average wind on the computational domain and atmospheric stability information within the scale of the domain (the temperature gradient dT/dz).
In a first step the domain-mean wind is adjusted for the kinematic effects of terrain (lift- ing and acceleration of the airflow over terrain obstacles), thermodynamically generated slope flows, and blocking effects.In a second step, wind observations, when available, are added to the first step field, and an objective analysis scheme is used to produce a new gridded field (U, V, W).The scheme is designed so that the observations are used to define the wind field within a user-specified radius of influence while the first step (U, V, W) field is used in subregions in which observations are unavailable.
Finally, a divergence-minimization procedure in terrain following coordinates is iteratively applied until the inequality (2.2) is satisfied (U=(U, V, W) and ε is an arbitrarily user defined small number).The final product of DWM is an approximately null-divergence wind field consistent with the observations (for further information see Douglas and Kessler, 1990).
The approximation of null-divergence wind field is generally applicable up to half a kilometer above ground level (Dutton and Fichtl, 1969) therefore to our study, since we treat the surface layer, i.e. the lowest part of the Planetary Boundary Layer (PBL).

Turbulence parameterization
Turbulent fluxes are given by the product of the fluctuation terms, and in agreement with the «K-theory», they are expressed in terms of gradients of average concentration where K=(Kx, Ky, Kz) is the turbulent diagonal diffusivity tensor.
Inside the atmospheric surface layer, the Monin-Obukhov similarity theory allows us to where l is the von Karman constant (l= 0.4), z is the distance from the ground and φ h is the atmospheric stability function for temperature (2.5) with β h =7.8, γ h =11.6 and Pr t ≈0.95.Therefore, evaluating K z requires knowledge of the friction velocity u * and the Monin-Obukhov length L. From a practical point of view, L is often a quantity difficult to estimate directly.Instead, it is operationally easier to evaluate this length using the «Bulk Richardson number» method.The Richardson number Ri is a measure of the dynamics on the buoyancy effects: instability when Ri<Ri crit and stability for Ri>Ri crit (Ri crit 0.2).The analytical relationship by Golder (1972) allows us to calculate the «Bulk Richardson number» Ri b as (2.6) where g is the gravity, θ and T the potential and the absolute air temperature (Jacobson, 1999), z m the geometric mean thickness of the considered layer and u m is the logarithmically interpolated wind speed at zm.
The knowledge of Rib permits us to analytically express z/L in terms of Rib and of the roughness length z0 in agreement with Byun (1990) (for an improvement of this method see also De Bruin et al., 2000).
The friction velocity is given by the similarity theory (2.7) where z1 is the height at which the wind veloc- ) is known and Ψ m (x) the atmospheric stability function for the momentum (see e.g., Dyer, 1974;Jacobson, 1999).Following the large eddy approach, the horizontal eddy diffusivity K x =K y =K h is simply assumed as discussed in Pielke (1974) (see also Pielke, 1984;Physick, 1988;Boybeyi and Raman, 1995;Azad and Kitada, 1998) (2.8) were α is a dimensionless constant empirically determined of order one and ∆x and ∆y the grid spacing (in this study we set α=3).A lower limit of K H =1 m 2 /s was imposed.

Model solutions and numerical aspects
Under the previously introduced assumptions and considering a null-divergence wind field, in a terrain following coordinate system (Douglas and Kessler, 1990) given by (2.9) (h(x, y) is the height of the topography), the eq.(2.1) for the mean scaled concentration C may be re-written in a generalized form as (Toon et al., 1988;Jacobson et al., 1996;Park and Kim, 1999) (2.10) where (U, V, W) are scaled wind speeds, Kh and Kz are scaled diffusion coefficients and Q * the source term in the new coordinates.In the case of the simple trasformations (2.9) the Jacobian is equal to the unit, J= ∂z/∂z * =1 , therefore c(x,y,z)=C(x,y,z * ) and Q=Q * (for a more general treatment see e.g., Toon et al., 1988;and Byun, 1999).In this work, a third-order upwind numerical scheme was used as given in Kowalik and Murty (1993) and tested for advection-diffusion equation in Sankaranarayanan et al. (1998), but we improved the original scheme by introducing different limiter methods.For example, in the x-direction we discretized as (2.11) for Ui,j,k≥ 0 and i=2,…, N x −1; instead for U i,j,k ≤0 and i=2,…, N x −1we used (2.12) Concerning the boundary points, we evaluated the convection term using a first order upwind.Near the boundaries, i.e. i =1 and i=N x (i ranges from i=0 to i=N x +1), we set (2.13) for U i,j,k ≥ 0, whereas (2.14) for U i,j,k ≤ 0.
The diffusion terms are evaluated using a central difference scheme for the general case with a not uniform turbulent diffusivity tensor, e.g., in the z-direction using (2.15)As shown by Sankaranarayanan et al. (1998), this scheme is numerically stable, when the following inequality is satisfied (2.16) Concerning the z-direction, we generalized the numerical schemes (2.11), (2.12) and (2.15) to the non-uniform vertical grid case.Finally, in order to improve the third-order numerical scheme, different limiter methods were introduced (e.g., Superbee, Van Leer, MinMod, etc.).Such limiters strongly reduced the numerical over-and under-shooting that commonly affect these schemes near discontinuities.In fact these methods preserve the monotonicity of the solution while the accuracy remain higher than the first order upwind methods (Sweeby, 1984).
In particular the «superbee» limiter was used in the applications described below.

Applications
Although this model can be applied to simulate gas dispersion on a large domain, this paper presented an example of application on a relatively small area (2000×1600×100 m).In fact, this study focuses on gas dispersion from the fumarolic vents inside the Solfatara crater which releases a large quantity of H2O, CO 2 and a small fraction of H 2 S into the surrounding area (Chiodini et al., 2001).The total amount of CO 2 discharged from the entire anomalous area is ∼1500 td −1 (Chiodini et al., 2001).

Input data
Because of the lack of meteorological stations around the crater at the present time (in the future there will be a fixed EC station), we chose a relatively small computational domain of 2 km× ×1.6 km ×100 m around the crater (see figs. 1 and 2) and used meteorological input data performed in the period between 8th and 25th June 2001 and published in Werner et al. (2003).Finally, the simulated concentrations were compared with the few available measurements of Werner et al. (2003).For all the simulated cases, flow rates at the ground were assumed to be the same as the real measured values over the entire degassing area (see fig. 1) as obtained during periodic monitoring by the Osservatorio Vesuviano (see e.g., Chiodini et al., 2001).
From the measured EC episodes (Werner et al., 2003) we chose data corresponding to two opposite cases of atmospheric stability.The first (8th June 2001) corresponds to a positive value of the Monin-Obukhov length L=2.2 m (stable case), the second (13th June The computational domain was discretized by a (40 ×32 ×11) grid a horizontal resolution of 50 m and a variable vertical grid spacing finer near the surface.
At the top and lateral boundaries, free flow conditions were assumed, whereas at the bottom boundary no diffusion flux was assumed (zero gradients).Moreover, we imposed different boundary conditions for outgoing or ingoing flux: for outgoing flux, zero derivative conditions were set, whereas for ingoing flux, null concentrations at oundaries were set.The choice of the lateral boundary conditions is important to avoid absorption or reflection from these boundaries.
These meteorological measurements include the Monin-Obukhov length L, friction velocity u * and wind speed and direction at 2.5 m 2001) to a small negative value L=−2.7 m (unstable case).As previously explained L is a measure of the relative importance of mechanical mixing to buoyancy turbulence.When L is positive, then the atmosphere is stable.When L is small and negative, mixing induced by buoyancy is more important than mechanical mixing and the atmosphere is highly unstable.Finally when L is negative and large, mechanical mixing overcomes mixing due to buoyancy and the atmosphere is weakly unstable.

Results
As previously said, we show the results of the simulations in two different cases of atmospheric stability.For each case we report a graph of the CO2 concentration at 2.5 m level and two different views of the integral concetration in the terrain following coordinates system: a planar and a lateral view respectively.Obviously these views can give an idea of the gas distribution but they should not be confused with the Table I.The column labeled with Cs represents the value at 2.5 m level at location 1 of (Werner et al., 2003) of the simulated gas concentrations.The measured values at 2.5 m level in the location 1 of Werner et al. (2003) are reported in the column labelled with Cm.A background CO2 concentration of 0.70× ×10 −3 kgm −3 was considered.corresponding views in a normal Cartesian system.
Table I compares both the measured and the simulated CO 2 concentrations at 2.5 m level.Considering the limitations due to the Eulerian approach to simulate the proximal range (Boybeyi and Raman, 1995), the obtained concentrations are in good agreement with the observed values (in both cases the obtained values above the background CO 2 concentration are within three times the measured values).
The obtained results reported in the figs.3, 4 and 5 and table I show a satisfactory performance of the presented model.Moreover, from figs. 3 and 4 it is evident that in the unstable case, gas is rapidly transported into the upper part of the atmosphere beyond the computational domain, while during the stable case most of the gas is transported in proximity of the surface.
As previously said, despite the relatively short distance of point 1 from the source and the simplicity of the model, in both cases the simulated concentrations are quite close to the measured values.Moreover the performance of the model should improve with increasing distance from the source.coupled with a K-type turbulent diffusivity and a mass consistent wind model.The model can be generally applied to estimate the gas concentration over complex and large domains in different atmospheric conditions.As input data it requires the heights of the topography, wind measurement by meteorological station, the atmospheric stability information and gas source fluxes data.Here, as an example, we tested the model simulating the distribution of the CO2 discharged from the Solfatara Volcano, Naples, Italy, comparing the results with the measured concentrations during the campaign of June 2001.The model simulation shows, in the Solfatara case, the need to have a reliable and realtime available measurement from an in situ meteorological station.Generally, the obtained results show a satisfactory agreement with the actual values.This indicates the usefulness of this simple model and its potential as a tool for interpretation of plumes data and for gas hazard assessment.
turbulent diffusivity K z in terms of friction velocity u * and the Monin-Obukhov length L (2.4)

Fig. 1 .
Fig. 1.Planar view of the investigated area around the Solfatara crater, Naples, Italy.Coloured part represents diffuse degassing area; the red rectangle encloses the periodically monitored area.

Fig. 2 .
Fig. 2. Shaded relief of the zone surrounding the Solfatara crater.The maximum height in the domain is about 203 m a.s.l.The resolution of the used topography model in our discretization was 50 m×50 m.

Fig. 3 ..Fig. 4 .
Fig. 3. Planar and lateral views (top and bottom respectively) of the integral concentration in terrain following coordinates for the 8th June 2001 event (stable).Time is in seconds, lengths are in meters and concentration values in kg/m 3 .Fig. 4. Planar and lateral views (top and bottom respectively) of the integral concentration in terrain following coordinates for the 13th June 2001 event (unstable).Time is in seconds, lengths are in meters and concentration values in kg/m 3 .

Fig. 5 .
Fig. 5. Simulated gas concentration values at 2.5 m level for the 8th June 2001 (top) and the 13th June 2001 (bottom) events.Time is in seconds, lengths are in meters and concentration values in kg/m 3 .