Introduction
The safety assessment of levees is currently carried out by dividing the levee into sections that are as much as possible homogeneous in terms of geometry and stratigraphy. Then, for each dike section, one dike profile and stratigraphy are chosen as the representative. This approach is time consuming and subjective in the way the representative cross-section is identified. Furthermore, the division of the dike into sections, whereby one cross section is calculated as the representative, can lead to assessing parts of a dike section unfavorably. The research topic of my bachelor thesis involved an automated calculation of cross sections at an interval of 25 meters along a 4 km long dike stretch. The high-density calculations were compared with the current approach, where only one representative cross-section is calculated.

The automated approach adopted in the research proved to be helpful in determining dike sections with more homogeneous characteristics than the ones that would result from the standard approach. Instead of determining the representative cross section based on assumptions on the most critical soil layering and geometry, the automated approach allows us to calculate all the cross sections exploiting all the available data. An eventual division of the dike in segments can be done at the end of the process, basing the decision on the calculated Factor of Safety, instead of on initial assumptions.

In this research I used GEOlib, a powerful geotechnical Python API that supports an automated process, to write the input file for D-Stability. I also used the D-Stability Console versions to execute the generated files.

In this article I will explain how Geolib can be implemented for automated dike calculation.

Data preprocessing
The geometries of the cross sections were retrieved from the AHN data using a python script that accesses the PDOK API. The location of the characteristic points (crest, berm, outer and inner toe) was visually determined for every cross-section and saved in a json file. The water levels per dike cross-section were generated using standpipe measurements and the previously determined characteristic points. The stratigraphy can be obtained by using the cpt-model and saving the output into a DataFrame. In total, the data relative to 117 cross sections were preprocessed.

Base model
Now let us start with the interesting part, we first import all the needed modules from geolib and Python package Shapely:

from geolib.models.dstability.dstability_model import DStabilityModel
from geolib.geometry.one import Point as DPoint
from geolib.soils import Soil, MohrCoulombParameters
from geolib.soils.soil import ShearStrengthModelTypePhreaticLevel
from geolib.models.dstability.states import DStabilityStatePoint, DStabilityStateLinePoint, DStabilityStress
from geolib.models.dstability.analysis import DStabilityBishopAnalysisMethod, DStabilityCircle, DStabilityBishopBruteForceAnalysisMethod, DStabilitySearchGrid, DStabilitySlipPlaneConstraints
from shapely.geometry import LineString, Point, Polygon
from pathlib import Path

We start defining our D-Stability object by initiating the model.

dm = DStabilityModel()

We can then continue setting the model data. We can begin by importing the soil types from a preprocessed file, in this case from a json file containing the soil data.

The example below shows how this can be done for the volumetric weight. The same structure is used for the other soil parameters.

for soil_name, soil_props in soil_properties_dict.items():
    soil = Soil()
    soil.name = soil_name
    soil.code = soil_name
    soil.soil_weight_parameters.saturated_weight.mean = soil_props["VolumetricWeightBelowPhreaticLevel"]
    soil.soil_weight_parameters.unsaturated_weight.mean = soil_props["VolumetricWeightAbovePhreaticLevel"]

    soil_id = Dm.add_soil(soil)

The next step is creating the layering for the dike and choosing the soil type for them. This can be done with the Python library Shapely by intersecting the geometry with the stratigraphy. After creating the layers, we can store the coordinates of each polygon into a list and add it to the model. When adding the polygons make sure that the starting point of a polygon that you add is not the same as the last point of the polygon. This can be achieved by looping until -1 and not to the end (credits to Rob van Putten).

layer_ids = []
for polygon in list_polygons:
    coords = list(polygon.exterior.coords)
    points = [DPoint(x=point[0], z=point[1]) for point in coords[:-1]]
    layer_id = dm.add_layer(points, polygon.layer)
    layer_ids.append(layer_id)

The water levels and state conditions can also be added to the base model with Geolib and a few lines of code.

phreatic_line_daily = list(waterlevels[0].coords)
phreatic_line_daily_id = dm.add_head_line(points=[ DPoint(x=point[0], z=point[1]) for point in phreatic_line_daily],
                                            label="Phreatic daily condition", 
                                            is_phreatic_line=True)

Add a calculation stage
In D-Stability, construction stages are applied, they represent sequential moments in time. The order of the stages is necessary for the yield stress administration. In this research, two construction stages were implemented: daily condition and high water. The base model that we create will be our daily condition. In Geolib we can create a new stage as follows:

dm.add_stage(label= 'High water', notes=None, set_current=True)

By adding a new stage, the geometry and materials are copied from the previous stage, in our case from the daily condition. Aspects that might affect the state like water pressures and loads are not copied to the new stage and need to be redefined.

Choosing the calculation method and defining the search Area
So far, we have created our model with two stages: daily condition and high water. To run a fullyautomated calculation, we first need to find a way to define the parameters for the search area to find the sliding surface with the minimum factor of safety. In my thesis I used the Bishop method. For this method, the search area is defined by a rectangular grid of centers and a tangent scope.

The location of the rectangular grid is determined automatically using a Python script that uses the characteristic points, the thickness of the soft layers and the width of the dike as an input. The scope of the tangent lines were also determined automatically by considering the thickness of the soft layers and the surface line. We can further use the characteristic points (the location of the inner and outer crest) to constrain the entry point of our slip plane in the x-direction. After we determine the parameters for the search area, we can save them into variables.

The Bishop analysis method can be set as follows:

dm.set_model(
DStabilityBishopBruteForceAnalysisMethod(
    search_grid=DStabilitySearchGrid(
        bottom_left= DPoint(x=bottom_left_x, z=bottom_left_y),
        number_of_points_in_x = points_horizontal,
        number_of_points_in_z = points_verticale,
        space=number_of_spaces,
    ),

slip_plane_constraints=
                    DStabilitySlipPlaneConstraints 
(is_zone_a_constraints_enabled = True, 
                                                 x_left_zone_a = BUK_x,
                                                 width_zone_a = width_crest),
    bottom_tangent_line_z= bottom_tangent,
    number_of_tangent_lines=number_of_lines,
    space_tangent_lines=space_tangent_lines,
))

After we have set the analysis method we can write the input file by serializing the model and saving it to a file.

dm.serialize(Path("directory/filename.stix"))

The model is now ready to be executed:

Calculating using D-Stability console / parallel calculation
The generated files can be executed manually by opening them in D-Stability. Since we want to fully automate the process, we used the D-Stability Console version to execute the files. By using the Python multiprocessing module we can fully leverage multiple processors on a given machine. This means that we can run several calculations at once. To do this we first make a function to execute our files:

def calculate(input_file):
    #Initialize the model
    dm = DStabilityModel()

    #Parse existing D-Stability file  
    dm.parse(Path(input_file))

    # execute the file
    dm.execute() 

Parallel calculations can be run like this:

from multiprocessing import Pool, cpu_count

Input_file_1 = "directory/filename_1.stix"
Input_file_2 = "directory/filename_2.stix"

Input_files = [Input_file_1, Input_file_2]
tasks = [*zip(input_files)]

with Pool(cpu_count()) as pool:                       
    pool.starmap(calculate, iterable=tasks)

This method of executing multiple files at once made it possible to run 117 files in just under 4 minutes (using all 12 cpu-cores of my laptop :). For the thesis, I further optimized the search area after one iteration of calculation. This optimization is done by moving the search area to the center of the calculated slip plane and by increasing the resolution of both the center's grid and tangent lines. The method of determining the search area automatically is in its infancy and therefore further research is needed.

Conclusion
This article has shown how we can use Geolib in automated levee stability calculations.

I hope this opens doors for more applications of Geolib and in general, further automations in geotechnical engineering. If you have any questions or want to know more about it leave you can send an email to aW5mb0BjZW1zYnYubmw=.