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 final script developed is linked at the bottom of this article. |
The model
Beach slopes were calculated to be dependent on the mass flow rate and the solids concentration of the tailings slurry.
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.
The distribution for the solids concentration is described below:
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).
1 2 3 4 5 6 7 |
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.
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.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 |
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.
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 |
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.
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 |
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.
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 solids_concentration and mass_flow_rate distributions.
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 |
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.
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.
Final script
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 |
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.