Atmospheric dispersion modelling of CO 2 emission in the Colli Albani volcanic district ( central Italy )

Carbon dioxide is a gas denser than air, and its point-source ground emission from natural systems or from areas impacted by CO2 injection underground may result in hazardous accumulation, especially in topographically-depressed sites. The use of atmospheric dispersion numerical models helps predicting the dispersion of the CO2-enriched gas plume once emitted from underground and allows an accurate map of hazard level through time under particular meteorological conditions. In this study, the accuracy of atmospheric dispersion simulations has been tested using a natural system of CO2 emission to atmosphere from underground in an area called Solforata di Pomezia, near the city of Rome in central Italy. This area is located in the Alban Hills, which underwent volcanic activity during the Quaternary, and is characterised by low permeability volcanic and sedimentary formations that allow the accumulation of gas at shallow depths and below surface. This site has been long investigated in terms of soil CO2 emission rates, which range from 44 to 95 ton∙day-1. Using the TWODEE2 numerical code, a number of simulations were performed considering a set of combined CO2 soil flux emission and meteorological (wind, temperature) from literature. The results fit well in the range of measured CO2 concentration in air at distinct heights in the site. The model does not predict lethal gas concentration at heights 1 and 2 m above the ground based on actual soil emission rate (95 ton∙day-1). Two probabilistic models were developed with emission rate five (500 ton∙day-1) and ten (1000 ton∙day-1 times bigger than nowadays but still no hazardous levels were predicted.


Introduction
The study of CO 2 emissions in natural systems provides a valuable information on the assessment and quantification of potential hazards related to underground carbon storage leakage.Emissions of CO 2 are studied in a large variety of geological environments, i.e.: sedimentary basins, active and non-active volcanic areas, seismically-active regions, and geothermal fields.Because of the physics of carbon dioxide gas, e.g., colourless, odourless, higher density compared to air, its accumulation may result hazardous and even lethal for life.The impact of CO 2 concentration in air on human beings has recently been reported by De Lary et al. [2012], who summarized the effects of its inhalation as a function of the concentration and duration of exposure.Being denser than air, CO 2 tends to accumulate on the ground filling depressions and forming the so-called CO 2 lakes, which could represent a serious hazard for human beings.This accumulation is governed by gravity, turbulence and dispersion.
This research article deals with the CO 2 emissions from a geological site located at the southern border of Rome municipality, central Italy.The site is called Solforata di Pomezia, taking its name by the typical smell of rotting eggs and bright yellow crystalline solids, and is known for a widespread degassing of mantle-derived CO 2 [Chiodini et al. 2004] (Figure 1).The objective of this research is to quantify the risk of CO 2 accumulation in leakage episodes in geological storage by using atmospheric model simulations.The robustness of these numerical calculations is based on the calibration with measurements of air composition available in the literature.The results of this numerical model may provide clues for the determination of environmental conditions under which the site and other CO 2 -emitting areas can host hazardous conditions for human life.
Atmospheric dispersion modeling is a tool that provides information to know how volatile compounds disperse in the atmosphere from a source located underground or on ground.It is used to predict future air concentration under specific risk scenarios, by developing maps of possible hazardous levels for human beings, and it can be useful in the risk assessment of geological storage of CO 2 .

Geological setting
The studied site is located in the Alban Hills region, in the Tyrrhenian sector of Central Italy (Figure 2).This area is characterised by low permeability volcanic and sedimentary formations that allow the accumulation of gas at shallow depths and surface.Also, this region underwent volcanic activity during the Quaternary, related to an extensional tectonic regime whose evolution is directly controlled by primary regional structures connected to the recent geodynamic evolution of the Tyrrhenian basin-Apennine Chain system [De Rita et al. 1995].The Alban Hills volcano is one of the main volcano craters of the peri-Tyrrenian belt, starting erupting 0.7 Ma ago [De Rita et al. 1988].Volcanic activity ended in the Holocene 2.7 ka with a phreatomagmatic phase that formed a large number of craters in the south-western part of the volcano [Selvaggi and D'Ajello Caracciolo 1998].The origin of CO 2 underground emission in the Alban Hills is related to buried highs of the carbonate basement, which act as traps for the gas of deeper origin migrating through a NE-SW fault extending from the coast to Albano Lake (Chiodini and Frondini 2001].The observable phenomenon of continuous strong degassing that is widespread in Latium district, suggests that a large amount of gas is being accumulated underground during long periods of time with short episodes of strong gas release [Voltaggio andBarbieri 1995, Funicello et al. 1992]. The Solforata di Pomezia is a natural degassing area located a dozen of km south-west of the city of Rome, being part of the quiescent Alban Hills Volcano district.This area is well known for a diffuse natural gas  leakage and the presence of sulphide-rich waters; also, it is affected by a significant uplift rate (1.3 mm•y -1 ) suggesting the presence at depth of a pressurized gas source [Tolomei et al. 2003].
The morphology of the Solforata di Pomezia has been modified in the last century, especially when was mined for sulphur.After mining, the site was abandoned and today it is periodically flooded during rainfall episodes.Presently, the area consists of a depression of 270,000 m 2 and is close to an industrial area.In the centre of the depression, there is a lake that typically dries out in summer.The site has a quite rectangular shape with an east-west orientation with a vertical slope between 20 and 30 metres.A drainage channel runs towards west at the southern border of the depression.Vegetation is present all around except for the area around the lake where mud is present.In Figure 3 the slope around the emission area, the location of the lake and the absence of vegetation around the lake are clearly observed.
During the summer period, due to the drying out of the soil, main gas emission points are easily recognisable (Figure 4).

CO 2 soil flux emissions
A first estimation of the CO 2 soil flux in the Solforata di Pomezia was made in 1996 by Chiodini and Frondini [2001] using an accumulation chamber, reporting 110 measurement stations over an area of 55,000 m 2 .In 2003 Carapezza et al. [2005] made another survey in the Solforata area covering 30,000 m 2 with 200 accumulation chamber measurements.Another soil flux survey was carried out in 2007 [Carapezza et al. 2012] extending the soil gas mapping to an area of 229,000 m 2 (356 measurements and 7 Tuneable Diode Laser, TDL, profiles).The latter set of data has been selected for this study; in particular, 278 of 356 accumulation chamber measurements were used (black dots in Figure 5), 7 TDL profiles (black segments in Figure 5), and the wind speed and direction data (Figure 8).The selection of only a part of the accumulation chamber survey is due to cross-check data with the TDL profiles, the wind speed and direction data, and the atmospheric dispersion outputs over the same emission area.The diffuse soil flux of CO 2 has been measured by two portable accumulation chambers [Chiodini et al. 1998, Carapezza andGranieri 2004], equipped with a Li-820 infrared CO 2 detector (0-2 vol.%).The final output of Carapezza et al. [2012] is the soil flux map shown in Figure 5.
The flux map clearly shows the main emission area around the pool, and almost 6 emission areas can be spotted.The main soil emission area is longitudinal to the valley axis.

Air CO 2 concentration
Air concentration has been measured by Carapezza et al [2012] using a Tuneable Diode Laser (TDL, see Figure 5, and Table 1) over the lake (lines 1, 3, 5 and 6) and the channel (lines 2, 4 and 7).This methodology allows the measurement of the average air concentration of CO 2 on a profile length every few seconds [Schiff et al. 1994, Tittel et al. 2006, Weber et al. 2005]; simultaneously, the wind direction and speed are recorded by a sonic anemometer [Carapezza et al. 2012].TDL profiles were recorded at two heights, 20 and 25 cm from surface; two of them were over the lake and one was along the channel.According to De Lary et al. [2012], air CO 2 concentration is admissible under 5% for human beings.This is considered a threshold beyond that breathing becomes extremely laboured with occurrence of headaches, sweating and bounding pulse.1) will be used for a comparison with soil flux data obtained in the same period.
TDL measures CO 2 concentrations as average values on the total length of the infrared laser; consequently, some areas along the laser path likely have higher concentrations than the average.Results shown in figure 6 are concentrations measured on the 68 m profile over the lake lasting a couple of hours, CO 2 reached a maximum concentration of 3026 ppm and a minimum of 139 ppm.Considering the Solforata di Pomezia in a steady state, concentrations recorded over the lake by TDL (Figure 6) are indicative of the level of CO 2 air concentration expected over the site.The data from a 20 hours TDL measurements along the main channel for a length of 118 metres are shown in Fig- ure 7.For this profile, maximum CO 2 concentration reached 3500 ppm in the early morning hours while in the night maximum concentration was around 600 ppm and 360 ppm before 8 p.m..These data will be compared with numerical outputs of model S1 in next section.As can be observed in both plots air concentration used to vary as wind speed changes.Values for the TDL in Figure 6 are higher on average than in Figure 7 because the first profile is over the main emission area (the lake) while the second profile is along the main valley axis, so a decrease of CO 2 in air due to dispersion is expected.

Metereological data
Carbon dioxide is a gas denser than air (1.976 kg•m -3 ) so that accumulates in topographic depressions like the bottom of a valley or in ditches, and under stable atmospheric conditions CO 2 concentrations can reach high values.Dispersion of gases denser than  air is governed by gravity and by the effects of lateral eddies which decrease the plume density through the incorporation of surrounding air.In the initial phase the buoyancy controls the gas dispersion and the cloud follows the ground (gravitational phase); in contrast, when the density contrast becomes less important, gas dispersion is mainly governed by wind and atmospheric turbulence (passive dispersion phase) [Costa et al. 2005, Folch et al. 2007].
Wind speed was recorded in the Solforata channel for 20 hours during May 14 th and 15 th 2007 using a sonic anemometer every 15 seconds.The recorded wind data at 1 m height are shown in Figure 8. Three wind direction families are recognised: the main wind direction blew from ENE to WSW reaching a maximum speed of 3.3 m•s -1 , and the other two blowing directions are from SSE and from W, is the latter being the weakest family for wind speed on average.Other wind directions were recorded but with lower wind speeds.
Unfortunately, atmospheric pressure could not be measured simultaneously with a similar frequency.This is a significant uncertainty for the input data implementation in the numerical model.As indicated in Gasparini et al. [2015], atmospheric pressure variations over a leakage gas site can move up or delay the emission from soil to atmosphere.Nevertheless, no abrupt change in pressure condition happened during measurements.
The first law of Gay-Lussac (also known as Charles's Law) says that at a constant pressure the volume of a gas is directly proportional to the absolute temperature (V/T=k).This law describes how a gas expands as the temperature increases.Soil respiration is responsive to temperature and increases when soil temperature are raised experimentally [Billing et al 1982, Van Cleve et al. 1990].Temperature is the single best predictor of the annual soil respiration rate of a specific location, but inclusion of precipitation in the regression does increase the predictive power of the model [Raich and Schlesinger 1992].
Air temperature was also recorded by the sonic anemometer during the TDL experiment, and ranged from 12.5 to 27ºC (Figure 9).The highest CO 2 concentration is coupled with periods with lower temperature.
It is clearly visible in the chart below how emission is concentrated during the night time.Generally looking, temperature starts from 25°C and decreases to 20°C at dusk.During the night the temperature still kept a descending trend till the early hours of the morning.In particular, have to observe that during the night there are many temperature oscillations with changes of 5° in half an hour.Comparing plots of CO 2 vs Temperature (Figure 9) and CO 2 vs Wind (Figure 7), temperature does not decrease when wind blows up, on the contrary CO 2 concentrations are greater when wind speed decreases.
In order to ascertain the relations amongst air CO 2 concentration and atmospheric conditions, we performed a forward-stepwise multiple regression of air (CO 2 ) and environmental data: this is a largely applied statistical method to identify the contribution of several variables on the total variance of a dependent variable [see Laiolo et al. 2016 and references therein).
Given the horizontal wind speed, the vertical wind speed, the wind direction and the air temperature measured by the anemometer, we could estimate the contribution of each on air (CO 2 ) and model CO 2 fluctuations due to variations in air temperature and wind conditions.
Results of the forward-stepwise multiple regression indicate that atmospheric variations are able to predict the 36.3% of the total CO 2 concentration variance, and that only air T and horizontal wind speed are statistically significant, and anti-correlated, on CO 2 fluctuations.The output model in Figure 10 shows how air T and wind speed variations can model part of the measured CO 2 .

Numerical code
The numerical simulations performed in this work were done using the TWODEE2 code [Costa et al. 2008].This code is based on shallow layer time-de-pendant Eulerian approach for dispersion of heavy gases, and it is an updated version of the TWODEE code [Hankin and Britter 1999a, b, c].TWODEE2 is based on depth-averaged equations obtained by interpreting mass, density and momentum equations over the fluid depth, from the bottom up to free surface.Such approach is able to describe the time evolution of gas concentration, depth-averaged velocity, and averaged cloud thickness.The input data include the topography, ground roughness, wind measurements and gas flow rate from the ground sources [Folch et al. 2007].In this study all simulations were developed using the meteorological processor DIAGNO [Folch et al. 2008] coupled with TWODEE2.DIAGNO uses wind data at a point of the domain and, assimilating terrain information, generates a zero-divergence wind field.

Numerical model set-up
The numerical model is based on a 2D mesh of 17500 regular quadratic elements of 20x20 m, and the topography of the site has been implemented from cartographic data with a scale of 1:25000 (see Figure 3).
No flow conditions at the boundary were imposed and an initial CO 2 concentration in the air of 400 ppm, the atmospheric background concentration, was considered.
All numerical simulations used wind data shown in Figure 8, and as wind has been recorded for 19 hours, the simulation time was 19 hours.
A first simulation, called S1, was based on mete- orological data and estimated emission rate (95 ton•d -1 ) from Carapezza et al. [2012] field data.The model S2 was built using data recorded by Chiodini and Frondini [2001], with an emission rate of 44 ton•d -1 .For this simulation, the atmospheric data was that recorded by Carapezza et al. [2012].Models S3 and S4 had respectively an emission rate of five and ten times bigger than model S1.
Table 2 reports all the model features of the performed simulations.

Results
A first significant prediction from the simulations is that the gas plume covers the emission area rapidly after short periods of time.In model S1, the gas would need 2.5 hours to cover the channel, 4.5 hours for model S2, an hour for S3 and half an hour for model S4, respectively (Figure 11).
The predicted plume in model S1 expands towards the west side of the area steering into the left channel (Figure 12).A small part of the plume moves to the east side remaining close to the emission area.Model S1 predicts a maximum air concentration of 9000 ppm in two small areas while the average air concentration is under 2500 ppm.
Similar behaviour to S1 is predicted for model S2 but with lower gas concentrations.The plume expands mainly over the emission area and it finds its preferential way to the west side, and only a small plume goes to the east part.The plume shape is narrower than model S1.Model S2 predicts a maximum concentration of 6000 ppm but air concentration on average does not raise over 1500 ppm (Figure 13).
Simulations S3 and S4 predict maximum air concentration of 15000 and 54000 ppm, respectively (Figure 14).For those models, the predicted plumes expand along the valley's axis to the north-west part of the channel covering an area much wider than the one predicted by S1 and S2.Interestingly, even the northeast area is predicted to be affected by CO 2 -enriched air.Two spots are visible for S3 with concentration under 5000 ppm, while S4 spotty concentrations are under 13000 ppm.
Model S4 covers an area that is quite the double of S1 at 25 cm height, without reaching a dangerous level (air CO 2 1.5%).Higher concentrations (5%) are reached in spot points over the emission area which is normally flooded and difficult to reach or to pass through, so that the Solforata di Pomezia is not a dan-   gerous area even with an emission rate of ten times bigger than actual rate but always with wind activity.In all four simulations, the CO 2 plume follows the longitudinal axis of the studied area with higher concentration predicted over depressions existing in the site.Plume horizontal dispersion follows the topography covering the lower area on the west side.Air concentration varied as wind intensity changes.On the other hand, areas with higher and lower CO 2 air concentration change location as wind changes blowing direction.
Comparing data from models S1, S3 and S4 at 1 m and 2m high from ground surface, a strong increase in the air CO 2 concentration is observed in the S3 and S4 cases (Figure 15, Figure 16 and Figure 17).The predicted plume of model S1 at 1 metre has an average concentration of 470 ppm with a maximum concentration of 3895 ppm.At 2 metres,very few spots are still predicted with a high concentration (maximum up to 3353 ppm) while the minimum level decreases down to 450 ppm on average.Models S3 and S4 show a larger plumes that cover areas never reached in model S1.Concentrations are higher especially above the emission area, and in the inner part of the plume for a length of 1 km till the bottleneck present on the middle west side; In the west side, concentrations decrease generally down to 2000/3000 ppm.Model S3 has an average value under 8000 ppm at 1 metre over the lake and 6000 ppm on the far west side (Figure 15).
Model S4 shows a similar behaviour as S3 with higher concentrations over the lake and a significant decrease on the west side.Interestingly, topography affects the plume dispersion in all simulations.The bottleneck area prevents plume dispersion sideways and allows the plume to reach higher elevations up to 90 metres a.s.l.Concentration in model S4 reaches 16000 ppm over the lake and decreases down to a maximum of 7000 ppm on the west side at 1 metre height.At 2 metres height concentration over the lake has an average value of 9000 ppm and 4500 ppm on the west side.

Discussion
The model predictions yields a significant difference on the plume size between simulations S1 and S2 outputs (see Figure 12 and 13).Also, a higher average air concentration is predicted in S1 due to the greater amount of emitted CO 2 .If the model data from S1 is compared to TDL data from measurements by Carapezza et al. [2012], it is concluded that the experimental data are always between the predicted minimum and maximum concentration and that the numerical model well reproduces the concentration changes caused by wind activity.
The mapping of hazard assessment in geological active regions and areas affected by underground storage of CO 2 is useful to forecast which zones would be more impacted in case of episodic emission.In this respect, simulations S3 and S4 were performed to predict possible hazardous scenarios with higher emission rate than nowadays, five and ten times greater than model S1, respectively.As reported above, the results obtained using present-day conditions (model S1) are comparable to real concentration data in the site.On that basis, the hazard map considering a more intense emission can be produced based on atmospheric dispersion models at different heights.
Model predictions show that under some emission conditions air concentration can reach values higher than 1%.However, even in the worst scenario the concentration is below 10 % at heights of 1m or higher so that major impact on human beings (emission rate of 1000 t•d -1 ) vertigo, headache, vomiting, shortness of breath, loss of mental ability, muscular weakness, drowsiness would not be expected.
Due to the lack of atmospheric pressure data each simulation has considered a constant atmospheric pressure of 990 hPa disabling the suction-pump effect.The barometric pumping effect relates to the diurnal changes of ambient pressure due to change in the temperature, as well as pressure changes due to weather systems such as mid-latitude cyclones.This variation in pressure at the surface acts as a suction-pump on the underground [Carrigan 2010].The wind effect drives to an enhanced mixing layer due to its turbulent character, which increases the rate of seepage into the atmosphere [Oldenburg et al. 2010].Low-frequency atmospheric motions are less effective than higher-frequency motions at trace gas transport [Waddington et al. 1996].According to Oldenburg and Under [2003] barometric pumping has a negligible effect on the time-averaged seepage flux and near-surface CO 2 concentration because of the cyclic nature of the pressure-induced flows.
From Oldenburg and Unger [2003], the unsaturated zone can attenuate until 96% of CO 2 after 100 years with a leakage rate of 4•10 4 kg•yr -1 , while with a rate of 4•10 6 kg•yr -1 the attenuation efficiency decrease to 19%.The attenuation efficiency of the unsaturated zone decreased with increasing leakage rate because the higher pressure surrounding the source zone caused more vertical migration of the CO 2 relative to lateral migration, which is more strongly affected by infiltration.According to Oldenburg and Unger [2003] for a CO 2 plume present in the unsaturated zone with no continuous replenishment, dissolution into infiltrating water causes relatively rapid attenuation.
Comparing Solforata di Pomezia data with Oldenburg and Unger   respectively 74675 and 40000 ton•yr -1 .That flux is recorded by Oldenburg and Unger [2003] with a leakage rate of a bigger order of magnitude (4000000 ton•yr -1 ).Water circulation plays a great role in the Solforata di Pomezia site but is not always present, as explained before in the summer period.Atmospheric pumping can favour the emission especially in the winter, due to the changes of high and low atmospheric pressure; while water content favour dissolution in the winter period due to amount of water circulating.
Comparing leakage rate used in Oldenburg and Unger [2003] with Solforata di Pomezia data it can be stated that Solforata di Pomezia has a low leakage rate so barometric pumping effect is negligible.
According to Norstadt [1984], Xu and Qi [2001] and Gasparini et al. [2015], an increase of ambient temperature increases CO 2 soil fluxes intensity.Instead, for this case study, an inverse correlation is found between CO 2 soil flux and temperature (Figures 9 and 10): the plot and the multiple regression model show that both air T and the horizontal component of wind speed are anti-correlated with carbon dioxide air concentration, and is clearly visible how CO 2 concentration increases when temperature and wind speed drop.

Conclusions
The comparison between experimental data and numerical calculations of gas dispersion for the CO 2 natural emission site of Solforata di Pomezia shows a good reasonable fit between both sets of data.The predicted data comply with experimental data both for air concentration and plume size.This suggests that numerical modelling of atmospheric dispersion of CO 2 plumes emitted from underground can be a valuable tool to develop hazard maps in natural emission sites or in areas affected by underground storage.
Atmospheric dispersion maps predicted from numerical modelling using sampling gas data recorded in 2007 do not show hazardous conditions at 25, 100 and 200 cm from surface.
Modelling of critical scenarios of large gas leakage in the same site (five and ten times larger than the present-day emission) does not predict lethal gas concentration at least for human beings at heights 1 and 2 m above the ground, concentrations are largely under 3 vol.%.Air temperature reduces or enhances CO 2 diffusion from soil, yielding respectively lower or higher air concentrations.The presence of wind is important because it favours dispersion avoiding CO 2 accumulation at dangerous levels.
Future study to conduct is develop a new soil flux survey and record meteorological data during different seasons in order to investigate climate effects on soil emission.Carry out a CO 2 air concentration survey at different heights to corroborate numerical results with experimental data.

Figure 1 .
Figure 1.Aerial view of the Solforata di Pomezia, along its main axis.Picture taken from Google Maps.The red segments outline the area mined in the past for sulphur; the green circle shows an episodic flooded area.

Figure 2 .
Figure 2. Tectonic map of the Tyrrhenian margin of the Central Italy Peninsula [from Acocella et al. 1999].The location of Solforata di Pomezia is highlighted by a red dot.

Figure 3 .
Figure 3. Various views of the Solforata di Pomezia site, vegetation, lakes and dry soil.

Figure 4 .
Figure 4. Example of gas emission point source at the Solforata di Pomezia in the summer period [from Carapezza et al. 2005].

Figure 5
Figure 5 shows the CO 2 soil flux map resulting from measurements and the location of 7 TDL profiles, obtained in February and May 2007, by Carapezza et al. [2012].For the present study, TDL profile on May 2007 (Table1) will be used for a comparison with soil flux data obtained in the same period.

Figure 5 .
Figure 5. CO 2 soil flux map in the Solforata di Pomezia from Carapezza et al. [2012].Black dots are the 278 sampling points and numbered segments are TDL profiles.White area denotes the lake, where measurements could not be performed.

Figure 8 .
Figure 8. Wind speed recorded at the Solforata di Pomezia during the days May 14 th and 15 th 2007 by Carapezza et al. [2012].
Simulation features of the four performed models.

Figure 11 .
Figure 11.Plume dispersion at 25 cm height for S1 after 2.5 hours of simulation, for S2 after 4.5 hours, for S3 after 1 hour and for S4 after half an hour.

Figure 12 .
Figure 12.Prediction of the gas plume dispersion in model S1 at 25 cm height for 6, 12 and 18 h of simulation time.

Figure 13 .
Figure 13.Prediction of the gas plume dispersion in model S2 at 25 cm height for 6, 12 and 18 h of simulation time.

Figure 14 .
Figure 14.Prediction of the gas plume dispersion in models S3 and S4 at 25 cm height after 18 h of simulation time.

Figure 15 .
Figure 15.Predicted CO 2 air concentration for model S1 outputs at 1 metre and 2 metres high.
[2003]  results, it enhances the equilibrium reached in the Italian volcanic area.The Solforata di Pomezia has a soil flux of 1500 g•m -1 d -1 with a leakage rate double than Oldenburg and Unger[2003],

Figure 16 .
Figure 16.Predicted CO 2 air concentration for model S3 outputs at 1 metre and 2 metres height.

Figure 17 .
Figure 17.Predicted CO 2 air concentration for model S4 outputs at 1 metre and 2 metres height.