Seb's repository

TephraProb v.1.7

2020-02-21T00:00:00+00:00

TephraProb v1.7 includes access to ERA-5 wind data, a different sampling of masses and some bug fixes.

Changes in TephraProb v1.7

Added functionalities

ERA-5

ERA-5 is now available for download. It follows the same procedure as ERA-Interim (refer to manual), but it uses a different Python API (cdsapi) - which also requires a new key.

Create a new account on the CDS website
Retrieve your UID and key and format it as such (where 1234 is your UID and the string following the column is the API key):
```
url: https://cds.climate.copernicus.eu/api/v2
key: 1234:abcdefghij-134-abcdefgadf-82391b9d3f
```
Enter this key in TephraProb under Input > Wind > Set ERA-5 API key
Install the cdsapi in TephraProb using Input > Wind > Install ERA-5 libraries

Note it is only required to perform steps 4-5 once.

Relevant ERA-5 documentation

Reference page for ERA-5
Good introduction to ERA-5 on retostauffer.org.
Introduction to CDS API

Loading wind datasets

It is now possible to load existing datasets in the Input > Wind > Download Wind module.

Logarthmic sampling of mass

For subplinian/Plinian scenarios, the sampling of the eruption mass when constrain=0 was linked to the distribution chosen for the plume height. A new variable mass_sample was introduced to control the shape of the distribution from which the mass is sampled:

mass_sample=0: the mass is sampled between min_mass and max_mass
mass_sample=1: the mass is sampled between log10(min_mass) and log10(max_mass)

Linking outputs to original ESPs

A new scheme was introduced to be able to link the content of the files dataProb.mat and dataT2_*.mat located in RUNS/runName/runNb/DATA/. This is explained in more details in

Bug fixes

Fixed a bug on the sampling wind profiles when seasonality=1

TephraProb troubleshoot

2018-12-20T00:00:00+00:00

This is an attempt to track issues that could not be solved within the code. Keep in mind that having a code that works on everybody’s computers is hard, especially on Windows where errors can arise because of loaaaaads of different factors. Also, please keep in mind I am not a personal hotline :)

Content

ECMWF-Related issues
- Pip
  - Case yey!
- SSL
- Setuptools
- ECMWF Python API
Tephra2 related issues
- Mac OS 10.14 Mojave

Pip

Pip is a package manager, see it as an AppStore for nerds, but the principle is the same. Most of the errors related to the access to ECMWF datasets can be solved by installing small bits of code, and Pip helps with that. The first thing you have to check is whether Pip is installed on your system. To do so, you need to open a terminal

Windows: Open the start menu and type cmd and open Command Prompt
Mac OS: Open Terminal.app
Linux: I assume you know how to do it

Now type:

pip

Case ney :(

Do you see anything that tells you that pip is not a recognized command? Then you need to install it. First, still in the terminal, let’s hope that the following command works

easy_install pip

Else, you unfortunately need to figure it out using this tutorial

Case yey!

Move on to installing the required package as described below!

SSL

Error:

urllib.error.URLError: <urlopen error [SSL: CERTIFICATE_VERIFY_FAILED] certificate verify failed: self signed certificate in certificate chain

Fix:

pip install pyopenssl

Setuptools

Error:

ImportError: No module named setuptools

Fix:

pip install setuptools

ECMWF Python API

Summary of troubleshooting issues here.

Mac OS 10.14 Mojave

If error during compilation, this thread helps fixing it.

Random tools

2018-12-13T00:00:00+00:00

Just a repository of cool and handy tools.

Github repositories

Screengif: Software to convert screencasts to gifs
Pandoc: Convert text files between a lot of formats

Python Resources

2018-12-04T00:00:00+00:00

Keeping track of interesting resources on Python.

Tutorials

Geology and Python: Cool tutorials for Python applications to Earth sciences
Open Geo Tutorials: Intro to GIS, remote sensing and bits of machine/deep learning

Github repositories

pyPhenology: Phenology toolbox
seasonal: Robustly estimate trend and periodicity in a timeseries
geog0133-practicals: Practicals on phenology modeling (especially this)

GEE Resources

2018-11-30T00:00:00+00:00

Keeping track of interesting resources on the Google Earth Engine.

Example scripts

Use of the 25th percentile to make cloud free mosaics link
Merging of Landsat products link

Toolbox

GEE implementation of the LandTrendr spectral-temporal segmentation algorithm link
GEE tile layers with Folium link

Tutorials

MyGeoBlog: Great tutorials

Papers

Roy et al. (2016): Scaling between Landsat 7 and 8

Tips

Assign value to masked regions link

Empty cells in table export [link](https://groups.google.com/forum/#!searchin/google-earth-engine-developers/reduce$20collection$20empty

sort:date/google-earth-engine-developers/qsoHZ9Dhqlw/kJncwcXUEAAJ)

Maps

Intro to Folium (french) link

Find wind profiles most similar to an average

2018-10-31T00:00:00+00:00

Ever needed to find the one wind profile most similar to an average wind profile? A new function in the Analyse wind profiles of TephraProb does just that!

Access the tool

The tool is available from Tools > Find similar profiles in the Analyse wind profiles window. Some specific parameters should be set for the code to work, namely:

Type must be Profiles
Average must be Averaged
Direction must be Direction

Method

The similarity between one wind profile and an average profiles is assessed following these steps:

The root mean square error of both direction and velocity is calculated between all wind profiles in the database and the average wind profile across all heights
The euclidian distance between the direction and the velocity is minimised
The similarity is expressed in terms of decreasing euclidian distance

Use the tool

Upon clickling the tool, keep in mind that the average data used in the process is the data that is plotted in the axis. So should you select the average over a specific range of months, the most similar wind profile will reflect that range.

The first option popping up allows to choose a range of elevations over which the previous calculation is applied. Having a plume height of 15 km asl above a volcano located at 5 km asl, the entire range of altitudes provided by reanalysis databases might not be relevant. This option allows to target where the similarity matters,

The next window that opens contains the 100 most similar wind profiles ordered from top to bottom in decreasing similarity. In other word, the first item in the list is the best fit. Each wind profile is named with the date - remember that the hours are UTC - and the corresponding file number as named by TephraProb. You can now explore which single profile is best for you by plotting one or multiple profiles.

TephraProb v.1.6

2018-10-02T00:00:00+00:00

There has been a change in the access to the Google Map API used by the function plot_google_map by Zohar Bar-Yehuda that can cause problems when TephraProb attempts to retrieve the map background. Google now requires to get an API Key, which needs to be associated with a Google Project and, tadaaaa, a billing account. Although this is likely to remain free for a small number of requests, some might be reluctant to add a billing account to get a key. A new function offering to export maps to Google Earth was therefore implemented, which relies on the Google Earth Toolbox by Scott Lee Davis.

Changes in TephraProb v1.6

TephraProb preferences

The preference pannel has been reworked. Regarding Google Maps, the options are:

Choose whether a basemap should be plotted. The option for Google Map remains, knowing that problems may occur if not used with an API Key;
Enter an API Key for Google Map, which is likely to solve the problems.

Regarding Google Earth, the options are:

Choose to export the set of plotted data to Google Earth. If selected, this option saves a .kmz file in RunName/RunNb/KML/;
Choose to display the .kmz upon completion of the generation of the maps. Note that this process will attempt to open the .kmz file with the default software assigned in your system for this file extension. If your computer opens another software, you must manually change the default software used for .kmz files.

Google Earth output file

If the option described above is selected, TephraProb will generate a Google Earth file upon completion of the Display / Probability maps and Display / Isomass maps functions. The Google Earth file generated by TephraProb is a KMZ file, which is an archive containing all the associated files required for display that can easily be shared amongst people. Once opened, the file contains:

Each map selected for plotting in TephraProb;
A color bar;
The extent of the original computation grid. Note that some re-interpolation is necessary to plot the maps in Google Earth and some data may appear outside the computation domain. This is an artefact and should be considered as such;
The vent.

Additionaly, if i) hazard curves have been computed for the working project , ii) the option Plot hazard curve location is activated in the preferences and iii) probability maps are being plotted, TephraProb adds the location of the hazard curves in the .kmz file. These points are clickable and will display the hazard curves.

Added dependency

Function Google Earth Toolbox by Scott Lee Davis.

Wind data for Tephra2

2018-07-27T00:00:00+00:00

I thought I’d share some problems I’ve recently encountered with the format of wind files used in Tephra2. Just as a reminder, wind file should be a 3-columns, tab-delimited ASCII file containing:

Column 1: Altitude (m asl)
Column 2: Wind speed (m/s)
Column 3: Wind direction (i.e. the direction the wind blows to; degrees from N)

Let’s review some aspects that can cause problems in Tephra2.

Null elevations

Tephra2 no like no 0 m elevation. Period. The reason is that wind directions and velocities are interpolated between levels. Entering a minimum elevation of 0 m asl will cause an error in the interpolation, which will not be output but will produce unrealistic results

End-of-the-file

It should not matter whether wind profiles are organised with ascending or descending altitudes as Tephra2 re-orders it anyway. BUT, because there is a but, one thing that can cause a problem is the definition of the End of the File - or EOF, which is:

…a condition in a computer operating system where no more data can be read from a data source.

This EOF can have an impact when i) wind files are organised in descending elevations and ii) the EOF is a new line (i.e. check if a file has an empty line at the end). If these two conditions occur, Tephra2 will reorder the wind data in ascending order, which will result in empty data as the first line and all hell will break loose.

Conclusion

The two main important points in the format of wind files for Tephra2 are:

Wind files should not contain elevations of 0 m asl. Enter 1 instead
Avoid empty lines in wind files, or at least make sure wind files are organised in increasing elevations

How to notice problems?

As already specified, should these problem occur, Tephra2 will not output an error but will produce unrealistic results. From experience, bogus results tend to display a conspicuous eastwards dispersal. So always double check the model output against the wind data, which is valid for both Tephra2 and TephraProb!

TephraProb v.1.5

2018-07-27T00:00:00+00:00

Version releases were never systematically attributed to TephraProb until now. This new version was labeled v1.5 as some important changes have been made.

Changes in TephraProb v1.5

Probability calculations

Changes on storage of probability data have been made to reduce disk space usage and improve speed of reading Tephra2 outputs during probability calculations. Namely, the probability_maker function, accessed from Post processing > Probability calculations in TephraProb, proceeds in 2 steps:

Read the Tephra2 output, sums them if needed and stores the output in RunName/RunNb/DATA/dataT2_*.mat. This is a 3D matrix, where the 3rd dimension equals the number of runs (i.e. each summed Tephra2 output is stored along the 3rd dimension). The 1st and 2nd dimensions vary with the grid type:
- If runs are performed on a full grid, dimensions 1 and 2 are equal to the dimensions of the grid northing and easting, respectively
- If runs are performed on a list of points, each output is in a 3 column format, with dimension 1 equal to the number of points and dimension 2 being [easting, northing, mass accumulation]
This new output is used for probabilit calculations, which are stored in RunName/RunNb/DATA/dataProb.mat. This file loads a Matlab structure containing the outputs of probability maps, probabilistic isomass maps and hazard curves

Note that to save disk space, the output of Tephra2 is now:

Rounded to a significant digit number
Accumulations below a thresholds are considered equal to zero

Both the significand digit number and the thresholds are controlled in the preferences of TephraProb, accessible from File > Preferences.

Export to ASCII files

ASCII files previously generated upon completion of the post-processing stage are not automatically generated anymore, which save a considerable amount of disk storage (i.e. the 3C/, GIS/ and MAT/ foldeds in the SUM/, PROB/ and IM/ folders). Instead, a new function accessible from Post processing > File management > Export to ASCII files allows the user to choose the type and the format of files to export.

Archive Tephra2 files

A new function available from Post processing > File management > Archive Tephra2 files offers to compress files generated by and for Tephra2 (i.e. the OUT/, CONF/ and GS/ folders). This was also implemented with concern for disk space usage. Note that this function should only be used once the Tephra2 modelling has successfully completed. This process takes a long time due to to the large amount of small files and individual folders created.

Probabilistic isomass maps

The method used to compute probabilistic isomass maps has been re-written and uses a percentile approach rather than an interpolation of iso-probability contours on probability maps. Note that a probability of occurrence of a given isomass is given by the inverse of the percentile, i.e. a percentile of 25% reflects a 75% probability of occurrence.

Hazard curves

It is now possible to import text files in the Input > Points dialogue box. Text files should be 3-columns tab-delimited ASCII files with columns containing:

Column 1: Point name
Column 2: Latitude (decimal degree, WGS84)
Column 3: Longitude (decimal degree, WGS84)

Map preferences

A more in-depth customization of output maps is now offered from File > Preferences. These parameters allow to:

Scale the colors of probability and isomass maps to the min/max values set in their respective contours. For instance, this is useful for plotting isomass maps for various probability threshold while preserving the value of a pixel color across maps
Control the display of contours labels. This is useful when the display gets jammed
Plot the points used for the calculation of hazard curves and their names
Chose the colormaps
Plot probabilistic isomass maps with a logarithmic colormap

Bug fix

A bug was corrected for Plinian scenarios, causing problems in the calculation of eruption duration for eruptions > 6h and the long-lasting option disactivated
A problem was reported by two users concerning an improper interpretation of ECMWF wind files in Tephra2. The problem is discussed here. Changes were made in TephraProb to avoid the occurrence of this bug, but please check any run that could have been affected

TephraProb on a cluster

2018-06-08T00:00:00+00:00

There is a relatively easy way to parallelise TephraProb on a computer cluster without having to struggle with Matlab’s Parallel Computing Toolbox. Matlab is not even needed on the cluster, and it is just a matter of sending single occurrences of Tephra2 on different nodes. Here is how to do so:

Generate your input files and eruption scenarios locally
Send the required files on the cluster
Run the scenario on the cluster
Retrieve the output files from the cluster to the local computer
Post-process the output files (e.g. probability calculations) locally

Setting up remote files

Following the procedure in the TephraProb manual, all tasks should be completed until section 5.3 - meaning you should be set to hit the Run Tephra2 function. The main file that will be used for the parallelisation is T2_stor.txt, located in RUNS/runName/runNumber/T2_stor.txt, which contains the Tephra2 commands for all single model runs of the scenario.

Transfer your run, grid and wind files on the cluster. Unless you decide to customize T2_stor.txt, the directory tree once on the cluster should look like the following. Note that not all files of the RUNS/, GRID/ and WIND folders need to be transferred.

ROOT
├── MODEL/
├── RUNS/
│   └── runName/
│       └── runNumber/
│           ├── CONF/*.*
│           ├── GS*.*
│           └── OUT/*.*
├── WIND/
│   └── windName/
│           └── pathToAscii/*.*
├── GRID/
│   └── gridName/
│           └── *.utm
├── T2_stor.txt
└── runTephraProb.sh

Note that T2_stor.txt should be at the root, and all commands specified in it should point to locations that have been uploaded on the cluster. Note that by default, all paths defined in T2_stor.txt are relative.

The rsync becomes handy here. To copy RUNS/ folders, you can use something like the following command, which will ignore all unecessary files and folders:

rsync -arvz --exclude *.mat --exclude FIG --exclude KML --exclude LOG --exclude SUM  RUNS/run_name host@server:~/TephraProb/RUNS/

Compiling Tephra2

Tephra2 needs to be recompiled for the cluster’s architecture. On the cluster, from the root of TephraProb, navigate to MODEL/forward_src/ and enter:

make

Now, navigate back to the root of TephraProb and type:

chmod 755 MODEL/tephra2-2012

That should get Tephra2 running on the cluster. Getting a File format not recognized error? See if that helps.

Running in parallel

The parallelization is achieved using job arrays. Conceptually:

T2_stor.txt
is cut in smaller files named T2_stor.txt00, T2_stor.txt01 … T2_stor.txtXX using the split Unix command
The subset of commands to Tephra2 in each of the sub-file is sent to a different node using a job array

On the cluster, split T2_stor.txt:

split -l 1000 -a 2 -d T2_stor.txt T2_stor.txt

where -l is the number of line of each sub-file (here T2_stor.txt will be split into subsets of 1000 lines) and -a is the number of digits appended to the name of the subfile (e.g. -a 2 produces 01, 02 etc…). The last argument is the generic name of the sub-file.

Let’s say that this created 10 files named T2_stor.txt00 to T2_stor.txt09. We need to adapt the job array to account for the range 0-9. We then use a handy little piece of code called GNU Parallel to send single CPU jobs to the nodes.

SLURM

On a SLURM cluster, the bash runTephraProb.sh might look like that:

module load GCC/4.9.3-2.25
module load OpenMPI/1.10.2
module load parallel
chunk=`printf "%02d" $SLURM_ARRAY_TASK_ID`
srun parallel -j 16 -a T2_stor.txt$chunk

The job can then be submitted using:

sbatch --array=0-9 runTephraProb.sh

OpenPBS

On an OpenPBS cluster, the bash runTephraProb.sh might look like that:

module load openmpi/1.4.5-gnu
module load parallel
cd $PBS_O_WORKDIR
chunk=`printf "%02d" $PBS_ARRAYID`
mpirun -np 12 -machinefile $PBS_NODEFILE parallel -j 12 -a T2_stor.txt.$chunk

The job can then be submitted using:

qsub -t 0-9 runTephraProb.sh

Post-processing

Once the modelling is finished, copy the remote version of RUNS/runName/runNumber/OUT/ on to your local RUNS/runName/runNumber/OUT/ and go on with the post-processing. Again, the rsync command can help you achieve something like that: The job can then be submitted using:

rsync -arvz --ignore-existing host@server:~/TephraProb/RUNS/run_name/ RUNS/run_name/

Tephra2 Inversion

2018-06-06T00:00:00+00:00

The current set of functions helps running the advection-diffusion model Tephra2 in inversion mode to estimate the best eruption source parameters (ESP) of a tephra deposit. The functions contain two sections. A first section contains a mixture of bash and python scripts for running the inversion of OpenPBS and SLURM clusters. A second section contains Matlab scripts for processing the inversion output designed to help its interpretation.

Updates - November 2018

The inversion scripts have been updated in November 2018. To find the scripts to post-process runs performed before this date, please download the previous version here.

Table of content

Introversion
- Single vs. batch runs
- Caveats
Installation
Getting started
- Case-study files
- Configuration file
Running the inversion
- Single vs. batch vs. seed runs
- OpenPBS
Post-processing
- Output figures
  - 1. Spaces of solutions
  - 2. Everything against everything!

Introversion

The inversion searches for the set of ESP that best reproduce observed values of tephra accumulations (kg/m²) with the computed values produced by Tephra2. As for any optimization problem, the many degrees of freedom require a critical interpretation of the inversion results. These scripts help building an empirical vision of the results in the perspective of the knowledge of the studied eruption.

For a purpose of illustration, let’s consider that the fit depends mainly on plume height and eruption mass. We will therefore attempt to find the combination of plume height and eruption mass that best reproduce the deposit. To do so, the Tephra2 uses a downhill simplex method, which searches for the minimum of a function - here the fit between observed and computed values - between user-defined ranges of [height_min - height_max] and [mass_min - mass_max].

Unfortunately, this method can suffer from numerical artifacts resulting in “fake” minima, especially when initial user-defined ranges are too large. For instance let’s consider two sets of plume heights and mass, height₁, mass₁ and height₂, mass₂, both giving a similarly good solution from a numerical perspective, but height₁, mass₁ being irrealistic from a volcanological point of view. We therefore need a method to subjectively discard this solution from the inversion result.

Single vs. batch runs

To avoid this problem, it is possible to run the inversion using three approaches, the single, the batch and the seed modes. The single mode searches for a minimum within the entire space defined by the ranges of plume heights and masses and consists of one single inversion run. The batch mode splits the space of plume heights and masses into smaller domains, and one inversion run is performed for each domain. The seed mode, for a lack of a better term, is pretty much the equivalent of the single mode repeated multiple times varying the seed. The seed is a number that controls the reproducability of stochastic runs, and in the case of the inversion it will control the starting point of the optimisation process. In case of unconstrained deposits, varying the seed is likely to affect the results. It is a good procedure to assess the sensitivity of the results to the input seed.

So what are the advantages and disadvantages of each method?

Single
- Pros: Provides the most accurate minimum
- Cons: Can get stuck in local minima
Batch
- Pros: Provides an overview of all minima on the height-mass space
- Cons: Small sub-spaces do not provide enough flexibility to accuratly define minima
Seed
- Pros: Advantages of the single approach without disadvantages
- Cons: Euuh… Not sure, it’s pretty much my prefered approach now

Caveats

Below are a few things to keep in mind when using the inversion:

The mass is always better constrained than the plume height
Constraining the plume height is more difficult because of the fact that in Tephra2, the plume height is intrinsically related to the computation of parameters controlling i) the grain size distribution, ii) the distribution of mass in the plume and iii) empirical parameters

Installation

Two components are needed, both available on GitHub:

Tephra2 source code, available here
Tephra2Utils, available here

For Tephra2, as far as we are concerned here, we only need the two executables created by the compilation of Tephra2 (tephra2012_inversion and tephra2-2012). For Tephra2Utils, we need the entire Inversion/ directory and the post processing files (plotBetaPlume.m, plotT2.m and processT2Inversion.m).

Organize your working folder like this:

ROOT
├── Inversion/
│       ├── _example_folder/                    -> Example of input files
│       │        ├── forwardGrid.utm
│       │        ├── inversionConfig.conf
│       │        ├── inversionInput.txt
│       │        ├── inversionWind.txt
│       │        └── runInversion.sh
│       ├── _scripts/                           -> Scripts used during inversion*
│       │        ├── genConfig.py
│       ├── _templates/                        -> Templates used during inversion*
│       │        ├── forwardConfTemplate.conf
│       │        └── inversionConfTemplate.conf
├── dependencies/                               -> Dependencies for post-processing*
│       └── xlwrite/
├── plotBetaPlume.m                             -> Used during post-processing*
├── plotT2.m                                    -> Used during post-processing*
├── processT2Inversion.m                        -> Main post-processing function*
├── plotT2Inversion.m                           -> Plot results of batch and seed runs*
├── tephra2012_inversion                        -> Tephra2 inversion exec
└── tephra2-2012                                -> Tephra2 forward exec

* Should not have to be edited

Getting started

Duplicate the folder Inversion/example_folder/ and rename it to whatever name makes you happy. I usually name this folder after the tephra deposit I am inverting. Note that you can have as many folders for as many deposits of as many volcanoes as you want here. Let’s now edit the required files.

Case-study files

Replace the following files by those for your own case-study.

inversionInput.txt:
This is the main file containing your field observations. It is a 4-columns, tab-delimited file containing:
- Column 1: Easting coordinates of the outcrop (UTM, WGS84)
- Column 2: Northing coordinates of the outcrop (UTM, WGS84)
- Column 3: Altitude (m asl). Note that Tephra2 can not accurately use the elevation. Instead, a mean elevation of all outcrops should be estimated and used as a unique value for all outcrops
- Column 4: The mass accumulation (kg/m²) estimated at each outcrop. Note that this is a mass accumulation and not a thickness
forwardGrid.utm:
The calculation grid used by the forward solution of Tephra2, and used during the inversion process to compare the deposit modelled by the inversion with field data. You can use TephraProb to create a grid following this tutorial.
inversionWind.txt:
If a wind observation is available for the studied eruption, it can be used to constrain the inversion. If no wind is available, ranges of wind speed and directions can be defined instead (described later), but please keep then the inversionWind.txt file in the folder. If the eruption is within the time interval provided by Reanalysis datasets, TephraProb can be used to retrieve it. The wind file should be a 3-columns, tab-delimited file containing:
- Column 1: Altitude (m asl)
- Column 2: Wind speed (m/s)
- Column 3: Wind direction (i.e. the direction the wind blows to; degrees from N)

Configuration file

The main inversion file is inversionConfig.conf, which must not be renamed. As a tip, replace the ERUPTION_NAME in the second line by the name of the eruption you are modelling. That my save you precious time when looking for data years from now!

Below is a summary of the different sections to fill in.

Inversion behaviour

This section controls the general behaviour of the inversion runs.

Variable	Description
`BATCH`	Enter 0 for single runs, 1 for batch runs and 2 for seed rund
`fixedWind`	Enter 0 to use ranges of wind direction/speed, 1 to use a wind profile
`SEED`	Used to control the seed for reproducibility. To randomely choose a seed, enter -1

Input files

This section defines the path to the input files previously described.

Variable	Description
`inputFile`	File containing the field observations
`windFile`	File containing the wind profile. Leave inversionWind.txt if fixedWind=0
`gridFile`	Path to the grid file

Vent properties

This section defines the vent geometry. Coordinates are in UTM and should be projected on the same zone as the field data.

Variable	Description
`ventE`	Vent easting (UTM, WGS84)
`ventN`	Vent northin (UTM, WGS84)
`ventA`	Vent elevation (m asl)
`ventZ`	UTM zone of the vent. Only the numeric part, negative in S hemisphere

Ranges to invert

This section contains all the variable parameters that are optimized during the inversion. Each parameter contains min/max values.

Variable	Description
`Ht`	Plume heights (m asl)
`Mass`	Total mass (log10, kg). Used only when BATCH=0
`Diff`	Diffusion coefficient for large particles (m²/s)
`FTT`	Fall-time threshold(s)
`MedPhi`	Median of the TGSD (Phi). Note that minMedPhi is the coarsest diameter
`SigPhi`	Standard deviation of the TGSD (Phi)
`Alpha`	Alpha parameter of the Beta function controlling the mass distribution in the plume
`Beta`	Beta parameter of the Beta function controlling the mass distribution in the plume
`WindSpeed`	Wind speed (m/s) - used when fixedWind=0
`WindDir`	Wind direction (degree from N) - used when fixedWind=0

Constant parameters

This section contains parameters that are kept constant during the inversion.

Variable	Description
`lithicDensity`	Density of the lithics (kg/m³)
`pumiceDensity`	Density of the pumices (kg/m³)
`minPhi`	Minimum bound (coarsest diameter; phi) of the TGSD
`maxPhi`	Maximum bound (finest diameter; phi) of the TGSD

Batch parameters

This section contains parameters that control the behavior of the batch inversion. It is still experimental, so it is recommended to keep the proposed parameters for now.

Variable	Description
`deltaMass`	Mass interval between each sub-run of the batch inversion (log10; kg)
`incrMass`	Mass increment between each sub-run of the batch inversion (log10; kg)
`deltaHt`	Height interval between each sub-run of the batch inversion (m asl)
`incrHt`	Height increment between each sub-run of the batch inversion (m asl)

Tephra2 advanced parameters

Variable	Description
`colSteps`	Number of integration steps along the height of the column
`partStep`	Number of integration steps along the TGSD
`fitTest`	Merit test: 0 (chi-squared test), 1 (root mean squared error), 2 (Tokyo log test)
`eddy`	Eddy constant
`plumeModel`	Mass distribution in the plume (2 is a beta distribution)
`windLevels`	Number of wind levels to consider. Only if fixedWind=0

Running the inversion

Tephra2 is parallelised using OpenMPI and commonly run on clusters. Different clusters have different architectures, which will probably require some measure of editing. If your files are prepared locally, now is the time to upload everything on the server. Be sure to keep the directory tree specified before.

The main script used to run the inversion is runInversion.sh. The _template/ folder contains templates for different cluster architectures, namely Slurm, OpenPBS and PBSPro. Edit the header of the script and copy it to the _example_folder/ folder. Navigate to _example_folder/ to start running the script.

Single vs. batch vs. seed runs

For parallelisation purposes, single runs, batch runs and seed runs are submitted in slightly different ways.

For single runs, the mass range is defined in the inversionConfig.conf file;
For batch runs use job arrays, the mass range is specified on the command line when submitting the job, which will allow each mass increment to be sent to a different node;
For seed runs, the mass is specified in inversionConfig.conf (as for single runs), but the number of seeds is specified on the command line when submitting the job.

OpenPBS

For a single inversion run, set BATCH=0 in inversionConfig.conf and type:

qsub runInversion.sh

For a batch inversion run with a mass range between 10⁹ and 10¹¹ kg, set BATCH=1 in inversionConfig.conf and use the -t flag:

qsub -t 9-11 runInversion.sh

For a seed inversion performing 100 runs with a seed number varying between 1-100, type:

qsub -t 1-100 runInversion.sh

Post-processing

Upon successful (yey!) completion of an inversion run, one (single run) or multiple (batch run) folders are created inside _example_folder/, named massM_htH, where M and H are the lower intervals of the mass (log10 kg) and plume height (km) ranges, respectively. For post processing, follow these steps:

In the _example_folder/ on your local computer, create a folder named n/, where n is the number of your inversion attempt (start at 1/ and name the subsequent folders in a continuous way). By doing so, each new inversion attempt doesn’t have to delete the previous one, and it becomes easier to track different attempts
Copy the content of _example_folder/ on the cluster into your local n/ folder
In Matlab, run the processInversion.m script. When asked, go and select your n/ folder and let the script work

This process will first create four new files in each folder:

obsVScomp.pdf: Plot of the computed values of tephra accumulation against their observed equivalents;
plume.pdf: The plume distribution
tephra2.out: Result of the forward solution
wind.pdf: Wind profile

Additionally, each inversion attempt is recorded in an Excel file at the root of the Inversion/ folder in a different sheet named after n/.

Output figures

1. Spaces of solutions

In the case of batch or seed runs, the function plotT2Inversion can be used to explore spaces of results. A plane of fit values is computed and interpolated between two variable of choice. In the case of batch runs, this visualisation method is useful to assess the location(s) of minima. In the case of seed runs, these figures are helpeful to assess the stability of the inversion output to the seed number.

2. Everything against everything!

In the case of seed runs, the script produces a matrix plot (type doc plotmatrix in the Matlab command line). I find this plot to be the most informative because:

It allows indentifying populations of solutions from the histograms
It allows indentifying relationships between parameters

Access sampled ESP in TephraProb

2018-02-21T00:00:00+00:00

This note shows how to access the ESPs sampled during Monte Carlo simulations and how to link them to a given output. This note has been updated on 2020-02-12 to work with v1.7.1 of TephraProb. Two files are needed:

The file located in RUNS/runName/runNb/ named runName_runNb.mat. This contains a structure called data, which contains a field named stor where all the ESPs are stored. data.stor contains sub-fields for the different parameter (e.g. plume height, duration, mass, MER, TGSD). These parameters are stored as vectors that can have different sizes. If the duration of each simulated eruption is <6 h, the size of all vectors should equal the number of simulated eruptions. If durations are >6 h, the size of those parameters that change at each 6 h increment (e.g. plume height, MER) represents the number of times each was sampled. For instance, in the case of a 15 h-long eruption, 3 plume heights are sampled. However, other parameters (e.g. TGSD) are constant over the entire duration of the eruption, so the resulting vectors will have a size equal to the number of simulated eruptions.
The file located in RUNS/runName/runNb/DATA/ named dataT2_all.mat, which contains two variables:
- dataT2 is a m by n matrix, where m is the number of points over which accumulations are computed and n is the number of runs. This variable contains all the summed Tephra2 outputs in kg/m2;
- runNb is a m by 1 vector, where m is the number of runs.

Now, there is something to keep in mind. Let’s say we have 10 text files named sequentially from 01.txt to 10.txt. Matlab reads sequential files as 01.txt, 10.txt, 02.txt … 09.txt, and not as 01.txt, 02.txt … 10.txt. In TephraProb post-processing, this means that:

ESPs in the variable data.stor are ordered in increasing run number
Runs outputs in dataT2 are ordered in the way Matlab reads the files

It is therefore necessary to re-index dataT2 to the order of the ESPs contained in data.stor.

Examples

Re-index ESPs with outputs

% Load the second sampling attempt of the run called Sakurajima
load RUNS/Sakurajima/2/Sakurajima_2 

% Load the isomass output
load RUNS/Sakurajima/2/DATA/dataT2_all.mat 

% Show the loaded variables
data.stor
dataT2
runNb

% Retrieve the selected ESPs
esp = [[data.stor.ht]'./1e3,log10([data.stor.mass]'),[data.stor.dur]'./3600,[data.stor.gs_med]',date(:,2)];

% Reorder the run numbers in ascending order and save their indices
[~,id] = sort(runNb);
% Reorder the isomass outputs
dataT2 = dataT2(:,id);

% dataT2 and data.stor are now indexed
random_run = 10;
run_output = dataT2(:,random_run);
run_ESPs = esp(random_run, :);

Plot a histogram of the plume height:

toPlot = [data.stor.ht]; % Need to vectorize the data

figure;
histogram(toPlot)
xlabel('Plume height')
ylabel('Frequency')
title([num2str(numel(toPlot)), ' occurrences'])

Manual compilation of Tephra2

2017-12-11T00:00:00+00:00

TephraProb compiles Tephra2 through a bit of Matlab hocus pocus to facilitate problems of OS dependancy, but the mighty Matlab can sometime fail us. In this case, it becomes necessary to manually compile Tephra2 from the command line to see what the error message is. I spent a fair amount of time writing this procedure in emails to people who encountered problems and I thought I would write it down here.

Just to clarify things, compiling is the process of transforming files containing lines of code into a program that can be executed by the user. The whole compilation process depends on such aspects as the computer’s architecture or OS, so it needs to be performed on each computer separately. Again, most of the time everything should be taken care of smoothly within TephraProb, but just in case…

Unix for non-nerds

This process requires to access Tephra2 through the command line (Terminal on Unix, Cygwin Terminal on Windows). Here are some tips on how to use it.

Navigating through folders

The main command to navigate through folders is cd:

cd folder 	# Enter the directory called folder
cd .. 		# Go back to the parent directory

Where to go?

Another useful command is ls, which lists the content of the current directory. When navigating directories with cd, note that the tab key auto-completes commands. Taking the example above, typing cd fold and hitting tab is likely to complete the command to cd folder/.

Locating Tephra2

The Tephra2 executable (tephra2-2012 on Unix, tephra2-2012.exe on Windows) should be located in the tephraProb/MODEL/ folder. Open the terminal and navigate to this location. If you are unsure in which directory you are when starting the Terminal, the command pwd will help. On Unix, it usually starts in your home folder (/home/username/ on Linux, /Users/username/ on Mac).

On Windows/Cygwin, things are slightly different. Let’s assume that TephraProb is located in E:/TephraProb/, the right way to navigate to this location on the Cygwin terminal is:

cd /cygdrive/e/TephraProb

Manual compilation

Once in the MODEL/ folder, enter the forward_src folder and type:

make

Hopefully, that should result in the last line printed to the command line being:

mv tephra2-2012 ../tephra2-2012

There might be warnings, but these should be ok.

Common errors

Here, I will try and list common errors that have been encountered.

File format not recognised

This error was encountered on Cygwin:

make[1]: Entering directory '/cygdrive/e/Seb/Documents/CODES/TephraProb/MODEL/forward_src'
`which gcc` -Wall -g -idirafter ../include -o tephra2-2012 new_tephra.o tephra2_calc.o -lm /libcygwin64/libgc.a -ldl
new_tephra.o: file not recognized: File format not recognized

This was solved by entering in the forward_src, removing the *.o files and recompiling:

cd forward_src
rm *.o 	# Remove *.o files
make

Manual Tephra2 run

If the compilation worked, the Tephra2 executable should now be located in the folder MODEL/. To start the command, type ./tephra2-2012 on Unix or ./tephra2-2012.exe on Windows, which will result in:

Missing comand line arguments,
USAGE: <program name> <config file> <points file> <wind file> <file of grain sizes>

You now have to specify all input files including:

A configuration file, typically a .conf file located somewhere in TephraProb/RUNS/runName/runNumber/CONF/
A grid file, typically a .utm file located in TephraProb/GRID/gridName/
A wind file, typically a .gen file located in TephraProb/WIND/windName/
A grain-size distribution file, typically a .gsd file located somewhere in TephraProb/RUNS/runName/runNumber/GS/

You can either specify the full path to these files or copy them in the same folder as the Tephra2 executable, in which case Tephra2 can be run by typing:

./tephra2 conf.conf grid.utm wind.gen tgsd.gsd > test.out # Unix
./tephra2.exe conf.conf grid.utm wind.gen tgsd.gsd > test.out # Windows

Note that the last argument > test.out specifies the output file to Tephra2.

Common errors

Permission denied

If a Permission denied error is returned upon using the ./tephra2-2012 command, try changing the access permissions:

chmod 755 tephra2-2012

Plot the output of Tephra2 in Matlab

2017-12-10T00:00:00+00:00

This function plots the output of Tephra2 runs on a map. Output files generated by Tephra2 are in a column format (i.e. easting, northing, elevation and mass accumulation), with the number of rows equal to the total number of points in the grid. For some obscure reason, most of Matlab’s plotting functions only accept inputs in a matrix format, i.e. separate m x n matrices for easting, northing, elevation and mass accumulation, where m and n are the numbers of northing and easting values, respectively. Matlab, unlike other plotting tools (e.g. GMT, Python), therefore requires a transformation the data.

Additionally, Tephra works with UTM coordinates. This function accepts optional arguments to convert the output of Tephra2 to lat/lon coordinates. For this, two additional functions published on the Matlab File Exchange are required:

utm2ll by François Beauducel
plot_google_map by Zohar Bar-Yehuda

The function

The function plotT2 can be downloaded here. One argument is required, which is a string containing the path to the file. Outputs are then matrices of easting (E), northing (N) and mass accumulation (M):

[E,N,M,Carea,Pmass,cont,varargout] = plotT2(pathToFile, varargin)

Required argument

Argument	Description
`pathToFile`	Path to Tephra2 output file

Optional name-value pairs of arguments

The optional pairs of name-value arguments are summarised below:

Argument	Description
`'zone'`	UTM zone of the Tephra2 output, used to convert projected to geographic coordinates and plot a Google Map background if both utm2ll and plot_google_map are available. Should be a double, positive in the N hemisphere, negative in the S hemisphere
`'plot'`	Defines if accumulation is plotted as log10 ('log10', default) or linear ('linear')
`'contours'`	Vector containing the values (kg/m²) of accumulation to contour. By default set to [1,10,100,1000]</td>
`'minVal'`	Minimum value to be represented on the continuous color surface (kg/m²)
`'vent'`	Plot the vent, entered as [easting, northing] coordinates.
`'points'`	Additional points to plot, entered as a 2-columns [easting,northing] matrix. If entered as a 3-columns [easting,northing,value] matrix, value is used as a label

Output arguments

Argument	Description
`E`	Easting matrix (m)
`N`	Northing matrix (m)
`M`	Mass accumulation matrix (kg/m²)
`Carea`	Area of each isomass contour (km²)
`Pmass`	Mass contained within each contour (kg)
`contour`	Contours (kg/m²)
`lat`	Optional - latitude matrix, if `'zone'` is specified
`lon`	Optional - longitude matrix, if `'zone'` is specified

Examples

When the argument 'zone' is passed, it is possible to request the matrix versions of latitude and longitude:

plotT2('example.out');
plotT2('example.out', 'zone', 18);
plotT2('example.out', 'plot', 'linear');
plotT2('example.out', 'contours', [10, 50, 100]);

Quantifying probabilities from ballistic models

2017-10-09T00:00:00+00:00

Here is presented an extended version of the Matlab post-processing function developed for the Great Balls of Fire model for the probabilistic hazard assessment of ballistic impact.

Table of content

tags: gbf
1. From ballistic modelling to probabilistic hazard assessment
- 1.1 The need of a reference area
- 1.2 What probabilities?
2 Using the function
- 2.1 Input arguments
- 2.2 Output argument
3 Outputs and results
- 3.1 ASCII grid format
- 3.2 Matlab format
  - .bal
  - .pixel
  - .concentric
  - .radial
4 Displaying results

Requirements

The initial function provided with the Great Balls of Fire model has been modified so output files of virtually any model can be loaded and processed through methods described here. Requirements are:

Impact coordinates must be defined in UTM WGS84
Files must be ascii and tab delimited
Columns in the output file should be organised as indicated in the table below. Note that adding extra columns is acceptable

Column	Data	Unit
1	Easting	m
2	Northing	m
3	Altitude	m asl
4	Mass	kg
5	Diameter	m
6	Kinetic energy	kJ

Table 1: Example of an input file for the processGBF function.

Source material

Source material includes:

Our JVGR paper
The GBF user manual
This 2017 IAVCEI presentation

Dependencies

The post-processing functions rely on two files published on the Matlab File Exchange:

utm2ll by François Beauducel
plot_google_map by Zohar Bar-Yehuda

Access the code

The post-processing function can be obtained here

1. From ballistic modelling to probabilistic hazard assessment

1.1 The need of a reference area

Unlike for tephra deposits, the discrete nature of ballistic impacts requires to average single ballistic impacts on a grid. Grids can fall in two main categories that are polar and cartesian. However, averaging impacts on a grid creates a dependancy of the final result to the spatial resolution. In other words:

Is it acceptable for a grid to have pixels areas that vary with distance from the vent (e.g., polar grids)?
Since increasing the pixel area increases the probability of impact, is there a threshold of grid resolution that provides stable results?
How does the grid shape/resolution interacts with the number of simulated ballistics and affect the stability of results?

Probabilistic hazard assessment for ballistic impacts should therefore address these issues, one possible way to do so being through sensitivity analysis. This new post-processing function can be accessed either through a GUI or programatically, which facilitates sensitivity analyses. Refer to our JVGR paper for more information.

Fig. 1 Illustration of the different approaches to average ballistics including i) polar grids with constant concentric distance from the vent, ii) polar grids with increasing distance from the vent, iii) cartesian grids and iv) area of interests

In this post-processing function, a true polar grid is not implemented. Instead, probabilities of impact are calculated using both a given distance from the vent (i.e., concentric rings) and at a given radial sector around the vent. The union probability of both approaches will result in a polar grid.

1.2 What probabilities?

In frequentist analyses, probabilities must represent a ratio between a number of ballistic impacts satisfying a given characteristics (e.g., a kinetic energy at impact that exceeds a threshold) and a total number of ballistics. Depending on our definition of the second term, let’s call it B_TOT, probabilities are going to express a different message.

If B_TOT is the total number of simulated ballistics, the probability expresses P(Z, E_T), that is:

What is the probability of a zone Z to be impacted by a ballistic impact with an energy threshold greater than E_T?

If B_TOT is the total number of ballistics impacting a given zone, the probability expresses P(E_T|Z), that is:

What is the probability of a ballistic impacting zone Z to exceed an impact energy greater than E_T?

These probabilities are significantly different, and both are considered in the post-processing approach presented here. For simplicity, we will refer to B_TOT as absolute probabilities and B_TOT as relative probabilities.

2 Using the function

2.1 Input arguments

Table 2 summarises the input parameters required by the processGBF function. Functions can be called using a GUI, arguments or a structure. The following steps assume that you are in the root of the folder containing the functions.

Variable	Description	Unit	Type
pth	Path to the output file	-	String
name	Run name	-	String
gridRes	Resolution of the cartesian grid, which can be enter either as a single value or a 1x2 vector. In the case of a 1x2 vector, the first value is the resolution of the grid used for the pixel approach and the second value is the resolution of the grid used to map the result of the concentric and radial approaches. If only one value is entered, the resolution of the grid for the concentric and radial approaches is arbitrarily set to be 10% finer than the grid used for the pixel approach. This is only done to ensure smooth display on maps.	m	Numeric
subset	Subset of the total population of VBPs. Default is 100%, but can be decreased for sensitivity analyses	%	Numeric
eT	Energy threshold(s) used to calculate exceedance probabilities. Multiple values can be defined as a vector	J	Numeric
pT	Probability threshold(s) used to calculate energies occurring at a given exeedance probability. Multiple values can be defined as a vector	%	Numeric
dI	Distance interval used to define concentring rings	m	Numeric
rI	Angle interval used to define radial sectors	Degrees	Numeric
vE	Vent easting	m	Numeric
vN	Vent northing	m	Numeric
vZ	Vent zone defined as a numeric value only (i.e., do not enter the letter) and positive in the N hemisphere and negative in the S hemisphere (numeric)	-	Numeric

Table 2: Input variables required by the processGBF function.

2.1.1 Using the GUI

In the Matlab command, type:

processGBF

This approach will open two consecutive GUIs that will prompt you to i) locate the input file and ii) define the additional variables described above.

2.1.2 Using arguments

This approach allows the user to enter all input argument separately from the command line:

processGBF(pth, name, gridRes, subset, eT, pT, dI, rI, vE, vN, vZ)

All variables refer to Table 1.

2.1.3 Using an input structure

The structure approach combines the arguments 2-9 of the argument approach into a structure. It is useful to use to output of a GBF post-processing run.

processGBF(pth, run)

where run is a Matlab structure named as the arguments 2-9 of the argument approach.

2.2 Output argument

Calling:

project = processGBF(...)

returns a Matlab structure containing all the output data… but stick around until the next section!

3 Outputs and results

Upon completion of the processGBF function, a folder named after the name variable is created, which contains results in a Matlab format (i.e., name.mat) and in a ESRI ASCII grid format. All probabilities are expressed in % and all energies are expressed in Joules.

3.1 ASCII grid format

Three main categories of ASCII files are written:

pixel, which use a cartesian grid
concentric, which use concentric zones around the vent
radial, which use radial sectors around the vent

Each of these categories have sub-categories defined by:

prob, which contain the exceedence probabilities (%) for ballistic impacts to exceed a given threshold of kinetic energy
en, which contain the typical kinetic energy (J) at impact for a given exceedance probabillity

The probability files are themselves differentiated between:

abs, which contain absolute probabilities and
rel, which contain relative probabilities

3.2 Matlab format

The Matlab output file, i.e. name.mat, contains the most comprehensive output. It is compiled as a Matlab structure containing:

inBal, which is a duplicate of the input arguments shown in Table 1
bal, which contains the data of all simulated particles
pixel, which contains the result of the processing performed using a cartesian grid
concentric, which contains the result of the processing performed using concentric zones around the vent
radial, which contains the result of the processing performed using radial sectors around the vent

.bal

Field	Description
x	Easting of each ballistic (n x 1 vector, where n is the total number of simulated ballistics)
y	Northing of each ballistic (n x 1 vector)
lat	Latitude of each ballistic (n x 1 vector)
lon	Longitude of each ballistic (n x 1 vector)
e	Kinetic energy (J) at impact of each ballistic (n x 1 vector)
n	Total number of simulated ballistics
d	Distance from the vent of each ballistic (n x 1 vector)
r	Angle from the vent of each ballistic (n x 1 vector)
xi	X coordinate of each ballistic on the cartesian grid (n x 1 vector)
yi	Y coordinate of each ballistic on the cartesian grid (n x 1 vector)
data	Original data contained in the user-defined input file

.pixel

Field	Description
east	Easting coordinates of the cartesian grid (i x j matrix, where i and j represent the northing and easting intervals of a grid of resolution gridRes covering the entire ballistic field, respectively)
north	Northing coordinates of the cartesian grid (i x j matrix)
lat	Latitude of the cartesian grid (i x j matrix)
lon	Longitude of the cartesian grid (i x j matrix)
Nt	Number of ballistics in each pixel of the cartesian grid (i x j matrix)
E	Kinetic energy (J) at impact for a given exceedance probability computed on the cartesian grid (i x j x pT matrix)
N	Number of ballistics with kinetic energies exceeding a given energy threshold in each pixel of the cartesian grid (i x j x eT matrix)
Pabs	Absolute probability (%) for ballistic impacts to exceed an energy threshold (i x j x eT matrix)
Prel	Relative probability (%) for ballistic impacts to exceed an energy threshold

.concentric

Field	Description
bin	Distances of the concentric rings around the vent (m)
Nt	Number of ballistics in each concentric ring (length(bin) x 1 vector)
N	Number of ballistics with kinetic energies exceeding a given energy threshold in each concentric ring (length(bin) x eT matrix)
E	Kinetic energy (J) at impact for a given exceedance probability computed in each concentric ring (length(bin) x pT matrix)
Pabs	Absolute probability (%) for ballistic impacts to exceed an energy threshold (length(bin) x eT matrix)
Prel	Relative probability (%) for ballistic impacts to exceed an energy threshold (length(bin) x eT matrix)
east	Easting coordinates of the cartesian grid used to represent the concentric rings on a map (i x j matrix). i and j represent the northing and easting intervals of a grid of resolution gridRes(2) (or 10% of gridRes(1) covering the entire ballistic field, respectively
north	Northing coordinates of the cartesian grid used to represent the concentric rings on a map (i x j matrix)
lat	Latitude of the cartesian grid used to represent the concentric rings on a map (i x j matrix)
lon	Longitude of the cartesian grid used to represent the concentric rings on a map (i x j matrix)

Note that binM, NtM, NM, EM, PabsM and PrelM contain the same data as bin, Nt, N, E, Pabs and Prel gridded on a matrix, which can be used along east/north or lon/lat for plotting on a map.

.radial

Field	Description
bin	Distances of the radial sectors around the vent (m)
Nt	Number of ballistics in each radial sector (length(bin) x 1 vector)
N	Number of ballistics with kinetic energies exceeding a given energy threshold in each radial sector (length(bin) x eT matrix)
E	Kinetic energy (J) at impact for a given exceedance probability computed in each radial sector (length(bin) x pT matrix)
Pabs	Absolute probability (%) for ballistic impacts to exceed an energy threshold (length(bin) x eT matrix)
Prel	Relative probability (%) for ballistic impacts to exceed an energy threshold (length(bin) x eT matrix)
east	Easting coordinates of the cartesian grid used to represent the radial sectors on a map (i x j matrix). i and j represent the northing and easting intervals of a grid of resolution gridRes(2) (or 10% of gridRes(1) covering the entire ballistic field, respectively
north	Northing coordinates of the cartesian grid used to represent the radial sectors on a map (i x j matrix)
lat	Latitude of the cartesian grid used to represent the radial sectors on a map (i x j matrix)
lon	Longitude of the cartesian grid used to represent the radial sectors on a map (i x j matrix)

Here again, binM, NtM, NM, EM, PabsM and PrelM contain the same data as bin, Nt, N, E, Pabs and Prel gridded on a matrix, which can be used along east/north or lon/lat for plotting on a map.

4 Displaying results

The displayGBF function offers an interactive GUI to plot the results of the processGBF function. Start it by typing

displayGBF

in the Matlab command. Then select the name*.mat file output of the processGBF function and you are good to go!

Manually importing ERA-Interim files in TephraProb

2017-09-15T00:00:00+00:00

When accessing ERA-Interim data in batch, the behaviour of the ECMWF Python API is often unpredictable, which sometime complicates the automatic downloading of wind data in TephraProb. If the API does not work for you, the option of importing NetCDF files that were manually downloaded from the ERA-Interim website was added to TephraProb. Ain’t that nice?

Note that you still will have to:

Create an account on the ECMWF website
Login to the ECMWF website
Retrieve an API key
Accept the license

1 Download the data

Select variables

Go to this page. On the top left corner of the page, be sure the type of level is pressure level
In the Select date box, select the date interval. You can either manually enter the dates or click in the grid. In any case, two things to keep in mind:
- Time ranges must be continuous
- Requesting longer time periods result in longer request times on their servers. With TephraProb, you can split your time range in smaller subsets.
In the Select time box, select the four time intervals.
In the Select level and parameter box, click on (which should select all levels):
- Geopotential
- U component of wind
- V component of wind
At the bottom of the page, click the Retrieve NetCDF button

Customise the grid

The next page appears and a few customisations need to be done:

Change the grid resolution to 0.25x0.25.
Enter a custom area. TephraProb interpolates the wind above the vent, so it is a good practice to define a grid that spans a few (typically ≥ 2) points around the vent, knowing the grid resolution defined above.
Click the Retrieve now button and the request will start.

Retrieving your request

Right, so depending on the size of your request, the processing on the ECMWF servers might take a while (a whole year usually takes me between 15-30 minutes). If you quit your browser’s page, no worries! This page provides you with an history of your requests.

2 Processing

Once you’ve got all your data ready, switch back to TephraProb and use the Input > Wind > Download wind data tool. Select the ECMWF ERA-Interim Offline Reanalysis dataset and fill the required inputs and here you go!

The effect of particle shape on the estimation of plume height

2017-09-15T00:00:00+00:00

The method of Carey and Sparks (1986) allows to estimate the height of strong plumes based on the geometry of isopleth maps. Although the method has its limitations, it is still the most used in the literature and remains important for the characterisation of eruption source parameters for eruptions that have no direct observations.

Commonly, the contouring of the lithic isopleths is preferred over the juvenile fraction because lithics are less likely to break and because the drag of smooth particles is easier to estimate. Computations behind the method of Carey and Sparks (1986) are somehow obscure and difficult to reproduce and is based on a fixed value of drag coefficient of 1.054 and the assumption of spherical particles.

The method relies on a manual plotting on downwind vs. crosswind ranges for four values of diameters/densities. For the needs of the TError package, we implemented the method of Carey and Sparks (1986) as a Matlab function in order to use it routinely and to be able to work with a continuum of particle sizes and densities (see the reference page for more details). But then, during a long flight back from Japan with my buddy Mohssen Bagheri in which the best movie Lufthansa could propose was the biopic on Jimi Hendrix (sacrilege!), we started wondering how much would the shape parameter of the particle influence the estimate of the plume height.

Method

Ok, remember that the two assumptions of the model of Carey and Sparks we would like to test are:

The use of spherical particles
The use of a fixed drag coefficient of 1.054 (note that even our method does not account for a variable drag coefficient).

Here, I won’t re-write Bagheri’s papers - which you can find at the bottom of the page - but just give an idea of the method taken.

For a given diameter and density, the corresponding plume height and wind speed are calculated with the method of Carey and Sparks (1986)
For the same diameter and density, we calculate the terminal fall velocity for a range of elongation and a flatness values as described in Bagheri and Bonadonna (2016). Note that the elongation is the ration between the intermediate and the longest axes (I/L) whereas the flatness is the ratio between the shortest and the intermediate axes (S/I) (Bagheri et al., 2015).
Equation 3 of Carey and Sparks is inverted to convert this terminal fall velocity to an equivalent diameter.
This equivalent diameter is in turn used to calculate the plume height and the wind speed through the method of Carey and Sparks (1986).

Voilà!

Results

This process is implemented in a Matlab function named CareySparks86. Figure 1 shows the output of the process for elongation and flatness values ranging from 0.1 to 1.0. The square box shows the range of shapes typically relevant for volcanic lapilli (Bagheri et al., 2015) and the red curve shows a fixed value of drag coefficient of 1.054.

Fig. 1 Illustration of the effect of particle shape on the estimation of the plume height (left) and wind speed (right). The red lines contour a drag coefficient of 1.054 specified in Carey and Sparks (1986) for a given clast diameter and density.

Recommendations

The plume height is a critical source parameter, and although recent method help estimating it more precisely when observed, older deposit still strongly depends on the method of Carey and Sparks (1986) and the characterisation of the maximum clasts. The work of Mohssen Bagheri really reduces the uncertainty on the estimation of the drag coefficient and here are the recommendations that come out of his papers for estimating an equivalent diameter:

For a given clast, look for the projection with the maximum area.
The largest axis from this projection defines L and the smallest axis defines I. Note that they do not need to be perpendicular.
Look for the projection with the smallest area. The smallest axis on this projection defines S.
Use the geometric mean of L, I and S as your equivalent diameter.

Manually importing NOAA Reanalysis files in TephraProb

2017-09-12T00:00:00+00:00

It is sometime necessary to import into TephraProb some NOAA Reanalysis NetCDF files downloaded manually. One of the reasons is that the NOAA server is accessed via FTP, and firewalls (such as those installed in University networks) might result in an error such as this:

Firewalls - such as those installed in University networks- can block the FTP access to NOAA servers and result in this kind of error:

Error using ftp/mget (line 109)
FTP error: 425.

In this case, it is necessary to manually download the NOAA Reanalysis NetCDF files. When working with the NOAA Reanalysis dataset (either the NCEP/NCAR Reanalysis 1 or NCEP-DOE Reanalysis 2 versions), TephraProb will look for NetCDF files located in the folder WIND/_Reanalysis1_Rawdata/ or WIND/_Reanalysis2_Rawdata/. Each file is named variable/year.nc, where variable can be hgt (i.e. geopotential height), uwnd (i.e. u wind component) and vwnd (i.e. v wind component). Each of this file contains data for an entire year for the Global Grid defined by NOAA.

When using the function Input > Wind > Download wind data from the TephraProb GUI, the code will first check whether files for the requested time period already exist in the _Rawdata/ folder in order to avoid multiple downloads. If they do not, well, then they are downloaded, else the downloading part is skipped and files are automatically processed, where temporal and spatial subsets are extracted.

What we will do is manually download files from the NOAA website from a usual web browser in order to bypass the use of FTP and manually place them in the _Rawdata/ to cheat the downloading part.

1 Download the data

Download the raw NetCDF files from the following links:

Reanalysis 1: hgt uwnd vwnd
Reanalysis 2: hgt uwnd vwnd

2 Processing

Move the downloaded NetCDF files to the folders WIND/_Reanalysis1_Rawdata/ or WIND/_Reanalysis2_Rawdata/. Now use the Input > Wind > Download wind data from the TephraProb GUI as you would normally do, specify the right dataset and voila!

Seb's repository

TephraProb v.1.7

Changes in TephraProb v1.7

Added functionalities

ERA-5

Relevant ERA-5 documentation

Loading wind datasets

Logarthmic sampling of mass

Linking outputs to original ESPs

Bug fixes

TephraProb troubleshoot

Content

ECMWF-Related issues

Pip

Case ney :(

Case yey!

SSL

Setuptools

ECMWF Python API

Tephra2 related issues

Mac OS 10.14 Mojave

Random tools

Github repositories

Python Resources

Tutorials

Github repositories

GEE Resources

Example scripts

Toolbox

Tutorials

Papers

Tips

Maps

Find wind profiles most similar to an average

Access the tool

Method

Use the tool

TephraProb v.1.6

Changes in TephraProb v1.6

TephraProb preferences

Google Earth output file

Added dependency

Wind data for Tephra2

Null elevations

End-of-the-file

Conclusion

How to notice problems?

TephraProb v.1.5

Changes in TephraProb v1.5

Probability calculations

Export to ASCII files

Archive Tephra2 files

Probabilistic isomass maps

Hazard curves

Map preferences

Bug fix

TephraProb on a cluster

Setting up remote files

Compiling Tephra2

Running in parallel

SLURM

OpenPBS

Post-processing

Tephra2 Inversion

Updates - November 2018

Table of content

Introversion

Single vs. batch runs

Caveats

Installation

Getting started

Case-study files

Configuration file

Running the inversion

Single vs. batch vs. seed runs

OpenPBS

Post-processing

Output figures

1. Spaces of solutions

2. Everything against everything!