# 2D beach slope stochastic modelling

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.

```
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.

```
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()
```

ust 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.

```
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.

```
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.

```
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.

```
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**.

```
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

```
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.