This tutorial runs through the implementation of a method of estimating TSF runout volumes as proposed by Fontaine & Martin (2015) presented at the Tailings & Mine Waste Conference in Vancouver in 2015.
The basic methodology presented in the paper will be turned into some Python code, and then we'll show how sensitivity analysis can be run and we'll recreate the graphs presented in the paper using Muk3D's scripting tools for plotting.
This example requires Muk3D v2019.4.1 or higher to run. 
Developing the basic script
Defining tailings deposit characteristics
First off, we'll defined the basic properties of the tailings deposit and calculate some basic characteristics.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 
def check(label, calculated_value, paper_val):
print '{}  difference={:.3f}%'.format(label, (calculated_value  paper_val)/paper_val)
mass_tailings_stored = 200000000 # t
average_dry_density = 1.4 # t/m3
specific_gravity = 2.7
tailings_saturation = 1.0 # % as dec fraction, so 1.0 == 100%
density_water = 1.0 # t/m3
volume_tailings_deposit = mass_tailings_stored / average_dry_density # eq. 2
check( 'Volume tailings deposit', volume_tailings_deposit, 143000000)
volume_tailings_solids = mass_tailings_stored / specific_gravity # eq. 3
check('Volume tailings solids', volume_tailings_solids, 74000000)
volume_voids = volume_tailings_deposit  volume_tailings_solids # eq. 4
check('Volume voids', volume_voids, 69000000)
porosity = volume_voids / volume_tailings_deposit # eq. 5
check('Porosity', porosity, 0.48)
void_ratio = volume_voids / volume_tailings_solids # eq. 6
check('Void ratio', void_ratio, 0.93)
volume_interstitial_water = volume_voids * tailings_saturation # eq. 7
check('Volume interstitial water', volume_interstitial_water, 69000000)
mass_interstitial_water = volume_interstitial_water * density_water # eq. 8
check('Mass interstitial water', mass_interstitial_water, 69000000)
moisture_content = mass_interstitial_water / mass_tailings_stored # eq. 9
check('Moisture content', moisture_content, 0.34)
density_bulk = (mass_tailings_stored + mass_interstitial_water) / volume_tailings_deposit # eq. 10
check('Bulk density', density_bulk, 1.9)

While we're developing this script, its a good idea to add some checks in. In this case, in Lines 1 & 2, a small function is used to print a message (characteristic name) and calculate the % error difference between the calculated value and the calculated value in the paper. There will obviously be rounding errors, but the difference should still be very small. This technique can help you identify if you've mistyped a calculation since you'll suddenly see a large error difference.
The various characteristics defined by equations 2  10 are calculated in the above script.
Defining the supernatent pond volume
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 
def check(label, calculated_value, paper_val):
print '{}  difference={:.3f}%'.format(label, (calculated_value  paper_val)/paper_val)
mass_tailings_stored = 200000000 # t
average_dry_density = 1.4 # t/m3
specific_gravity = 2.7
tailings_saturation = 1.0 # % as dec fraction, so 1.0 == 100%
density_water = 1.0 # t/m3
volume_supernatent = 10000000 # m3
volume_tailings_deposit = mass_tailings_stored / average_dry_density # eq. 2
check( 'Volume tailings deposit', volume_tailings_deposit, 143000000)
volume_tailings_solids = mass_tailings_stored / specific_gravity # eq. 3
check('Volume tailings solids', volume_tailings_solids, 74000000)
volume_voids = volume_tailings_deposit  volume_tailings_solids # eq. 4
check('Volume voids', volume_voids, 69000000)
porosity = volume_voids / volume_tailings_deposit # eq. 5
check('Porosity', porosity, 0.48)
void_ratio = volume_voids / volume_tailings_solids # eq. 6
check('Void ratio', void_ratio, 0.93)
volume_interstitial_water = volume_voids * tailings_saturation # eq. 7
check('Volume interstitial water', volume_interstitial_water, 69000000)
mass_interstitial_water = volume_interstitial_water * density_water # eq. 8
check('Mass interstitial water', mass_interstitial_water, 69000000)
moisture_content = mass_interstitial_water / mass_tailings_stored # eq. 9
check('Moisture content', moisture_content, 0.34)
density_bulk = (mass_tailings_stored + mass_interstitial_water) / volume_tailings_deposit # eq. 10
check('Bulk density', density_bulk, 1.9)
mass_free_water = volume_supernatent * density_water # eq. 11
check('Mass supernatent', mass_free_water, 10000000)

In Line 10, the volume of supernatent was defined and then in Line 40, the mass of supernatent (free water) is calculated.
Estimating the solids content of the breach outflow
The next step is to estimate the solids content of the tailings leaving the breach. In reality, that's not something that can be predicted easily, so an assumption is made, to be consistent with the failure scenario being modeled.
In this example, the %solids of the tailings lost from the breach is estimated at 55%.
Predict the Dam Breach Initial Flood Wave Outflow Volume
The final step is to calculate the total volume of the flow from the breach based on the solids density
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 
def check(label, calculated_value, paper_val):
print '{}  difference={:.3f}%'.format(label, (calculated_value  paper_val)/paper_val)
mass_tailings_stored = 200000000 # t
average_dry_density = 1.4 # t/m3
specific_gravity = 2.7
tailings_saturation = 1.0 # % as dec fraction, so 1.0 == 100%
density_water = 1.0 # t/m3
volume_supernatent = 10000000 # m3
breach_solids_content = 0.55 # % as decimal fraction so 0.55 = 55%
percent_supernatent_lost = 1.0 # 100% as a decimal fraction
volume_tailings_deposit = mass_tailings_stored / average_dry_density # eq. 2
check( 'Volume tailings deposit', volume_tailings_deposit, 143000000)
volume_tailings_solids = mass_tailings_stored / specific_gravity # eq. 3
check('Volume tailings solids', volume_tailings_solids, 74000000)
volume_voids = volume_tailings_deposit  volume_tailings_solids # eq. 4
check('Volume voids', volume_voids, 69000000)
porosity = volume_voids / volume_tailings_deposit # eq. 5
check('Porosity', porosity, 0.48)
void_ratio = volume_voids / volume_tailings_solids # eq. 6
check('Void ratio', void_ratio, 0.93)
volume_interstitial_water = volume_voids * tailings_saturation # eq. 7
check('Volume interstitial water', volume_interstitial_water, 69000000)
mass_interstitial_water = volume_interstitial_water * density_water # eq. 8
check('Mass intersitial water', mass_interstitial_water, 69000000)
moisture_content = mass_interstitial_water / mass_tailings_stored # eq. 9
check('Moisture content', moisture_content, 0.34)
density_bulk = (mass_tailings_stored + mass_interstitial_water) / volume_tailings_deposit # eq. 10
check('Bulk density', density_bulk, 1.9)
mass_free_water = volume_supernatent * density_water # eq. 11
check('Mass supernatent', mass_free_water, 10000000)
mass_solids_mobilised = mass_free_water / ((1.0 / breach_solids_content)1  moisture_content) # eq. 12 reorgnaised to solve for mass of solids mobilised.
check('Mass solids mobiised', mass_solids_mobilised, 21000000)
check_moisture = moisture_content < (1.0/breach_solids_content)  1
if not check_moisture:
# not valid result
print 'Moisture content in the runout solids is too low'
end()
# mass mobilised < mass_tailings_stored
if mass_solids_mobilised > mass_tailings_stored:
print 'More tailings mobilised than stored'
end()
mass_interstitial_water_mobilised = ((mass_solids_mobilised / average_dry_density)  (mass_solids_mobilised/specific_gravity)) * tailings_saturation * density_water # eq. 14
check('Mass interstitial water mobilised', mass_interstitial_water_mobilised, 7000000)
volume_tailings_mobilised = (mass_solids_mobilised + mass_interstitial_water_mobilised) / density_bulk # eq. 15
check('Volume tailings mobilised', volume_tailings_mobilised, 15000000)
breach_outflow_volume = volume_supernatent * percent_supernatent_lost + volume_tailings_mobilised # eq. 16
check('Beach outflow volume', breach_outflow_volume, 25000000)
percent_lost = breach_outflow_volume / (volume_supernatent + volume_tailings_deposit) # eq. 17
check('Percent lost', percent_lost, 0.16)

The breach solids content of 55% is added as a variable in Line 11. Another variable, percent_supernatent_lost was added in Line 12. While not in the paper, this allows to decide how much (from 0  100%) of the pond supernantent water flows through the breach.
The mass of solids mobilised is then calculated in Line 44.
Some checks are done in Lines 48 to 55 to make sure that the volumes don't use more water or solids than are available. If one of these conditions is True, then the script will terminate at the end function.
The end function is a Muk3D specific command to terminate scripts gracefully. The common way to do this in normal Python is to use the sys.exit function. Running this function in a Muk3D Python script may cause the entire program to terminate. 
In Lines 60  69 the remaining volumes are calculated and the percentage of tailings lost from the TSF is calculated.
Enhancing the function
Now that we have a working script, we can make some modifications to it to make it a little more useful:
 Wrap the code in a function to make it easier to execute from different places.
 Run the code for a range of values and plot the results.
 Run the code varying 2 input values at a time and plotting the results.
Creating a function
By taking our script and creating a function around it, it makes it easier to use the code. If we want to call it multiple times in a script, then we just call the function, rather than duplicating the entire script. It also means that if we want to change something in the calculation, we update the function. If we had copied the code around, we'd have to make the change/fix in each segment of code, leading to a greater chance that stuff will get messed up/
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 
mass_tailings_stored = 200000000 # t
average_dry_density = 1.4 # t/m3
specific_gravity = 2.7
tailings_saturation = 1.0 # % as dec fraction, so 1.0 == 100%
density_water = 1.0 # t/m3
volume_supernatent = 10000000 # m3
breach_solids_content = 0.55 # % as decimal fraction so 0.55 = 55%
percent_supernatent_lost = 1.0 # 100% as a decimal fraction
class NotEnoughWaterError(Exception):
"""
Exception to be raised when the volume of water in the runout exceeds facility free+interstitial water
"""
pass
class NotEnoughTailingsError(Exception):
"""
Exception to be raised when the volume of tailings released is greater than TSF volume.
"""
pass
def estimate_runout_volume(mass_tailings_stored,
average_dry_density,
specific_gravity,
tailings_saturation,
density_water,
volume_supernatent,
breach_solids_content,
percent_supernatent_lost):
"""
Estimates the runout volume from a TSF based on the method presented by
D.Fontaine from KP.
"Tailings mobilization estimates for dam breach studies", Daniel Fontaine, P.Eng, Violeta Martin, Ph.D., P.Eng,
Proceedings Tailings and Mine Waste 2015
Vancouver, BC, October 26 to 28, 2015
"""
volume_tailings_deposit = mass_tailings_stored / average_dry_density # eq. 2
volume_tailings_solids = mass_tailings_stored / specific_gravity # eq. 3
volume_voids = volume_tailings_deposit  volume_tailings_solids # eq. 4
porosity = volume_voids / volume_tailings_deposit # eq. 5
void_ratio = volume_voids / volume_tailings_solids # eq. 6
volume_interstitial_water = volume_voids * tailings_saturation # eq. 7
mass_interstitial_water = volume_interstitial_water * density_water # eq. 8
moisture_content = mass_interstitial_water / mass_tailings_stored # eq. 9
density_bulk = (mass_tailings_stored + mass_interstitial_water) / volume_tailings_deposit # eq. 10
mass_free_water = volume_supernatent * density_water # eq. 11
mass_solids_mobilised = mass_free_water / ((1.0 / breach_solids_content)1  moisture_content) # eq. 12 reorgnaised to solve for mass of solids mobilised.
check_moisture = moisture_content < (1.0/breach_solids_content)  1
if not check_moisture:
raise NotEnoughWaterError()
# mass mobilised < mass_tailings_stored
if mass_solids_mobilised > mass_tailings_stored:
raise NotEnoughTailingsError()
mass_interstitial_water_mobilised = ((mass_solids_mobilised / average_dry_density)  (mass_solids_mobilised/specific_gravity)) * tailings_saturation * density_water # eq. 14
volume_tailings_mobilised = (mass_solids_mobilised + mass_interstitial_water_mobilised) / density_bulk # eq. 15
breach_outflow_volume = volume_supernatent * percent_supernatent_lost + volume_tailings_mobilised # eq. 16
percent_lost = breach_outflow_volume / (volume_supernatent + volume_tailings_deposit) # eq. 17
return breach_outflow_volume, percent_lost
runout_volume, percent_lost = estimate_runout_volume(mass_tailings_stored,
average_dry_density,
specific_gravity,
tailings_saturation,
density_water,
volume_supernatent,
breach_solids_content,
percent_supernatent_lost)

In Line 11 and Line 18, two new Exceptions are created by subclassing Pythons Exception class. These will be raised if any of the check conditions are met when calculating the runout volume.
In Line 25 a function is defined that wraps the code from the previous iteration of the script. The check statements have been removed for clarity. The function returns the breach output volume and the percentage of tailings lost. Other values could be returned if desired.
In Line 56, a check of the moisture in the runout tailings vs the impoundment is done, and if it fails, the NowEnoughWater exception is raised. In Line 60 a total volume check is done and the NotEnoughTailingsError is raised.
Running sensitivity
Now that we can run the volumetric analysis using a function, its easy to now wrap the function call in a loop to do sensitivity analysis.
Supernatent pond volume
The first sensitivity will be for supernantent pond volume. For a series of pond volumes, the breach outflow volume will be calculated.
In Line 73, we create an empty Python list in which to store results.
In this example, in Line 75, the variable volume_supernatent is replaced with successive values from the list (defined with square braces). In Line 87, the results list is printed to Muk3D's output window.
72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 
results = []
for volume_supernatent in [0, 5000000, 7500000, 10000000, 12500000, 15000000]:
runout_volume, percent_lost = estimate_runout_volume(mass_tailings_stored,
average_dry_density,
specific_gravity,
tailings_saturation,
density_water,
volume_supernatent,
breach_solids_content,
percent_supernatent_lost)
results.append(runout_volume)
print results

Plotting the results
We can create a plot of the results using the ui.plot.PlotWindow class from Muk3D's Python api. This plot window provivdes an interface to the Matplotlib plotting library for Python. A basic tutorial on using this library within Muk3D can be found in this tutorial.
In the following code, the function developed earlier is omitted from the code snippets, but is still in the executed Python scripts. 
To use the PlotWindow, it must be imported from the Muk3D Python API (from muk3d.ui.plot) and a new Python class is created subclassing this PlotWindow. Within the PlotWindow class, there is a method called create_plot which is overridden with code used to create the plot.
73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 
from muk3d.ui.plot import PlotWindow
class BasicPlot(PlotWindow):
def create_plot(self, figure):
results = []
pond_volumes = [0, 5000000, 7500000, 10000000, 12500000, 15000000]
for volume_supernatent in pond_volumes:
runout_volume, percent_lost = estimate_runout_volume(mass_tailings_stored,
average_dry_density,
specific_gravity,
tailings_saturation,
density_water,
volume_supernatent,
breach_solids_content,
percent_supernatent_lost)
results.append(runout_volume)
axes = figure.subplots()
axes.plot(pond_volumes, results, marker='x')
axes.set(xlabel='Pond volume (m$^3$)', ylabel='Breach outflow volume (m$^3$)',
title='Sensitivity of Breach Outflow Volume to Pond Volume')
axes.grid()
figure.tight_layout()
axes.set_xlim(0, max(pond_volumes))
axes.set_ylim(0, max(results))
plot = BasicPlot()
plot.show()

Line 75: The class is created and is subclassing the PlotWindow class.
Line 76: The class method create_plot is overridden, and takes the arguments self (an instance of the BasicPlot class) and an argument called figure, which is a Matplotlib figure instance. A figure is a blank page on which plot axes can be added.
Line 77: An empty list is created. This will be used to store the runout volume results for each trial.
Line 78: A list of pond volumes to loop through is created.
Line 79: A for loop is used to loop through each pond volume and create a variable called volume_supernatent.
Line 80: The estimate_runout_volume function is called with the basic parameters and the overridden volume_supernatent to calculate the runout_volume.
Line 88: The resulting runout_volume value is added to the results list using the append method.
Line 90: A Matplotlib Axes object is created and added to the figure using figure.subplots()
Line 91: A line plot is created using the axes.plot method. The xvalues are the pond_volumes, yvalues are the results, and an x is used as the point marker.
Line 92: The axes.set method is used to set the plots xaxis label, yaxis label, and title. The odd looking text $^3$ is Matplotlib's method of creating formulas/formatting text.
Line 94: A grid is added to the plot.
Line 95: The figure.tight_layout adjusts the plot layout so the plot and labels are clearly visible.
Line 96: The xaxis range is set to be between 0 and the largest pond volume.
Line 97: The yaxis range is set to be between 0 and the largest runout volume.
Line 99: An instance of the BasicPlot class is created.
Line 100: The plot is shown. At this time the create_plot method is called.
Running different cases
Now that we can create a single plot, we can run sensitivity varying 2 inputs at a time.
75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 
class BasicPlot(PlotWindow):
def create_plot(self, figure):
axes = figure.subplots()
pond_volumes = [0, 5000000, 7500000, 10000000, 12500000, 15000000]
pond_volumes_Mm3 = [v/1000000.0 for v in pond_volumes]
max_volume = 0
for breach_solids_content in [0.25,0.35,0.45,0.55,0.65]:
results = []
for volume_supernatent in pond_volumes:
runout_volume, percent_lost = estimate_runout_volume(mass_tailings_stored,
average_dry_density,
specific_gravity,
tailings_saturation,
density_water,
volume_supernatent,
breach_solids_content,
percent_supernatent_lost)
results.append(runout_volume/1000000.0)
max_volume = max(max_volume, runout_volume)
axes.plot(pond_volumes_Mm3, results, marker='x', label='{:.0f}%'.format(breach_solids_content * 100.0))
axes.set(xlabel='Pond volume (Mm$^3$)', ylabel='Breach outflow volume (Mm$^3$)',
title='Sensitivity of Breach Outflow Volume to Pond Volume')
axes.grid()
axes.legend()
figure.tight_layout()
axes.set_xlim(0, max(pond_volumes_Mm3))
axes.set_ylim(0, max_volume/1000000.0)

Line 80: Since we want to display the volumes in Mm3 instead of m3, we've taken the volume array and divided each element by 1000000 to use as the labeling array when we plot the data.
Line 81: A variable has been created to track the maximum runout volume so we can use this for scaling the yaxis later on.
Line 83: A for loop is used to loop through a range of breach solids content values.
Line 85: A list is created to store the results, each time we change the breach_solids_content.
Line 86: A for loop is used to loop through all the pond volumes.
Line 95: The runout_volume is added to the results, but divided by 1000000 since we're presenting the results in Mm3.
Line 97: The max_volume variable is set to the largest of the current max_volume and the current runout_volume value.
Line 98: The results for a particular breach solids content are plotted. The label keyword sets the series label to the breach solids content value.
Line 103: A legend is added to the plot.
Lines 105/106: The axis limits are set.