< All Topics

More 2D stochastic slope modelling

In the previous article, a simple statistical model was created to generated different beach slopes on a daily basis.  These beach slopes were deposited using a simple 2D model, and based on many deposition runs, a concave beach slope was found to emerge.

In this tutorial, we’ve going to implement the statistical model outlined by Seddon et al in the paper below.

Seddon, KD, Pirouz, B & Fitton, TG 2015, ‘Stochastic beach profile modelling’, in R Jewell & AB Fourie (eds), Proceedings of the 18th International Seminar on Paste and Thickened Tailings, Australian Centre for Geomechanics, Perth, pp. 455-465.

The model

Beach slopes were calculated to be dependent on the mass flow rate and the solids concentration of the tailings slurry.

mceclip0.png

The distribution for the mass flow rate is described below.  The authors of the paper describe an asymmetric distribution, with the lower and upper tail SD’s described.  For the sake of simplicity in this article, we will just work with a normal distribution using the lower SD.

mceclip1.png

The distribution for the solids concentration is described below:

mceclip2.png

Creating the Normal Distributions

Line 1-4: Importing some of the modules we’ll be using.

Line 6: Creating a normal distribution using stats.norm for the solids concentration.

Line 7: Creating a normal distribution using stats.norm for the mass_flow_rate (using the lower SD – we’re not making an asymmetric distribution for now).

import numpy
import scipy.stats as stats
import scipy.optimize as optimize
import scipy.interpolate as interpolate

solids_concentration = stats.norm(loc=62.8, scale=3.14)
mass_flow_rate = stats.norm(loc=1250, scale=150)

Creating the slope lookup table

The next step is to create a 2D lookup table to interpolate the beach slope based on the solids_concentration and mass_flow_rate.  To do this the data from the source paper will be parsed to separate out the row labels (solids concentration), column labels (mass flow rate) and the 2D array of slope data.

Line 10-17: The values from the Table 4 from Seddon et al is entered as tab separated text (it was copy/pasted into Notepad++ from Excel, which by default will be tab separated).

Line 19: All the rows of the Beach_Slopes text is split by row.  In Python, this is done using the .split method of the text and using ‘\n’ as the split token.  \n is the newline character and is hidden at the end of each row of text.

Line 20: The x indices are contained in the first item in slope_rows.  Note that there is a blank column in the top-left (if you look at the data in Notepad++ and show all characters you’ll see a tab symbol at the start of the line). The slope_rows[0] is the first line in the split Beach_slopes, and its split by a tab character. All values after the first character (denoted by the [1:] notation) are then passed to the map function, which converts each token (number as text) into a floating point number.

mceclip3.png

Line 21: The first value in each row (except the first row) is extracted and converted into values for the solids concentration. The notation [r.split(‘\t’)[0] for r in slope_rows[1:]] is a shorthand way for creating a new list in Python and is called list comprehension.

Line 22: The slope data is extracted by taking the data from the second row in slope_rows and ignoring the first entry in each row (which is the solids concentration). This is a 2D array.

import numpy
import scipy.stats as stats
import scipy.optimize as optimize
import scipy.interpolate as interpolate

solids_concentration = stats.norm(loc=62.8, scale=3.14)
mass_flow_rate = stats.norm(loc=1250, scale=150)


Beach_Slopes = """	0.39	0.38	0.36	0.35	0.33	0.3	0.27
56.65	0.58	0.6	0.61	0.62	0.64	0.67	0.71
59.19	0.73	0.75	0.77	0.78	0.81	0.84	0.9
61.15	0.87	0.89	0.91	0.93	0.96	1.01	1.07
62.8	1.02	1.05	1.07	1.09	1.13	1.18	1.25
64.45	1.2	1.24	1.26	1.29	1.33	1.39	1.48
66.41	1.49	1.54	1.57	1.6	1.66	1.73	1.84
68.95	2.03	2.08	2.13	2.17	2.25	2.36	2.51"""

slope_rows = Beach_Slopes.split('\n')
slope_x_indices = map(float, slope_rows[0].split('\t')[1:])
slope_y_indices = map(float, [r.split('\t')[0] for r in slope_rows[1:]])
slope_data = [map(float, r.split('\t')[1:]) for r in slope_rows[1:]]

Now that we’ve got the solids concentration and mass flow rate indices and the 2D data for the lookup table, we can create a 2D lookup table using a class from scipy.interpolate called interp2d.

Line 24: The slope_interpolator is created with the x lookup values (mass flow rate), y lookup values (solids concentration) and the 2D slope data.

import numpy
import scipy.stats as stats
import scipy.optimize as optimize
import scipy.interpolate as interpolate

solids_concentration = stats.norm(loc=62.8, scale=3.14)
mass_flow_rate = stats.norm(loc=1250, scale=150)


Beach_Slopes = """	0.39	0.38	0.36	0.35	0.33	0.3	0.27
56.65	0.58	0.6	0.61	0.62	0.64	0.67	0.71
59.19	0.73	0.75	0.77	0.78	0.81	0.84	0.9
61.15	0.87	0.89	0.91	0.93	0.96	1.01	1.07
62.8	1.02	1.05	1.07	1.09	1.13	1.18	1.25
64.45	1.2	1.24	1.26	1.29	1.33	1.39	1.48
66.41	1.49	1.54	1.57	1.6	1.66	1.73	1.84
68.95	2.03	2.08	2.13	2.17	2.25	2.36	2.51"""

slope_rows = Beach_Slopes.split('\n')
slope_x_indices = map(float, slope_rows[0].split('\t')[1:])
slope_y_indices = map(float, [r.split('\t')[0] for r in slope_rows[1:]])
slope_data = [map(float, r.split('\t')[1:]) for r in slope_rows[1:]]

slope_interpolator = interpolate.interp2d(slope_x_indices, slope_y_indices, slope_data)

print slope_interpolator(0.38, 59.19), 0.75
print slope_interpolator(0.3, 62.8), 1.18

Now we can try to lookup values using mass flow rates and solids concentration drawn from the normal distributions.

Line 27: Use a for loop and run it 10 times.

Line 28: Draw a value from each normal distribution using the .rvs() method.  These are then used to look up a slope value from the slope_interpolator.

import numpy
import scipy.stats as stats
import scipy.optimize as optimize
import scipy.interpolate as interpolate


solids_concentration = stats.norm(loc=62.8, scale=3.14)
mass_flow_rate = stats.norm(loc=1250, scale=150)


Beach_Slopes = """	0.39	0.38	0.36	0.35	0.33	0.3	0.27
56.65	0.58	0.6	0.61	0.62	0.64	0.67	0.71
59.19	0.73	0.75	0.77	0.78	0.81	0.84	0.9
61.15	0.87	0.89	0.91	0.93	0.96	1.01	1.07
62.8	1.02	1.05	1.07	1.09	1.13	1.18	1.25
64.45	1.2	1.24	1.26	1.29	1.33	1.39	1.48
66.41	1.49	1.54	1.57	1.6	1.66	1.73	1.84
68.95	2.03	2.08	2.13	2.17	2.25	2.36	2.51"""

slope_rows = Beach_Slopes.split('\n')
slope_x_indices = map(float, slope_rows[0].split('\t')[1:])
slope_y_indices = map(float, [r.split('\t')[0] for r in slope_rows[1:]])
slope_data = [map(float, r.split('\t')[1:]) for r in slope_rows[1:]]

slope_interpolator = interpolate.interp2d(slope_x_indices, slope_y_indices, slope_data)

for i in range(10):
    print slope_interpolator(mass_flow_rate.rvs(), solids_concentration.rvs())

The output window should the interpolated slopes.

mceclip0.png

Running 2D deposition

The next step is to take the 2D deposition code from the previous tutorial and integrate this new slope model into it.

Line 27 & below: This was copy/pasted from the previous tutorial script.

Line 68: Beach slope is set to interpolate from the 2D lookup table using values drawn from the

import numpy
import scipy.stats as stats
import scipy.optimize as optimize
import scipy.interpolate as interpolate


solids_concentration = stats.norm(loc=62.8, scale=3.14)
mass_flow_rate = stats.norm(loc=1250, scale=150)


Beach_Slopes = """	0.39	0.38	0.36	0.35	0.33	0.3	0.27
56.65	0.58	0.6	0.61	0.62	0.64	0.67	0.71
59.19	0.73	0.75	0.77	0.78	0.81	0.84	0.9
61.15	0.87	0.89	0.91	0.93	0.96	1.01	1.07
62.8	1.02	1.05	1.07	1.09	1.13	1.18	1.25
64.45	1.2	1.24	1.26	1.29	1.33	1.39	1.48
66.41	1.49	1.54	1.57	1.6	1.66	1.73	1.84
68.95	2.03	2.08	2.13	2.17	2.25	2.36	2.51"""

slope_rows = Beach_Slopes.split('\n')
slope_x_indices = map(float, slope_rows[0].split('\t')[1:])
slope_y_indices = map(float, [r.split('\t')[0] for r in slope_rows[1:]])
slope_data = [map(float, r.split('\t')[1:]) for r in slope_rows[1:]]

slope_interpolator = interpolate.interp2d(slope_x_indices, slope_y_indices, slope_data)
    
# copy/pasted from previous tutorial       
terrain_length = 2000

weights = numpy.ones((terrain_length,))
weights[0] = 0.5
weights[-1] = 0.5

base_grade = 0.01  # 1% slope 

ground_array = numpy.zeros((terrain_length,))
x_coordinates = numpy.linspace(0,terrain_length-1,terrain_length)
gradient = x_coordinates * base_grade
ground_array = ground_array - gradient

from muk3d.geometry.polyline import PolyLine
from muk3d.view import add_to_scene

polyline = PolyLine()
polyline.new_line_from_points(zip(x_coordinates,  [0,] * terrain_length,  ground_array))

add_to_scene('tailings profile', polyline, overwrite_existing=True)

def pour_tailings(discharge_elevation, beach_slope, current_ground):
    tailings_surface = discharge_elevation - x_coordinates * beach_slope
    
    new_tailings_surface = numpy.maximum(current_ground, tailings_surface)
    area = numpy.sum((new_tailings_surface- current_ground) * weights)
    
    return area, new_tailings_surface

current_ground = ground_array.copy()
    
def optimise_z(discharge_elevation, beach_slope, current_ground, target_area):
    area, new_ground = pour_tailings(discharge_elevation, beach_slope, current_ground)
    return (area - target_area)/target_area

target_area = 100

deposited_area = 0
for i in range(100):
    z_guess = max(current_ground) + 1
    slope = slope_interpolator(mass_flow_rate.rvs(), solids_concentration.rvs())
    final_z = optimize.newton(optimise_z, z_guess, tol=0.001, args=(slope, current_ground, target_area))

    deposition_2d_area, tailings_surface = pour_tailings(final_z, slope, current_ground)

    polyline.new_line_from_points(zip(x_coordinates,  [0,] * terrain_length,  tailings_surface))
    deposited_area += deposition_2d_area
    current_ground = tailings_surface
    
print 'deposited area', deposited_area

Run the script.  You’ll see something odd happen.

mceclip2.png

The tailings beach is incredibly sleep.  This is because the beach slopes are as percentages (not as a decimal fraction).  To fix this in Line 68 when we get the slope, we need to divide it by 100 to get percent as a decimal fraction.  Make this change and run again. It should look much more reasonable.

mceclip3.png

Final script

import numpy
import scipy.stats as stats
import scipy.optimize as optimize
import scipy.interpolate as interpolate


solids_concentration = stats.norm(loc=62.8, scale=3.14)
mass_flow_rate = stats.norm(loc=1250, scale=150)


Beach_Slopes = """	0.39	0.38	0.36	0.35	0.33	0.3	0.27
56.65	0.58	0.6	0.61	0.62	0.64	0.67	0.71
59.19	0.73	0.75	0.77	0.78	0.81	0.84	0.9
61.15	0.87	0.89	0.91	0.93	0.96	1.01	1.07
62.8	1.02	1.05	1.07	1.09	1.13	1.18	1.25
64.45	1.2	1.24	1.26	1.29	1.33	1.39	1.48
66.41	1.49	1.54	1.57	1.6	1.66	1.73	1.84
68.95	2.03	2.08	2.13	2.17	2.25	2.36	2.51"""

slope_rows = Beach_Slopes.split('\n')
slope_x_indices = map(float, slope_rows[0].split('\t')[1:])
slope_y_indices = map(float, [r.split('\t')[0] for r in slope_rows[1:]])
slope_data = [map(float, r.split('\t')[1:]) for r in slope_rows[1:]]

slope_interpolator = interpolate.interp2d(slope_x_indices, slope_y_indices, slope_data)
    
# copy/pasted from previous tutorial       
terrain_length = 2000

weights = numpy.ones((terrain_length,))
weights[0] = 0.5
weights[-1] = 0.5

base_grade = 0.01  # 1% slope 

ground_array = numpy.zeros((terrain_length,))
x_coordinates = numpy.linspace(0,terrain_length-1,terrain_length)
gradient = x_coordinates * base_grade
ground_array = ground_array - gradient

from muk3d.geometry.polyline import PolyLine
from muk3d.view import add_to_scene

polyline = PolyLine()
polyline.new_line_from_points(zip(x_coordinates,  [0,] * terrain_length,  ground_array))

add_to_scene('tailings profile', polyline, overwrite_existing=True)

def pour_tailings(discharge_elevation, beach_slope, current_ground):
    tailings_surface = discharge_elevation - x_coordinates * beach_slope
    
    new_tailings_surface = numpy.maximum(current_ground, tailings_surface)
    area = numpy.sum((new_tailings_surface- current_ground) * weights)
    
    return area, new_tailings_surface

current_ground = ground_array.copy()
    
def optimise_z(discharge_elevation, beach_slope, current_ground, target_area):
    area, new_ground = pour_tailings(discharge_elevation, beach_slope, current_ground)
    return (area - target_area)/target_area

target_area = 100

deposited_area = 0
for i in range(100):
    z_guess = max(current_ground) + 1
    slope = slope_interpolator(mass_flow_rate.rvs(), solids_concentration.rvs())/100.0
    final_z = optimize.newton(optimise_z, z_guess, tol=0.001, args=(slope, current_ground, target_area))

    deposition_2d_area, tailings_surface = pour_tailings(final_z, slope, current_ground)

    polyline.new_line_from_points(zip(x_coordinates,  [0,] * terrain_length,  tailings_surface))
    deposited_area += deposition_2d_area
    current_ground = tailings_surface
    
print 'deposited area', deposited_area

Next steps

The next tutorial will look at how this can be simulated with 3D deposition.

Table of Contents