At the Paste 2015 conference in Cairns, we saw a paper presented called Stochastic beach profile modelling. In this paper, the authors developed a statistical model of thickener performance and equated this to beach slopes. On a daily basis, a volume of tailings was deposited drawing a beach slope from this statistical distribution. A total of 365 days of tailings were modelled, each day with a different linear slope, until a stack had been created that effectively exhibited non-linear slopes. At the time, this would have been a monumental undertaking given that it was done by manually entering each days worth of deposition.
The results of this modelling is shown in the figures below (taken from the paper).
In this article, we'll develop a 2D modelling tool for taking a range of linear slopes drawn from a statistical slope model and see how a non-linear beach appears. In subsequent articles we'll implement the models outlined by Seddon et al, and then apply it to 3D deposition modelling.
 |
The final script developed is linked at the bottom of this article. |
Simple stochastic slope model
The first step will be to create a simple model for beach slopes that draws a slope from a normal distribution. This will then be deposited on a simplified 2D terrain for a given volume (or in this case, cross section area) and repeated to form a 2D beach slope.
To do this, we'll use a few Python libraries:
- numpy: Array manipulation
- Scipy: Scientific libraries that includes statistical models and optimisation tools (for goalseeking beach elevation).
Importing libraries
Create a new script and import the following libraries.
1 2 3 |
import numpy
import scipy.stats as stats
import scipy.optimize as optimize
|
Creating the slope model
The next step is to create the stochastic slope model. This is done by creating a new Python class called StochasticTailings. This has two methods: get_slope() and get_volume(). In the StochasticTailings class these are not implemented and need to be redefined by derived classes.
Next we have defined a new class called NormalDistribution which is derived from the StochasticTailings base class. It takes a slope mean and standard deviation, and a discharge rate (/area). In the __init__ method, a normal distribution is created in Line 15.
The get_slope method of the base class is then overridden in Line 20, and it draws a value from the normal distribution in Line 21. In Line 22, we make sure that the slope doesn't fall outside the optional maximum and minimum bounds. We do this just to ensure that there is no possibility of a slope being returned that is negative (which would indicate uphill flow).
In Line 24, the get_volume method is overridden and returns the discharge_rate value provided in the __init__ method.
In Line 27, an instance of NormalDistribution is created with a mean slope of 2% (0.02) and a standard deviation of 0.5%. A discharge rate of 100 is used.
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 |
import numpy
import scipy.stats as stats
import scipy.optimize as optimize
class StochasticTailings:
def get_slope(self):
raise "Not implemented"
def get_volume(self):
raise "Not implemented"
class NormalDistribution(StochasticTailings):
def __init__(self, mean, sd, discharge_rate, min_slope=0.0001, max_slope=100):
self._pdf = stats.norm(loc=mean, scale=sd)
self._min = min_slope
self._max = max_slope
self._discharge_rate = discharge_rate
def get_slope(self):
p = self._pdf.rvs()
return min(max(p, self._min), self._max)
def get_volume(self):
return self._discharge_rate
tailings = NormalDistribution(0.02, 0.005, 100)
for i in range(10):
print tailings.get_slope()
|
Just to test this, in Line 29 we draw 10 slopes from the NormalDistribution instance and print it.
Creating the base 2D terrain
Next we will create a numpy array to use as the base terrain for deposition. To keep it simple we'll just create a linear terrain with a user defined gradient.
 |
Note that the previous segment of the script is not shown in the code below. |
Line 29: Create a variable to hold the terrain length in m/ft. This terrain segment is 2km long.
Line 31-33: An array is created with weights that are used when calculating the area between 2 lines (e.g. the ground and a new tailings surface.) The value of all but the end values of this array is 1, with the ends set to 0.5. The numpy function ones creates an array that is terrain_length long and all values initialised to 1.
Line 35: Create a variable with the base grade for the terrain. This is entered as a percent as a decimal fraction.
Line 37: A variable called ground_array is created to hold the elevation of the terrain. It is initially set at elevation 0, but is modified in the following lines to reflect the base_grade.
Line 38: An array of x-coordinates (distance along the ground) is created using numpy's linspace function. It takes a start distance, end distance, and number of points.
Line 39: The x coordinate array is multiplied by the gradient to create a sloping line.
Line 40: The gradient is subtracted from the ground array to create our base ground.
29 30 31 32 33 34 35 36 37 38 39 40 |
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
|
Displaying the line
If we run the script, we won't see anything. So now we need to add the terrain line to the 3D window. This is done using the PolyLine geometry object from the Muk3D API.
Line 42: The PolyLine class is imported from the muk3d.geometry.polyline module.
Line 43: The function add_to_scene is imported from the muk3d.view module. This is used to add a geometry object to the 3D window.
Line 45: An instance of PolyLine is created.
Line 46: A new line is added to the PolyLine instance. The x-values are the x_coordinates, y-values are 0, and z-values are the ground elevation array. The zip function takes these 3 arrays and turns them into a list of 3D points.
Line 48: The polyline is added to the 3D window, with the name tailings profile. The keyword argument overwrite_existing is set to True to the user isn't asked to overwrite the layer each time the script is run.
29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 |
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)
|
Pour tailings
Now we'll create a new line to represent some poured tailings. A function called pour_tailings is defined that takes a beach slope and discharge elevation, and creates a 2D line at the beach slope and merges it with the supplied base_ground array.
Line 50: The function pour_tailings is defined with the arguments of discharge_elevation, beach_slope, and the current ground array.
Line 51: An array representing the tailings surface is created by subtracting the x_coordinates x beach_slope from the discharge_elevation.
Line 53: The tailings_surface is merged into the current_ground by taking the maximum elevation at each point from the 2 arrays.
Line 54: The area between the current_ground and the new_tailings_surface is calculated.
Line 56: The cross section area and the new_tailings_surface array are returned.
Line 58: A variable called current_ground is created which contains a copy of the ground_array.
Line 60: Tailings is poured at an elevation of 5, with a slope drawn from the StochasticTailings model, and the current_ground array is supplied. The two return values are assigned to the deposition_2d_area, and tailings_surface variables.
Line 62: The new tailings surface is added as a new line to the PolyLine instance and will be visible in the 3D window.
48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 |
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()
deposition_2d_area, tailings_surface = pour_tailings(5, tailings.get_slope(), current_ground)
polyline.new_line_from_points(zip(x_coordinates, [0,] * terrain_length, tailings_surface))
|
The results are shown below.
Optimise tailings elevation for a given cross section area
In the example above the tailings is deposited from a given elevation, but we want to be able to deposit given areas of tailings instead. To do this, we need to optimise the discharge elevation until the target area is calculated. To do this we'll use the Scipy.optimize module. Within this module there is a function called newton which will optimize a function using Newton-Raphson method.
The newton function takes the following arguments:
- a function to minimise. The function needs to take at least one argument (the optimisation variable, in this case the discharge_elevation) and it must return the error between the current value of the function and the target value. The discharge_elevation will be modified until the error is below a tolerance.
- a guess for the initial value of the variable (discharge_elevation).
- a tolerance to optimise to.
- an optional list of arguments to pass to the optimise function (e.g. target area, beach slope, ground array)
Line 58 - 60: A function called optimise_z is created that will be passed to the newton function. This function takes a discharge_elevation, beach_slope, current_ground, and target_area, calls the pour_tailings function, and then calculates the error between the returned area and the target_area. This error value is returned and newton uses this to optimise the discharge_elevation so the error value is less than the required tolerance.
Line 62: A target area is defined.
Line 63: A beach slope is drawn from the NormalDistribution tailings model. We do this instead of using the get_slope() function when we want a slope since each time the get_slope() function is called it will return a different value.
Line 64: The newton function is called with the following arguments: optimise_z function, 5 (an estimate of the discharge_elevation), the keyword argument tol with the desired tolerance, and the keyword argument args that has the remaining arguments that the optimise_z function accepts - the slope, current_ground array and the target_area. The return value for this is the final discharge_elevation that gave an error less than the desired tolerance.
Line 66: Since we only know the required discharge elevation and don't have the tailings line, we run the pour_tailings function with the discharge elevation of final_z.
Line 68: Print the deposition_2d_area to the output window to verify that its close to the desired target.
Line 69: The new tailings surface is added to the PolyLine instance and shown in the 3D window.
48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 |
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
slope = tailings.get_slope()
final_z = optimize.newton(optimise_z, 5, tol=0.001, args=(slope, current_ground, target_area))
deposition_2d_area, tailings_surface = pour_tailings(final_z, slope, current_ground)
print 'Deposited area', deposition_2d_area
polyline.new_line_from_points(zip(x_coordinates, [0,] * terrain_length, tailings_surface))
|
Running a sequence
Now that we can optimise elevation for a given cross section area, we can run a sequence of runs, each time drawing a new beach slope from the NormalDistribution tailings model.
Line 66: A variable to track deposited area is created and initialised to zero.
Line 67: A for loop is created and set to execute 100 times.
Line 68: A guess for the initial tailings elevation is created based on the maximum elevation of the current_ground and adding 1 to it. This ensures that the initial guess will always be above ground. If it was lower than the ground the newton function may not locate the correct discharge_elevation.
Line 69: The beach slope is drawn from the NormalDistribution.
Line 70: The newton function is used to optimise the final_z value.
Line 72: The new tailings_surface is created using the final_z value.
Line 74: The tailings_surface is added to the PolyLine instance.
Line 75: The incremental deposition area is added to the deposited_area variable.
Line 76: The current_ground is redefined to be the current tailings_surface.
59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 |
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 = tailings.get_slope()
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
|
The resulting polyline is shown below.
The 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 78 |
import numpy
import scipy.stats as stats
import scipy.optimize as optimize
class StochasticTailings:
def get_slope(self):
raise "Not implemented"
def get_volume(self):
raise "Not implemented"
class NormalDistribution(StochasticTailings):
def __init__(self, mean, sd, discharge_rate, min_slope=0.0001, max_slope=100):
self._pdf = stats.norm(loc=mean, scale=sd)
self._min = min_slope
self._max = max_slope
self._discharge_rate = discharge_rate
def get_slope(self):
p = self._pdf.rvs()
return min(max(p, self._min), self._max)
def get_volume(self):
return self._discharge_rate
tailings = NormalDistribution(0.02, 0.005, 100)
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 = tailings.get_slope()
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 step is to implement the 2D lookup table outlined by Seddon et al in this tutorial.