Numerical integration for depth-dependent problems

In geotechnical engineering, we are often interested in how certain properties accumulate over depth. Examples are the calculation of vertical effective stress from effective unit weight, pile shaft resistance from unit skin friction, ...

These problems can be written as integral expressions:

When the soil is uniform, evaluating the integral analytically is still possible and the solution is straightforward. However, when the soil is layered, or the function f(z) shows complex variations with depth, this problem becomes hard to solve analytically. Furthermore, we are often interested to know the value of the integral as the depth increases (e.g. for knowing when the pile tip is deep enough to have sufficient shaft resistance).

Theory

We this need to have a calculation method for evaluating the integral which returns the numerical value of the integral at a large number of depths. To achieve this, we can approximate the function f(z) in a given interval (a, b) with a linear function as shown in Figure 1. The integral in that interval can then be expressed as:

When this approximation is performed over a large number (N) of intervals (i.e. a small depth increment), we get the following expression for the integral:

Note that we take into account in the notation that the size of the increments may vary.

This approximation becomes more accurate as the number of increments increases or, in other words, when the increment size decreases.

Implementation

We can implement the recipe from the equation above using a Python loop but we can also use SciPy's scipy.integrate.cumulative_trapezoid function (https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.cumulative_trapezoid.html ).

This function will apply the trapezoid rule specified array of f(z) and z values to return the cumulative quantity with depth.

When looking at the function documentation, we can see that the following arguments are expected:

- y: These are the values of f(z)

- x: These are the values of z

- dx: This can be used instead of x if the increment size is equal for all depths

- axis: For one-dimensional arrays, this is just left at the default

- initial: If an initial value (e.g. non-zero vertical effective stress) is needed, this can be specified here.

Basic example: Vertical effective stress profile

Let's consider a 10m thick layer of soil with a uniform effective unit weight of 9kN/m3 and a hydrostatic pore pressure distribution. When we want to calculate the vertical effective stress, we need to evaluate the following integral:

It is easy to see that the vertical effective stress function is 9z and that the vertical effective stress at 10m depth is this 90kPa. We can also calculate this with the cumulative_trapeziod function.

Let's create a Numpy array for the depths with 50 equally spaced depths between 0 and 10m (using np.linspace). We can also create an array with the same length for the effective unit weight (using np.ones).

import numpy as np
# 50 equally spaced depths between 0m and 10m
z = np.linspace(0, 10, 50) 
# Numpy array with 50 elements, each with value 9
gamma_eff = 9 * np.ones(50)         

Calculating the vertical effective stress can then be done using the analytical expression. Note that NumPy performs the calculation in a vectorized manner, multiplying each element of effective unit weight by the corresponding depth and storing the results in a Numpy array. This is very effective when performing computations with arrays!

sigma_vo_eff = gamma_eff * z
sigma_vo_eff
OUTPUT: array([ 0.        ,  1.83673469,  3.67346939,  5.51020408,  7.34693878,
        9.18367347, 11.02040816, 12.85714286, 14.69387755, 16.53061224,
       18.36734694, 20.20408163, 22.04081633, 23.87755102, 25.71428571,
       27.55102041, 29.3877551 , 31.2244898 , 33.06122449, 34.89795918,
       36.73469388, 38.57142857, 40.40816327, 42.24489796, 44.08163265,
       45.91836735, 47.75510204, 49.59183673, 51.42857143, 53.26530612,
       55.10204082, 56.93877551, 58.7755102 , 60.6122449 , 62.44897959,
       64.28571429, 66.12244898, 67.95918367, 69.79591837, 71.63265306,
       73.46938776, 75.30612245, 77.14285714, 78.97959184, 80.81632653,
       82.65306122, 84.48979592, 86.32653061, 88.16326531, 90.        ])        

The same can be achieved using the Scipy function cumulative_trapezoid by filling in the arguments appropriately.

from scipy.integrate import cumulative_trapezoid        

cumulative_trapezoid(y=gamma_eff, x=z)
OUTPUT: array([ 1.83673469,  3.67346939,  5.51020408,  7.34693878,  9.18367347,
       11.02040816, 12.85714286, 14.69387755, 16.53061224, 18.36734694,
       20.20408163, 22.04081633, 23.87755102, 25.71428571, 27.55102041,
       29.3877551 , 31.2244898 , 33.06122449, 34.89795918, 36.73469388,
       38.57142857, 40.40816327, 42.24489796, 44.08163265, 45.91836735,
       47.75510204, 49.59183673, 51.42857143, 53.26530612, 55.10204082,
       56.93877551, 58.7755102 , 60.6122449 , 62.44897959, 64.28571429,
       66.12244898, 67.95918367, 69.79591837, 71.63265306, 73.46938776,
       75.30612245, 77.14285714, 78.97959184, 80.81632653, 82.65306122,
       84.48979592, 86.32653061, 88.16326531, 90.        ])        

We can see that the result is exactly the same, except that the first element (zero vertical effective stress for z=0) is missing. This can be mitigated by specifying the initial keyword argument. In the absence of surcharge, the vertical effective stress is 0kPa at z =0.

cumulative_trapezoid(y=gamma_eff, x=z, initial=0)
OUTPUT: array([ 0.        ,  1.83673469,  3.67346939,  5.51020408,  7.34693878,
        9.18367347, 11.02040816, 12.85714286, 14.69387755, 16.53061224,
       18.36734694, 20.20408163, 22.04081633, 23.87755102, 25.71428571,
       27.55102041, 29.3877551 , 31.2244898 , 33.06122449, 34.89795918,
       36.73469388, 38.57142857, 40.40816327, 42.24489796, 44.08163265,
       45.91836735, 47.75510204, 49.59183673, 51.42857143, 53.26530612,
       55.10204082, 56.93877551, 58.7755102 , 60.6122449 , 62.44897959,
       64.28571429, 66.12244898, 67.95918367, 69.79591837, 71.63265306,
       73.46938776, 75.30612245, 77.14285714, 78.97959184, 80.81632653,
       82.65306122, 84.48979592, 86.32653061, 88.16326531, 90.        ])        

Now, the direct evaluation using the analytical formula and the numerical integration using the trapezoid rule give the same result!

Gridding for geotechnical calculations

Numerical integration for geotechnical problems involves using the trapezoid rule on a depth axis which includes a given layering. To prevent writing functions with many if statements, we can make use of the functionality in the groundhog library.

groundhog contains the SoilProfile class which defines a layering. This is done by specifying the depths of the tops of each layer and the depths of layer bottoms. Soil parameters (both numerical and text-based) can then be added as additional columns. SoilProfile inherits from the Pandas DataFrame so all DataFrame functionality can also be used on SoilProfile objects. In addition to that, the SoilProfile class includes a number of functions which are specific for geotechnical profiles.

Here, we will declare a three-layer soil profile (sand overlying clay, with silt below). We can assign an effective unit weight to each layer.

from groundhog.general.soilprofile import SoilProfile
sp = SoilProfile({
    'Depth from [m]': [0, 5, 10],
    'Depth to [m]': [5, 10, 20],
    'Soil type': ['SAND', 'CLAY', 'SILT'],
    'Effective unit weight [kN/m3]': [10, 7, 8]})
sp        

The soilprofile contains all the information for the layering but to actually perform numerical integration, we need a finer grid. In groundhog this is achieved by creating a CalculationGrid from the SoilProfile object.

When creating the grid, nodes are added in between the layer interfaces using a regular node spacing (dz argument). When layer interfaces do not correspond to the nodal coordinates, the layer interfaces are added to the nodal coordinate list.

from groundhog.general.soilprofile import CalculationGrid
grid = CalculationGrid(soilprofile=sp, dz=0.3)        

The CalculationGrid has a nodes attribute which contains the nodal coordinates. When creating the CalculationGrid object from a SoilProfile object, the properties of the SoilProfile object are interpolated to the CalculationGrid object. The nodes attribute is a Pandas DataFrame.

We can print the nodes attribute.

grid.nodes        

For certain operations, it is useful to work with the elements between individual nodes. To achieve this, the CalculationGrid object also has an elements attribute which provides a listing of all elements with the top and bottom coordinate and the interpolated properties.

The elements attribute of a CalculationGrid object is a SoilProfile object so all methods defined for SoilProfile objects can also be used on this attribute.

We can print the elements attribute:

grid.elements        

The calculation of the vertical effective stress from effective unit weight can be done on the grid.nodes dataframe. We just need to call cumulative_trapezoid with the appropriate column names to supply the y and x arguments to the routine. To ensure that the output has the correct length, we also specify the initial value.

The result is assigned to the column 'Vertical effective stress [kPa]'.

grid.nodes['Vertical effective stress [kPa]'] = \
    cumulative_trapezoid(
        y=grid.nodes['Effective unit weight [kN/m3]'],
        x=grid.nodes['z [m]'], initial=0)        

To visualise the results, we can make use of the LogPlot in groundhog. This class allows creation of plots of geotechnical parameters with a miniature log of the stratigraphy alongside it.

from groundhog.general.plotting import LogPlot        

The LogPlot is defined by supplying a SoilProfile and determining how many panels there are. Here, we will create a panel for the effective unit weight and one for the vertical effective stress. The colors for the stratigraphic log are defined in a Python dictionary (fillcolordict) which maps each entry in the 'Soil type' column to a color.

Data is added to the plot using the .add_trace method. This method takes an array (or dataframe column) for the x and z arguments and the user needs to decide on which panel to plot the trace (panel_no argument). A name is also defined for each trace which is used in the legend.

The axes are set with the .set_xaxis (which takes the panel_no as argument) and .set_zaxis methods. Here, we set the title for each axis and define the range for the effective unit weight axis.

resultplot = LogPlot(
    soilprofile=sp, no_panels=2,
    fillcolordict={
        'SAND': 'yellow', 'CLAY': 'brown', 'SILT': 'orange'})
resultplot.add_trace(
    x=grid.nodes['Effective unit weight [kN/m3]'],
    z=grid.nodes['z [m]'],
    name='Effective unit weight',
    panel_no=1)
resultplot.add_trace(
    x=grid.nodes['Vertical effective stress [kPa]'],
    z=grid.nodes['z [m]'],
    name='Vertical effective stress',
    panel_no=2)
resultplot.set_xaxis(
    title='Effective unit weight [kN/m3]', panel_no=1,
    range=(0, 12))
resultplot.set_xaxis(
    title='Vertical effective stress [kPa]', panel_no=2)
resultplot.set_zaxis(title='z [m]')
resultplot.show()        

When inspecting the plot, we can see that the effective unit weight changes at layer transitions are not immediate. This is because a layer interface node cannot have two values of effective unit weigth. It will always take the value of the underlying layer.

This is why we can also calculate the vertical effective stress based on the element definitions. For this, we can use the depth_integration method of the SoilProfile object. We need to supply the column name a the parameter which will be integrated ('Effective unit weight [kN/m3]') and provide a name for the output. As the integration of a constant value will lead to a linear variation over the element, the output will be two columns, one giving the value at the top of the element ('Vertical effective stress from [kPa]') and one with the value of the bottom of the element ('Vertical effective stress to [kPa]').

grid.elements.depth_integration(
    parameter='Effective unit weight [kN/m3]',
    outputparameter='Vertical effective stress [kPa]')
grid.elements.head()        

We can use LogPlot again to visualise the results but instead of using add_trace, we use add_soilparameter_trace which just needs the name of the soil parameter to be plotted and the panel number. If a parameter shows a linear variation, we supply the name without 'from' or 'to'.

resultplot = LogPlot(
    soilprofile=grid.elements, no_panels=2,
    fillcolordict={
        'SAND': 'yellow', 'CLAY': 'brown', 'SILT': 'orange'})
resultplot.add_soilparameter_trace(
    parametername='Effective unit weight [kN/m3]', panel_no=1)
resultplot.add_soilparameter_trace(
    parametername='Vertical effective stress [kPa]', panel_no=2)
resultplot.add_trace(
    x=grid.nodes['Vertical effective stress [kPa]'],
    z=grid.nodes['z [m]'],
    name='Nodal vertical effective stress ',
    line=dict(dash='dot'),
    panel_no=2)
resultplot.set_xaxis(title='Effective unit weight [kN/m3]',
    panel_no=1, range=(0, 12))
resultplot.set_xaxis(title='Vertical effective stress [kPa]',
    panel_no=2)
resultplot.set_zaxis(title='z [m]')
resultplot.show()        

We can see that results are identical between the nodes and elements. The transitions between unit weights at layer interfaces are now also sharp because an element is either fully above or fully below a layer interface.

Currently, groundhog leaves it up to the user to decide whether to operate on nodes or elements.

When a non-zero vertical effective stress is necessary (e.g. when a surcharge is present), the initial value can be specified using the start_value argument. Here, we can work with a surcharge of 40kPa:

grid.elements['Vertical effective stress [kPa]'] = \
    grid.elements.depth_integration(
        parameter='Effective unit weight [kN/m3]',
        outputparameter='Vertical effective stress [kPa]',
        start_value=40)        

As expected, the vertical effective stresses are shifted to the right by 40kPa.

Conclusions

Being able to apply the trapezoid rules to arrays of depth values and corresponding geotechnical parameters at those depths is very useful in geotechnical calculations. While the standard SciPy syntax can always be used, the CalculationGrid in groundhog uses the same functionality to quickly apply numerical integration to soil profiles. Using this syntax, advanced geotechnical workflows can be created with just a few lines of code.

In the next post, I will show you how the groundhog library can be used to create nice and instructive charts for grain size and plasticity results (thanks Joe Seery for the suggestion).

Although I'm using LinkedIn articles to reach the people in my network, the editor is quite clunky, especially when you're used to all the goodness that the markdown editor offers in the Jupyter Notebook ;-). I will keep writing the LinkedIn articles, but I am also making the Jupyter Notebooks on which they are based available as part of the groundhog package https://github.com/snakesonabrain/groundhog/tree/DEV/articles . Check the DEV branch via the link above for the latest articles.