Subsurface Modeling with GemPy

6 min readOct 20, 2021

Creating Soil Profiles and 3D models with Python

Determining what’s under our feet and effectively mapping it out is is one of the first and most important parts when designing any kind of engineered structure. It comes with a significant margin of error as even the most advanceed subsurface investigation methods we use today can’t fully map out the soils, rocks, water elevations, etc, over larger areas. Whether you’re using borings (SPT or CPT), seismic refraction or some other method, you will have to end up doing a lot of interpolation to build your subsurface model.

GINT Subsurface Model using Cross-Sectional Fences

Software such as Bentley’s GINT are great for creating 2D and 3D soil profiles and their precision is really only limited by the data they’re given and the engineers experience. These subsurface models can also be created using Python, and although it’s not as straightforward as GINT , it gives us some more freedom to manipulate our models as we see fit. Once you get the hang of it you could generate good soil profiles faster than some softwares and definitely faster than plotting it in CAD.

Below is a tutorial on how to create a subsurface model using Python, with a primary focus on GemPy which is an open-source geomodeling library. It is capable of constructing complex 3D geological models of folded structures, fault networks and unconformities.

Starting our Model

Install GemPy

It can be installed using :

pip install gempy

Setting Up

Data we will need from the projects site and subsurface investigation (CSV format as we will use pandas to turn them into dataframes):

Surface_points →Site surveyed surface elevevations for the area (x,y,z)
Orientations →Assigns a formation value to surface points (Gx,Gy,Gz or azimuths). Essentially it’s angle associated with a surface point representing faults passing throught the layers
Grid →Size of your model (number of cols x rows x layers over some cubic length )
Surfaces→soil/rock layers
Series →Groups Surfaces and Faults into seperate series
Faults →Location and Information on faults being modeled (have to be inserted independently)

For this tutorial let’s calle the CSV file containing these 6 parameters as “MyProject.csv”.

The following imports are needed to model with GemPy:

Making our Geological Model

To begin we need to extract the data from the CSV files that contain our site information into a dataframe and create a model with GemPy.

This will generate our model with our data in 100 rows x 100 columns x 100 layers extending across a 2 km X 2 km X 2 km volume which is our “resolution” and “Grid size” respectively. The larger you make these, the longer it’ll take to run the model. The model takes in both the Surface_point data (path_i) and the Orientation (path_o) which it will use for interpolation and generate our 2D and 3D plots.

Modeling

Using .surface() we see that the geo_model that was made with “MyProject” data has 5 surfaces as such (basement is automatically generated,and acts as a the lower bound)

geo_model.surfaces

The surface layers need to be split into into two seperate categorical series, straigraphical surface layers (Strat_Series) and faults (Fault_Series). GemPy is a bit finicky when it comes to the order these surfaces are in. Faults always need to be moved to the top of the pile, are independent of the other, and their order determines which one cuts through the surfaces first (tectonic relationship). The Strat_Series surface should be in an order representative of the real world location we are modeling, with basement alkways being at the botttom. We can accomplish that with the following:

gp.map_stack_to_surfaces(geo_model,
                         {"Fault_Series": 'Main_Fault',
                          "Strat_Series": ('Sandstone_2',               
                          'Siltstone', 'Shale', 'Sandstone_1',  
                          'basement')},
                          remove_unused_series=True)
geo_model.surfaces

All we did was put the one Fault_Series on top but now the models surfaces are in order and Gempy should handle them without issue.

Graphing initial data

We can create 2D and 3D plots of our initial data for visualization which also shows us the interfaces and orientations.

# 2-Dgp.plotta(geo_data, 
direction='y') 
  # 3-Dgp.plot_3d(geo_model, 
plotter_type='basic')

Interpolation

The above plots don’t look like much of anything especially in a 3-dimensial space. We want to fill up this grid up and essentially make it a solid 2km x 2km x 2 km block almost as if you’d have ripped out a giant block of earth to examine it. To get this we will populate the plot by interpolating the data we currently.

GemPy interpolates using a method called Kriging, also known as Gaussian process regression, which is a method governed by covariance that gives the best linear unbiased prediction at unsampled locations. We interpolate the data as such:

interpol_data = gp.InterpolatorData(geo_data, output='geology',
                                    compile_theano=True,               
                                    theano_optimizer='fast_compile')

The kriging parameters can be changed as needed but personally I wouldn’t mess around with them.

gp.get_data(geo_model, 'kriging')

Model Computations

After interpolating we can use GemPy to create our model using

gp.compute_model(interpol_data)

The resulting model will results in two arrays being generated, one for the geological surfaces and one for the faults that cuts through them.

interpol_data.geo_data_res.get_formations()

Surfaces Array: [Main_Fault, Sandstone_2, Siltstone, Shale, Sandstone_1]

Fault Array:(5, object) [Main_Fault, Sandstone_2, Siltstone, Shale, Sandstone_1]

These two arrays also have 2 subarrays as the block model solutions.

Plotting

Using Matplotlib we can plot the model we just generated and create the soil profile we’ve been after. The following code essentially will take out 2km³ block (represented by a 100 row x 100 column x 100 layer grid) and slice it at the y=50 mark giving us the cross section from the middle of our model.

gp.plotting.plot_section(geo_data,
                         gp.compute_model(interpol_data)[0], 
                         cell_number=50, direction='y',           
                         plot_data=False)

3D Plots

We can visualize our model in 3d by using VTK which is an open-source software system for 3D computer graphics, image processing and scientific visualization. It’s already incormporated into GemPy so no need to import anything new.

We first need to get the vertices and triangulations from our surfaces and fault with the function get_surface . We can then use these to plot our 3-D model:

ver, sim = gp.get_surfaces(interpol_data, original_scale=True)

gp.plotting.plot_surfaces_3D(geo_data, ver, sim,  
                             plotter_type='basic')

Getting Fancy with some Topography

Adding topography always make a model look better than just some block. The topography can be highly customizable as you can use a real world raster image, create your own surface overlay, color change based on elevation, plenty of choices. We’ll just add topographty based on the depth and layer color.

This basic walkthrough was to create a basic geological 3D model from which you can extract 2D cross sectional soil profiles and use them for design. It definitely needs to be built up and incorporate things such as groundwater elevation, soil/rock parameters for each of the layers to truly become useful. In the end you will wind up with a very interactive model that you can easily extract information from.

References