### Modeling and animating the evolution of a cinder cone in Python

Numerically modeling of the evolution of a cinder cone where individual rock particles or "bombs" are launched from a central point with "random" trajectories from a central location and deposition/accumulation. Each process is explicit, eventually this process develops the shape of a classic cinder cone.

I wrote this as part of a graduate course and added comments throughout so it is easy to follow what is happening. Also included is an example of how to save images of the process and combine them into a gif animation.

Scientific and numerical modeling in Python is best accomplished with the Numpy library. This book (affiliate link) is a great resource if you are interested in improving your numerical methods in Python using Numpy, SciPy, and matplotlib.

%matplotlib inline import numpy as np import matplotlib.pyplot as plt import imageio from pathlib import Path import re

Define physical parameters of bomb size and trajectories, determined from literature and trial and error.

## Bomb velocities, m/s vel_min = 32 vel_max = 76 ## trajectory angles, degrees ang_min = 345 ang_max = 375 ## diameter range, in meters dia_min = .1 dia_max = 0.9

Functions for creating, launching, and depositing bombs and defining the spatial domain.

###### ## bomb creation ###### def makeBomb(vMin,vMax,aMin,aMax,dMin,dMax,num=1): """Inputs are velocity (m/s), angle (degrees), and diameter (m) min and max values outputs are particle resultant velocity, angle(rad), and area (sq m)""" vel = (vMax - vMin)*np.random.random(num) + vMin ## velocity in m/s ang = (aMax - aMin)*np.random.random(num) + aMin ang *= np.pi/180 ## return angle in radians dia = (dMax - dMin)*np.random.random(num) + dMin are = np.pi*(dia/2.)**2 ## return area in square meters return {'velocity':vel,'angle':ang,'area':are} ## returned as dict for readability ## x spatial grid for simulation and plotting x_grid = np.linspace(-1000,1000,2001) # unit meter grid ###### ## bomb launching ###### def launchBomb(vel,ang,z_0=0,x=x_grid): """ Inputs: velocity (vel), angle of launch in rads (ang), launch elev (z_0), and x grid that is defaulted to use a predifined grid in main program return x and z position vectors """ g = 9.81 ## m/s^2 x_vel = vel * np.sin(ang) z_vel = np.abs( vel * np.cos(ang)) z = np.zeros_like(x) z = (z_vel * x)/x_vel - (1/2.) * g * (x/x_vel)**2 + z_0 return z ## elevation position array ###### ## bomb depsiting ###### def depositBomb(x_traj,z_traj,elev_profile,bombArea): """Inputs: x distance grid, z (elevation) coordinates for a trajectory, previous land elevation profile and bomb area. Finds where bomb hits ground and adds accumulated mass to the elevation profile and returns it""" vent_ind = np.where(x_traj == 0)[0][0] ## get index of vent (middle of x grid) ## check if lauch went to the right if np.where(z_traj == np.max(z_traj)) > vent_ind: x = np.where(z_traj[vent_ind:] <= elev_profile[vent_ind:])[0] ## where trajectory below elevation deposit_ind = np.min(x[1:]) + vent_ind -1 ## get the minimum index of occurence (first intersection with land) ## redeposit on cliffs while elev_profile[deposit_ind] - elev_profile[deposit_ind + 1] > np.sin(np.pi/6): ## forward roll deposit_ind += 1 while elev_profile[deposit_ind] - elev_profile[deposit_ind - 1] > np.sin(np.pi/6): ## backward fall deposit_ind -= 1 ## prevent indefinate rolling if steep in both directions if elev_profile[deposit_ind] - elev_profile[deposit_ind+1] > bombArea: break ## check if launch went to the left of vent elif np.where(z_traj == np.max(z_traj)) < vent_ind: deposit_ind = np.max(np.where(z_traj[:vent_ind] <= elev_profile[:vent_ind])[0]) ## ##redeposit particles on cliffs while elev_profile[deposit_ind] - elev_profile[deposit_ind - 1] > np.sin(np.pi/6): # forward deposit_ind -= 1 while elev_profile[deposit_ind] - elev_profile[deposit_ind + 1] > np.sin(np.pi/6): # backwards deposit_ind += 1 ## Prevent infinite movement if steep on both directions if elev_profile[deposit_ind] - elev_profile[deposit_ind-1] > bombArea: break ## check if bomb landed in the vent itself else: deposit_ind = vent_ind ## need to fix redistribution in the vent, prevent huge buildup if elev[deposit_ind] > elev[deposit_ind+1] or elev[deposit_ind] > elev[deposit_ind-1]: deposit_ind += np.random.choice(np.arange(-15,15,1)) ## create some chaos to redistribute vent deps elev_profile[vent_ind-5:vent_ind+5] -= 0.1 ## cause meltdown of vent occasionally to prevent filling in ## increment elevation profile at the correct index by area of bomb elev_profile[deposit_ind] += bombArea return elev_profile

Combining semi-random bomb creation, ejection, and deposition in a loop. This is also an example of how you can plot and save a series of images and then combine them to make an animation of the physical simulation.

###### ## Launch Zone ###### ## initialize elevation profile using zero grid as input bomb = makeBomb(vel_min,vel_max,ang_min,ang_max,dia_min,dia_max) traj = launchBomb(bomb['velocity'],bomb['angle']) elev = depositBomb(x_grid,traj,np.zeros_like(x_grid),bomb['area']) ## number of eruptions and how often to show ndiv = 5000 nerupts = 400000 ## for final figure, not movie fig = plt.figure(figsize=(12,9)) ax = fig.add_subplot(111) ## loop through number of eruptions simulating cinder cone development for i in range(nerupts): bomb = makeBomb(vel_min,vel_max,ang_min,ang_max,dia_min,dia_max) traj = launchBomb(bomb['velocity'],bomb['angle'],z_0=elev[np.where(x_grid == 0)[0]]) elev = depositBomb(x_grid,traj,elev,bomb['area']) if (i % ndiv == 0): plt.plot(x_grid,elev) ## plot updated elevation profile plt.xlim(-200,200) plt.show() plt.pause(0.2) ## will plot current trajectory plt.plot(x_grid[traj>elev],traj[traj>elev], 'r.', ms = 1.5) plt.savefig(f'fig{i}.png') # for animation ## plot final elevation profile with every ndiv profile and trajectory overlain plt.plot(x_grid,elev) plt.xlim(-400,400) plt.ylim(0,650) plt.ylabel('Height (m)') plt.title('Cinder cone evolution after %i eruptions,\n deposits and trajectory shown every %i th eruption' % (nerupts,ndiv), fontsize = 16) plt.savefig('Cinder cone evolution after {} eruptions.png'.format(nerupts), dpi = 600)

Gather the saved images and use a natural (human) sorting to order them from start to finish then save as a gif animation with the imageio Python library.

def natural_key(string_): """See https://blog.codinghorror.com/sorting-for-humans-natural-sort-order/""" return [int(s) if s.isdigit() else s for s in re.split(r'(\d+)', string_)] images = sorted(list(map(str, Path().glob('*.png'))), key=natural_key) images = [x for x in images if not 'Cinder' in x] images = [imageio.imread(i) for i in images] imageio.mimsave('eruption.gif',images,fps=3)

And now for the fun part, click on the video below to watch the resulting animation of cinder cone evolution and bomb trajectories plotted every 5000 eruptions:

## Comments

## Post a Comment