Ensemble simulations

Let's simulate a volcanic eruption at Fuego using GEFS forecasts. The following input files are required:

  • fuego.inp: FALL3D configuration file
  • fuego.gfs_pXX.nc: A list of GEFS files in netCDF format
  • GFS.tbl: Dictionary file with GFS variable names (included in the FALL3D distribution under the folder Other/Meteo/Tables)
  • JOB.cmd: A job script with a series of directives to inform the batch system about the characteristics of the job
  • Fall3d.r8.x: FALL3D executable program

Copy the required files to your projects folder:

cd /gpfs/projects/nct01/$USER
cp -r /gpfs/projects/nct00/nct00014/FALL3D_material/hands-on-3 .

Enter to the new folder and create a symlink to the FALL3D executable file. In other words, if the FALL3D installation path is $FALL3D_PATH, run the commands:

cd hands-on-3
ln -s $FALL3D_PATH/bin/Fall3d.r8.x .

Now you can run FALL3D!

Introduction and motivation

  • Atmospheric dispersion models can provide realistic distributions of airborne volcanic ash and gases or tephra deposits
  • Traditionally, operational forecast systems rely on volcanic ash transport and dispersal (VATD) models to produce deterministic forecasts

Why ensemble modelling?

  • Uncertainty in model input parameters: Deterministic models are highly sensitive to uncertain model input parameters (e.g. eruption source parameters) and meteorological fields. We can take into account these uncertainties using ensemble modelling
  • Quantification of model output uncertainty: Ensemble-based modelling allows one to characterise and quantify model output uncertainties. In addition to traditional forecasting products, the associated errors can be provided
  • Improvement of forecast skill: Real observations can be incorporated into dispersal models using ensemble-based data assimilation techniques
  • Source inversion: Different techniques for source term inversion have been proposed based on ensemble modelling

Running ensemble runs

Ensemble simulations can be performed as a single parallel task. In order to perform ensemble runs, FALL3D must be executed with the optional argument -nens to define the ensemble size. For example, the following command will generate a 12-member ensemble and perform the FALL3D task for each ensemble member:

mpirun -np 12 ./Fall3d.r8.x FALL3D name.inp -nens 12

A new folder structure will be created and the results for each ensemble member will be organized in different sub-folders.

Submitting jobs

The job script JOB.cmd contains a series of directives to inform the batch system about the characteristics of the job:

#!/bin/bash
#SBATCH --job-name=FALL3D
#SBATCH --output=%x_%j.out
#SBATCH --error=%x_%j.err
#SBATCH --nodes=1
#SBATCH --ntasks=48
#SBATCH --time=00:15:00
#SBATCH --qos=training
#SBATCH --reservation=Computational24

module purge
module load intel/2017.4
module load impi/2017.4
module load netcdf

INPFILE="fuego.inp"
FALLTASK="all"
NX=2
NY=2
NZ=1
NENS=12

if [ "${NENS}" -gt 1 ] ; then
    for i in $(seq ${NENS})
    do
        ENSDIR="$(printf "%04d" ${i})"
        IENS="$(printf "%02d" ${i})"
        echo "Creating folder ${ENSDIR}"
        mkdir -p ${ENSDIR}
        ln -sfr fuego.gfs_p${IENS}.nc ${ENSDIR}/fuego.gfs.nc
    done
fi

srun ./Fall3d.r8.x ${FALLTASK} ${INPFILE} ${NX} ${NY} ${NZ} -nens ${NENS}

You can submit the job with sbatch:

sbatch JOB.cmd

In this case, we are requesting 48 tasks to run 12 (ensemble size) instances of FALL3D.

Visualizing and analyzing model outputs

Once the model has run, the task PosEns can be executed to merge and post-process the outputs from individual ensemble members (fuego.res.nc) in order to produce a single netCDF file containing ensemble-based deterministic and/or probabilistic outputs for all variables of interest (e.g. concentration at native model levels or at flight levels, cloud column mass, ground deposit load, etc...). Run the command:

mpirun -np 12 ./Fall3d.r8.x PosEns fuego.inp -nens 12

to generate a single ensemble output file: fuego.ens.nc. The content of this file depends on the ENSEMBLE_POST block definition in the configuration file:

 --------------------
 ENSEMBLE_POSTPROCESS
 --------------------
   !
   POSTPROCESS_MEMBERS      = yes
   POSTPROCESS_MEAN         = yes
   POSTPROCESS_LOGMEAN      = no
   POSTPROCESS_MEDIAN       = no
   POSTPROCESS_STANDARD_DEV = no
   POSTPROCESS_PROBABILITY  = no
   POSTPROCESS_PERCENTILES  = no
   !
   IF_POSTPROCESS_PROBABILITY 
      CONCENTRATION_THRESHOLDS_(MG/M3) = 2
      COLUMN_MASS_THRESHOLDS_(G/M2)    = 1
      COLUMN_MASS_THRESHOLDS_(DU)      = 100
      GROUND_LOAD_THRESHOLDS_(KG/M2)   = 1 
   !
   IF_POSTPROCESS_PERCENTILES
      PERCENTILE_VALUES_(%) = 50 

where you can enable/disable different deterministic and probabilistic outputs.

Visualization using Python

The model results can be plotted using the python package Cartopy. First, activate the following Anaconda environment:

module load anaconda
source activate volcanology

and run the following python script under the folder POSTPROCESSING:

python plot_map.py

to plot the ensemble mean of SO2 column mass for every time.

Content of file plot_map.py. Click to expand!
import numpy as np
import xarray as xr
import matplotlib
import matplotlib.pyplot as plt
from matplotlib.colors import BoundaryNorm
import cartopy.crs as crs
import cartopy.feature as cfeature

###
### Parameters
###
minval          = 10
key             = "SO2_col_mass_mean"
fname           = "../fuego.ens.nc"
levels          = np.arange(0.0,3000,200)
vlon, vlat      = -90.88, 14.473
cmap            = plt.cm.RdYlBu_r

###
### Set mininmum level
###
if minval>0: levels[0] = minval

###
### Read file
###
ds = xr.open_dataset(fname)

###
### Generate map
###
proj = crs.PlateCarree()
fig, ax = plt.subplots( subplot_kw={'projection': proj} )

###
### Add map features
###
BORDERS = cfeature.NaturalEarthFeature(
        scale     = '10m',
        category  = 'cultural',
        name      = 'admin_0_countries',
        edgecolor = 'gray',
        facecolor = 'none'
        )
LAND = cfeature.NaturalEarthFeature(
        'physical', 'land', '10m',
        edgecolor = 'none',
        facecolor = 'lightgrey',
        alpha     = 0.8
        )

ax.add_feature(LAND,zorder=0)
ax.add_feature(BORDERS, linewidth=0.4)

###
### Add grid lines
###
gl = ax.gridlines(
	crs         = crs.PlateCarree(),
    draw_labels = True,
	linewidth   = 0.5,
	color       = 'gray',
	alpha       = 0.5,
	linestyle   = '--')
gl.top_labels    = False
gl.right_labels  = False
gl.ylabel_style  = {'rotation': 90}

###
### Add vent location
###
ax.plot(vlon,vlat,color='red',marker='^')

###
### Plot contours
###
ax.set_title("Ensemble mean", loc='left')
cbar = None
for it in range(ds.time.size):
	time_fmt = ds.isel(time=it)['time'].dt.strftime("%d/%m/%Y %H:%M").item()
	ax.set_title(time_fmt, loc='right')
	fc = ax.contourf(
		ds.lon,ds.lat,ds.isel(time=it)[key],
		levels    = levels,
		norm      = BoundaryNorm(levels,cmap.N),
		cmap      = cmap,
		extend    = 'max',
		transform = crs.PlateCarree()
		)

	###
	### Generate colorbar
	###
	if not cbar:
		cbar=fig.colorbar(fc,
			orientation = 'horizontal',
			label       = 'SO2 column mass [DU]',
			)

	###
	### Output plot
	###
	fname_plt = f"map_{it:03d}.png"
	plt.savefig(fname_plt,dpi=300,bbox_inches='tight')

	###
	### Clear contours
	###
	for item in fc.collections: item.remove()

You can create an animated gif using the following command:

convert -delay 10 -loop 0 *.png animation.gif