2018-volcanic-ash-plume
Authors:
Christoph von Matt, 13-105-358, [email protected]
Benjamin Schuepbach, 14-100-564, [email protected]
This project is our homework submission for the seminar 438745-FS2018-0- Geodata analysis and modeling at the University of Bern in Switzerland. We built a model of the ash plume dispersion after the eruption of Eyjafjallajökull in 2010.
To build our model we worked with Python v2.7.0 by using a Conda Environment.
The model dynamics is inspired by the many different particle-transport-models that exist and which often base on the advection-diffusion equation (as described in e.g. Folch 2012). We decided however to implement simplified transport-diffusion dynamics. The transport of the particles is determined by wind direction and wind speed. Therefore, we classified the surrounding cells according to the wind-direction. Depending on the wind speed, different percentages of the cells concentration are transported.
For further information on the transport and diffusion process, please consider the explanations below.
Our model can be run with different settings:
Test or Simulation: Choices related to wind field initialisation
The Test mode initializes constant wind fields with constant wind speed and wind direction. In this mode the wind field does not change during the whole modeling process. As indicated this mode is primarily for (functionality) testing purposes. Wind speeds for U-wind and V-wind components can be specified by the user. As the wind fields are created artificially, the user can specify the model resolution too.
The Simulation mode initializes wind fields out of NetCDF-wind files provided by the user. The user has to provide two NetCDF-files, one for each wind component. It is very important that all variables, except the ones for U- and V-wind components are identical!
Attention: The only supported data format is the NetCDF-format!
Manual Parameterisation or Eyjafjallajökull Parameterisation: Choices related to eruption characteristics
The Manual mode allows the complete specification of all parameters needed to calculate the erupted ash concentration. This includes the geographic location of the volcano, the plume height, the durance of the eruption, the mass fraction of the particles (<63 micrometers), as well as tephra mass and tephra volume rates. This results in one single concentration value which is erupted for the specified durance. Again, this mode is primarily for testing purposes.
In the Eyjafjallakökull mode the erupted ash concentration is calculated with predefined values for all needed parameters. We therefore use literature based values for the Eyjafjallajökull eruption event back in 2010. The result is an available eruption concentration sequence of 156 values (6-hourly-resolution) starting from the 14th April 2010 until the 22th May 2010.
We obtained the Eyjafjallajökull 2010 concentration by using the following literature:
- Plume height, eruption durance, (tephra) mass rate, (tephra) volume rate → (Gudmundsson et al. 2012)
- Ash mass fraction → (Mastin et al. 2009)
In the following all needed input parameters are listed and shortly explained. Depending on the chosen model settings the user may be prompted to provide input values.
Required parameters and components are:
-
Wind Data
In the Simulation mode the user has to provide 2 NetCDF-filenames (ending in ".nc") for both U-wind and V-wind components. It is very important that all variables, except the ones for U- and V-wind components are identical!
In the Test mode the user can specify the constant wind speeds for both wind-components.
In both cases the initialised wind fields consist of two rasters with one containing the values for the U-wind component and the other for the V-wind component.
-
Timesteps
Time is defined as steps, which in turn represent iterations.
In the Simulation mode the user can choose for which available time period (in days) the simulation should run. The amount of days is calculated from available timesteps and temporal resolution. The number of timesteps is then adjusted according to the chosen period length.
In the Test mode the hourly resolution is by default 1 hour. Thus, the user input equals directly the timesteps executed.
-
Resolution
-
Spatial Resolution
In the Simulation mode the user has to specify the spatial resolution such that it equals the one of the provided wind datasets.
In the Test mode the spatial resolution can be chosen arbitrary. To not obtain non-sense results both, the resolution in kilometers as well as the resolution in degrees should be chosen interdependently.
-
Temporal Resolution
In the Simulation mode the user has to specify the temporal resolution (in hours) such that it equals the one of the provided wind datasets.
In the Test mode the temporal resolution cannot be changed and is by default 1 hour.
-
-
Location of Volcano
In the Manual mode the user can specify arbitrary longitude and latitude values for the volcano (point source) location.
In the Eyjafjallakökull mode the location of the volcano is specified as follows:
→ Longitude: -19.625
→ Latitude: 63.625.Attention: Be aware that the model uses decimal coordinate values!
-
Particle fall-out
The particle fall-out parameter describes what amount of particles falls out of the system by each timestep.
The particle fall-out can be specified by the user in all model settings and should be given as 1 - percent value [0-1].The default fall-out was determined by the idea that the finer particles can remain in the air for several days without major fall-out. The default fall-out assumes that this durance is 6 days.
-
Diffusion
-
Diffusion coefficient
The diffusion coefficient was calculated according to a Gaussian plume model and was calculated as follows (assuming a neutrally stable atmosphere):
Diffusion coefficient = 68 x resolution^0.894 → (Cimbala 2018) -
Diffusion Percentage
The diffusion percentage determines how much of the ash concentration is diffused.
For the calculation of the default diffusion percentage the diffusion coefficient is used. As the Gaussian Plume model may not be accurate for synoptic scale dispersion / diffusion the user can also choose an arbitrary diffusion percentage in the range between 0 and 1.
-
Diffusion Type
We implemented two different diffusion types in our model (see Section "What our Model does (How to Volcano)").
The choices are as follows:
1 - Diffusion in all directions (Choice: 0)
2 - Gradient-dependent diffusion (Choice: 1)
3 - no diffusion (Choice: any other number)
-
-
Ash Plume height
The ash plume height specifies the height of the erupted ash cloud in kilometers.
As already declared above, in the Manual mode the user can specify the ash plume height, while in the Eyjafjallajökull mode predefined values (according to Gudmundsson et al. 2012) are used.
-
Mass- and Volume Rates
The mass- and volume rates specify the erupted rates of tephra mass (in g/m^3) and tephra volume (in m^3/s).
As with the plume height, the user can specify the mass- and volume rates manually in the Manual mode while predefined values are used in the Eyjafjallajökull mode (according to Gudmundsson et al. 2012).
-
Ash Fraction
The ash fraction defines the percentage of particles smaller than 63 micrometers.
The erupted tephra mass is multiplied by this factor to obtain only the fine ash fraction, which potentially can be transported over large distances.
In the Manual mode the user can specify the mass fraction value arbitrary between 0 and 1.
In the Eyjafjallakökull mode the ash mass fraction is predefined as 0.5 (according to Mastin et al. 2009).
-
Ash Concentration
The erupted ash concentration is calculated as follows:
concentration = mass rate * ash fraction / volume rate / resolution^2
The final division with the squared resolution accounts for the issue that the spatial resolution of the model is almost all cases much larger than volcanic eruption area (point source!). The consideration of this areal difference therefore slightly reduces the ash concentration in the model.
After the initialisation of all required input parameters the model performs more or less raster calculations.
The two most important rasters are:
- U- and V-wind speed rasters
- Particle concentration raster
The particle concentration raster stores the different ash concentration distributions for each timestep.
1) Eruption is happening!
By its initialisation the particle concentration raster is empty. Thus, before transport and diffusion can happen an eruption has to happen. This is realized by adding the calculated concentration to the raster cell which is closest to the specified geographic location of the volcano.
2) Transport of particles
After the eruption has happened, the ash can be transported by the wind.
The transport depends on both wind speed and wind direction.
For each cell of the particle concentration raster the same procedure is performed:
2.1) First, the wind direction is calculated. This is done by calculating the angle resulting from both components.
2.2) Second, the maximal wind speed between the individual components and also the resulting wind speed of both wind components together is calculated.
The transport mechanism works as follows (see Figures below):
Say the model is processing the cell of the particle concentration raster with the given indices i (row) and j
(column).
Figure 1:
Each surrounding cell is classified in the model starting by the index i - 1, j (Cell 1) proceeding clockwise till
i - 1, j - 1
(Cell 8). Further, the cells are classified with respect to the 360° circle. We therefore divided the
360° by the 8 surrounding cells, resulting in a classification where each surrounding cell falls between a 45° range.
Figures 2 and 3:
Figure 2 shows an example wind field with constant wind field with only U-wind components as it could typically be running the model in the Test mode.
The calculated wind direction angle in this example would be 90° as there is only a U-wind component present in the two wind rasters. By the means of this angle and the explained cell classification (see Figure 1) the receiving wind transport direction cell would therefore be Cell 3 (see Figure 3).
As mentionned before, the particle transport depends not only on wind direction, but on wind speed too. That is were the maximal wind speed and the spatial model resolution comes into play:
Let's assume that the spatial model resolution is 80 km. If the calculated maximal wind speed exceeds the spatial model resolution (here 80 km/h) then 100% (1.0) of the ash concentration in the Cell ( i , j ) is transported to Cell 3. Hence, if the maximal wind speed is lower than 80 km/h, then only a predefined proportion of less than 100% (< 1.0) of the ash concentration will be freighted to Cell 3.
That is actually not the whole story. What is with the diffusion ash concentration part?
The amount of ash concentration dedicated to diffusion is subtracted before the transport process takes place.
So, in the example only 1.0 * (total ash concentration - diffusion concentration amount) is transported to Cell 3.
Attention: Please remember that in the model the U- and V-wind component are stored in separate rasters!
3) Diffusion of particles
After the particle transport mechanism has been completed, the diffusion part begins.
Our model provides two different diffusion types: one in all directions and a gradient-dependent one. The amount of ash concentration which is diffused, depends on the initially defined diffusion percentage.
Continuing the above example, the Cell i , j contains now the ash concentration of the diffusion part (which was not
transported) and a possible transported part from neighbouring cells.
The current concentration in Cell i , j is thus:
current concentration = 1 - 1.0 * (total ash concentration - diffusion concentration amount) + transported concentration
which is equal to
current concentration = (remaining) diffusion concentration amount + transported concentration (in this cell).
Before diffusion can take place - just like in the transport mechanism - the diffusion concentration amount has to be specified again. The available ash concentration for diffusion is therefore current concentration * diffusion percentage.
Let's now look at the different diffusion types!
3.1) Diffusion in all directions
In this diffusion type the diffusion concentration amount is equally distributed to all 8 surrounding cells. Each cell
receives therefore an ash concentration of diffusion concentration amount / 8 (see Figure 4).
3.2) Gradient-dependent diffusion
The second diffusion type is slightly more complex.
As a first step, the concentrations of all surrounding cells are retrieved. Then the
difference with respect to the Cell i , j is
calculated. Then, these differences are divided by the distance separating each cell from the Cell i j.
This is done with the spatial model resolution. The differences then become gradients.
An example calculation for Cell 3 looks like:
gradient = (concentration Cell 3 - concentration Cell i, j ) / 80 km
Note: The diagonally situated cells 2, 4, 6, 8 require a corrected distance. This is realized using simple trigonometry.
In the next step the available ash concentration for diffusion is equally distributed to the cells showing a negative gradient with respect to Cell i , j (see Figure 5).
Thus each receiving cell gets the following concentration:
diffusion concentration amount / Number of cells with negative gradient.
With this diffusion type, the diffusion distributes the diffusion concentration amount only towards lower ash concentrations.
4) Particle Fall-out
As last step, before the whole transport diffusion loop starts all over again, the particle fall-out takes place.
To model the particle fall out, the initially specified fall-out percentage is applied over the whole
particle concentration raster. Every cell of this raster looses the same percentage of ash concentration with each
timestep. Thus, cells with higher concentrations lose a proportionally higher amount of ash particles.
Finally, our model outputs three fully drawn maps for each iteration of the transport and diffusion loop. So in the end, the model results can easily be displayed as a GIF file.
To guarantee that the model does correct calculations and that no mass is created or disappearing mysteriously, we implemented a control mechanism in our model.
Principle:
With each timestep the fall-out is summed-up and stored. After the transport and diffusion procedure is complete, all
concentrations of the very last particle concentration raster are summed-up and added to the summed-up fall-out.
Simultaneously, each eruption of ash concentration is summed-up as well.
In the end sum of eruptions has to be equal to the sum of Sum final particle raster + Sum fall-out.
As mentionned in the "What Our Model does (How to Volcano)" section, our model generates three different maps. All maps are realized with the Basemap package.
The colouring and spacing of the colorbars is strongly inspired by the EUMETRAIN's Volcanic Ash Training Module.
All the specified (and relevant) parameters used in a specific model run are visible on the individual plots and are updated automatically if different model settings and input parameters are used.
For the Eyjafjallajökull simulation mode we use 500hPa ERA-Interim wind data of April 2010 which we directly obtained from the ECMWF website.
In this section we present some of our results we got from our model runs.
A validation was not yet performed and thus the adequacy of the model for the Eyjafjallajökull eruption 2010 will not be
discussed here.
A possible validation strategy could be a comparison with other transport-dispersion models, such as the visualisation
on the EUMETRAIN's Volcanic Ash Training Module.
This figure shows the model simulation on a global map.
The parameters for this model run are specified on the plot.
This figure shows a version of Figure 1 with limited extent. Here default extent of this figure is Europe.
With this figures we try to give an impression of how the individual components (transport, all directions diffusion, gradient-dependent diffusion) work.
The last plot generated by our model accounts for the ash concentration risks for air traffic. We used here the concentration breaks defined by the Civil Aviation Authority (CAA).
We are fully aware that our model has several limitations and issues which are not solved yet. They shall be discussed in this section.
Limitations related to the model dynamics consists of issues related to spatial and temporal model resolution as well as to fall-out and dispersion.
The first dynamics-related problem considering the model resolution is, that with each timestep the maximal particle transport distance is limited to the spatial resolution of the model. No matter whether maximum wind speeds are much higher than the distance to traverse a model cell, the particles will maximally be transported as far as the model resolution.
Example:
With a spatial resolution of 80 km and maximum windspeeds of 160 km/ h the particles won't be transported further than 80 km although the potential transport is twice as far.
Another issue is the temporal resolution. The temporal resolution of the wind fields we used in the "Results" section is
six-hourly. This means we have only one new wind field every six hours.
Using one wind field per timestep would however not make sense as the spatial issue even aggravates as the particle transport
distance is heavily underestimated.
Our model attempts to reduce this temporal resolution issue by iterating with the unchanged wind fields for the durance of the model resolution.
In addition to model resolution related issues, there are also limitation with the fall-out parameterisation
First, the ash can sometimes stay in air for several days without any major fall-out. This is not considered in our model as fall-out begins from the first iteration on. As we model only one specific level we neglect this issue.
Another fall-out issue is its constant parameterisation. Although we limited our simulation to particles smaller than 63 micrometers, the ash particle size can be still quite diverse. The constant fall-out does therefore not account for the fall-out velocity of different particle size fractions. This is clearly a limitation of our model.
General issues consider some code functionality problems. These are not discussed in this documentation.
As our model operates almost exclusively with raster data we can imagine an implementation of this model in a geographic information system. A combination with air space territories could be worth a try to detect flight zone restrictions due to ash concentrations.
Thanks go to our instructors and supervisors during the semester, Dr. Andreas Zischg (University of Bern), Dr. Jorge Ramirez (University of Bern) and Dr. Pascal Horton (University of Bern) for their strong commitment to make the "Geodata analysis and modeling" lecture in its first realisation a full success. Thanks for the very informative and interesting lecture!
Cimbala, J. M. (2018): Gaussian Plume Model and Dispersion Coefficients. Pennsylvania State University. Latest revision: 30th January 2018. Available under: https://www.mne.psu.edu/cimbala/me433/Lectures/Tables_for_Gaussian_Plume_Model.pdf
Folch, A. (2012): A review of tephra transport and dispersal models: Evolution, current status, and future perspectives. In: Journal of Volcanology and Geothermal Research 235-236: 96-115. DOI: 10.1016/j.jvolgeores.2012.05.020.
Gudmundsson, M.T., Thordarson, T., Hoeskuldsson, A., Larsen, G., Bjoernsson, H., Prata, F.J., Oddsson, B., Magnusson, E., Hoegnadottir, T., Petersen, G.N., Hayward, C.L., Stevenson, J.A., Jonsdottir, I. (2012): Ash generation and distribution from the April-May 2010 eruption of Eyjafjallajoekull, Iceland. In: Scientific Reports 2:572: S.1-12. DOI: 10.1038/srep00572.
Mastin, L.G., Guffanti, M., Servranckx, R., Webley, P., Barsotti, S., Dean, K., Durant, A, Ewert, J.W., Neri, A., Rose, W.I., Schneider, D., Siebert, L., Stunder, B., Swanson, G., Tupper, A., Volentik, A., Waythomas, C.F. (2009): A multidisciplinary effort to assign realistic source parameters to models of volcanic ash-cloud transport and dispersion during eruptions. In: Journal of Volcanology and Geothermal Research 186: 10-21. DOI: 10.1016/j.jvolgeores.2009.01.008.