Building Pit The Hague

Requirements

  • Be sure QGIS version 3.28.00 or higher is installed.
  • Be sure the gistim Python package is installed (see installation for instructions).
  • Download the tutorial material. Follow this link.
  • Installation of the QGIS-Tim plugin and the MOD plugin is part of this Tutorial. The necessary ZIP files are included in the tutorial material.
  • Internet connections is optional during this Tutorial. It is only required for installation of additional plugins and the use of an online topographic background map.

Description

In this tutorial, you will learn how to:

  • install and use the QGIS-Tim plugin;
  • use the basic of QGIS for pre- and postprocessing of Tim;
  • create several steady state models (TimML) and a transient model (TTim);
  • analyse the results;
  • export your model to a Python script.

Objective

Calculation of a pumping well extraction.

Introduction case The Hague

In the old city centre of The Hague a building plan with a parking basement is constructed. The dimensions of the building pit are 34.2 m * 83.2 m. The area is shown in Figure 1.

Figure 1: Schematization of the building area in The Hague city centre (orange: building pit)

To carry out the construction in dry conditions, the groundwater level needs to be lowered by 3.5 m. The required drainage time is 4 months.
The environment is vulnerable with valuable buildings and monuments (see Figure 2) that have a shallow foundation. The subsoil consists of sandy layers which locally contain shallow peat/clay/silt layers and deeper layers of clay/silt. The builder contractor must use a sheet pile wall that needs to be installed in a water-retarding layer (aquitard).

Figure 2: The monumental environment

The questions to be elaborated in the hydrogeological advice is:
“How deep should a sheet pile wall be placed to prevent damage? Can risks arise from a leak?”

Applications of a sheet pile wall (with a hydraulic resistance of 100 days) and a concrete cut off wall (with a hydraulic resistance of 1000 days) need to be studied.

The soil profile has been investigated with borings and CPTs. An example of a CPT is shown in Figure 3.

Figure 3: CPT to a depth of NAP -20 m

The interpretation of the CPT and borings gives the following soil layer scheme:

  • The soil mainly consists of fine sand in Holocene layers (max depth NAP -15 m) and coarse sand in deeper Pleistocene layers.
  • Just below the groundwater level (NAP -0.5 m), there is an approximately 1.5 m thick peat/clay/silt layer (not in CPT).
  • The layer between NAP –5.5 m and –8.5 m consists of clay lenses. There is little certainty about the resistance of that layer.
  • Between NAP –13 m and –15 m there is a clay/silt layer.

Based on experience in this region, the following geohydrological schematization is known (Table 1). The layer depth is indicated relative to the groundwater level of NAP-0.5m, because only that part is relevant for our model. In Table 1 c = vertical hydraulic resistance, k= horizontal permeability, S’= specific storage: :

Table 1: Geohydrological model of the subsoil
Cluster Layer type Layer Bottom c k S’
[m-REF] [days] [m/day] [1/m]
1 Aquitard -1 3000 0.01
Aquifer -5 10 0.00025
2 Aquitard -8 40 0.0002
Aquifer -13 10 0.00005
3 Aquitard -15 100 0.001
Aquifer -60 25 0.00001

According to DINO loket the average groundwater level on site is NAP -0.5 m. There are strict demands from authorities for allowable groundwater lowering. The maximum allowed drawdown outside the construction pit is 0.1 m.
To perform the calculation, dewatering is schematized to 8 wells. Tip: Don’t place the wells in the model too close to the sheet piles.

Tasks:

  • Which element from the Tim package would you use to place the wells and sheet piles?
  • Do you need to install sheet piles deep (in sand layers of aquifer 1 and 2) or shallow (only in the first sand layer of aquifer 1)?
  • Investigate the effect of the presence of a lower or higher hydraulic resistance in the layer between NAP -5.5 m and -8.5 m.
  • Estimate the extraction rate to keep the building pit dry.
  • What are the risks? Look at the groundwater drawdown outside the construction pit in situations with different sheet pile depth.
  • Also check for point leaks. Model these leaks with a hole in the sheet pile in the southwest corner of the area.

Getting Started

  1. Launch QGIS from your START menu, from your desktop or click on …\QGIS3.28.0\bin\qgis-bin.exe.

Intermezzo: QGIS language settings

Perhaps your QGIS was installed in another language than English. Because the Tutorial refers to the English version, let’s change to English.

  1. From the main menu click on Settings and select Options (e.g. in Dutch Extra and Opties).
  2. In the new window go to the General section (Dutch: Algemeen) on the left.
  3. Check the box to allow Override System Locale (Dutch: Landinstellingen negeren) and expand this sub menu.
  4. From the drop-down menu “User interface translation” (Dutch: Vertaling gebruikers-interface) select American English and click OK.
  5. Close QGIS and open it again to activate your language change.

We start with the creation of a new QGIS project.

  1. From the main menu click on Project and select New.

The case in this tutorial is located in The Netherlands, so next we select the appropriate projection.

  1. From the main menu click on Project and select Properties.
  2. In the Properties window select the category CRS, search for “EPSG:28992” and you find “Amersfoort / RD New”. Select this option and click the Apply button, followed by the OK button to close the window.

In case your work is mostly in The Netherlands and in the “Amersfoort / RD New” projection, consider making this your default projection.

  1. From the main menu click on Settings and select Options….
  2. In the section CRS and Transforms select CRS (handling), pick the radio button Use a default CRS and select “EPSG:28992 -Amersfoort / RD New”.
  3. Click OK.
  4. Close this window.

Install plugins

This is the moment to import four plugins needed for this tutorial:

  • the QGIS-Tim plugin (developed by Deltares).
  • the iMOD plugin (developed by Deltares).
  • the Value Tool. A plugin to display in table or plot values from raster layers (or mesh layers) at the current mouse position.
  • the PDOK plugin. This plugin gives access to a large database from which we will load the topographic maps and use the navigation option.

To install the plugins from QGIS you need an internet connection!

No internet connection? Follow the next steps to import the two Deltares’ plugins from a ZIP file provided in the Tutorial Dataset.

  1. Go to Plugins from the main menu and select Manage and Install Plugins… to open the plugin window.
  2. On the left section select Install from ZIP.
  3. Click the Browse button () and select the ZIP file “QGIS-Tim_Tutorial\QGIS-iMOD-plugin.zip”.
  4. Click Install Plugin.
  5. In the same way, install the QGIS-Tim plugin using the ZIP file “QGIS-Tim_Tutorial\QGIS-Tim-plugin.zip”.

If you have an internet connection, all four plugins can be downloaded within QGIS.

  1. At the top, find the Plugins menu (~sixth object in the menubar).
  2. Find "Manage and Install plugins" (~first object in drop-down).
  3. Find "All" (~first in left section).
  4. Search for "Qgis-Tim".
  5. Click "Install Plugin".
  6. Search for “iMOD” and install it.
  7. Search for “Value Tool” and install it.
  8. Search for “PDOK services plugin” and install it.
  9. Make sure that under Plugins > Manage and Install Plugins > Installed now the 4 added plugins are checked.
  10. Close the Plugins window.

See in the toolbar section of QGIS that the plugins are installed:

  • iMOD Toolbar
  • QGIS-Tim
  • Value Tool
  • PDOK Services Plugin

Further in this Tutorial we will use some default toolbars that might be hidden at the moment. Let’s check that and unhide if necessary.

  1. Select View from the main menu and choose Panels and be sure these two toolbars are checked:
    - Layers;
    - Browser.
  2. Select View from the main menu and choose Toolbars and be sure these three toolbars are checked:
    - Advanced Digitizing Toolbar;
    - Snapping Toolbar;
    - Attributes Toolbar.

Prepare your project

For navigation purposes, let’s load a topographic map for The Netherlands from the online PDOK database.

No internet connection? Follow the next steps to import a simple PNG file as a background.

  • Go to Layer in the main menu, go to Add layer and select Add Raster layer.
  • Use the Browse button () and from the tutorial material select “…\QGIS-Tim_Tutorial-TheHague\TopographicMapTheHague.png”.
  • Click on Add and Close the window.
  • If you do not see the map, select the layer “TopographicMapTheHague”, click your right mouse button and select “Zoom to Layer(s)”.
  • Continue after step 23.
  1. If you do have an internet connection click on the PDOK plugin button () to open the “PDOK Services Plugin” window.
  2. From the tab PDOK Services search for “pastel” and you will find a WMTS type layer called “BRT Achtergrondkaart WMTS”.
  3. Select the layer.
  4. In the section “laag toevoegen” click the button Onder.
  5. Close the PDOK window.

Our project area is in the centre of the city of The Hague so let’s navigate to that city using the PDOK plugin.

  1. Type “Parkstraat” in the PDOK search field, near the PDOK button ()
  2. One of the locations PDOK will find is “Parkstraat, ’s-Gravenhage”. Click on it and QGIS will fly you to the project location.

Let’s now open a shape file containing the circumference of the building location.

  1. Go to Layer in the main menu, go to Add layer and select Add Vector layer.
  2. Use the Browse button () and from the tutorial material select “…\QGIS-Tim_Tutorial-TheHague\building_pit.shp”.
  3. Click on Add and Close the window.

Tip: a fast alternative for adding layers: from the menu View > Toolbar add the Manage Layers Toolbar and use the button .

  1. In the Layers panel on the left, select the layer “building_pit”.
  2. Click your right mouse button and from the menu select Properties.
  3. In the new window go to the section Symbology on the left and try to pick a polygon style with only a contour color.
  4. Click on OK a to save and close the window.

Let’s save this project to be able to return to it later or in case of a crash of QGIS.

  1. Go to Project in the main menu, select Save As and select a folder and a file name for your project, e.g. “…\QGIS-Tim_Tutorial-TheHague\TheHague.qgz”

Open the QGIS-Tim panel

Now we are ready to activate the QGIS-Tim plugin.

  1. Click on the QGIS-Tim plugin button () and the QGIS-Tim panel appears.
  2. Go to the tab Model Manager.
    Here we will create an empty database (geopackage) to store all elements and parameters for the model.
  3. Click the New button to create the GeoPackage and save it for instance in the folder with your tutorial data, e.g. “..\QGIS-Tim_Tutorial-TheHague\case-TheHague.gpkg”.

Your window looks like in Figure 4.

Figure 4: QGIS-Tim panel
  1. Check in the Layers panel on the left that your new geopackage is added as a group.
    A sub group timml for the steady state model input, the sub group ttim for the transient model input and a series of output formats (vector/mesh/raster).

If you had no introduction to the Tim plugin, read the Intermezzo below for a general explanation of the components.

Intermezzo: introduction Tabs on the Tim panel

  • Model Manager: an overview of the elements in your geopackage. In case you switch to transient modelling, an extra column with ttim elements is added.
  • Elements: a list of at least 16 Tim elements from which you can build your model.
  • Results: here you can define your domain and cell size, decide if your model is transient or not and manage the output files.

Let’s save this project to be able to return to it later or in case of a crash of QGIS.

  1. Go to Project in the main menu, select Save As and select a folder and a file name for your project, e.g. “…\QGIS-Tim_Tutorial-TheHague\TheHague.qgz”

Start your Tim model

To get an idea of the parameterization of our Aquifer we can use the layer information available from the LHM database (Landelijk hydrologisch model). That model provides 3D geological layering within The Netherlands and the corresponding geological parameterization (permeability, resistance, thickness).

An example of parameter values near our building pit area is presented in Table 2. In each LHM layer the Aquitard information is positioned on top of Aquifer.

Table 2: The characteristics of the subsoil based on LHM data (example, your values may differ)
fid
-		 layer
-		 aquitard
c		 aquitard
npor		 aquitard
s		 aquifer
k		 aquifer
npor		 aquifer
s		 aquifer
top		 aquifer
bottom		 semiconf
top		 semiconf
head
0 0
1 1 287 62.37 -18.97 -38.30
2 2 1 37.94 -38.30 -60.54
3 3 3056 13.21 -66.62 -77.50
4 4 708 9.25 -84.59 -88.38
5 5 3143 10.94 -106.06 -120.63
6 6 8092 5.68 -151.11 -254.94
7 7 17752 3.10 -307.82 -390.71

A choice should be made at which depth the hydrogeological base is chosen for this project. Where high aquifer k is followed by a large aquifer c the base is chosen. For this case the base is chosen at -60.5 m (all depth data in Table 2 are m NAP). As you can see, values for Holocene layers are not available here and thus not displayed in LHM results.

Aquitard resistance in the city is quite high, so a value of 3000 days is best guess from experience.

For marine deposits at -12 a value of 50 to 100 days per meter layer thickness is a proper choice. The permeability of fine Holocene sand is chosen at 10 m/d.

We are now ready to define our first steady state model by parameterizing our Aquifer.

  1. From the Layers panel, select the layer timml Aquifer:Aquifer.
  2. Click your right mouse button and from the menu select Open Attribute table (). Search for the same icon somewhere on the Attributes Toolbar.
  3. Switch to the editing mode in the table with a click on () Toggle Editing.
  4. Then by clicking 3 times on () Add feature you will see new rows are added.
  5. Now you are able to fill in the desired values shown in Figure 5.
    NB! Do not fill in column FID. ‘Autogenerate’ will take care.
Figure 5: The characteristics of the aquifer

NB! Closing this table window with the will not save the filled in data or stop the editing mode!

  1. Click on the Save-Edits button () to save you data during the process or click on the Toggle Editing Mode button () to stop edditing and QGIS askes you if your changes should be saved.
  2. Now you can close the Attribute Table window.

In the QGIS-Tim Window we now introduce the following elements for ground water modelling in the QGIS-Tim tab Element:

  • Leaky Line Doublet
  • Well

Adding a Leaky Line Doublet

  1. In the QGIS-Tim window go to tab Elements and select the button Leaky Line Doublet.
  2. Fill in the name of the layer in the pop-up panel, e.g. ‘sheet_pile’, and click OK and this new Layer is added to the Layers panel.

First, we will draw the location of the sheet pile before adding the parameter values.

  1. In the Layers panel right click on the layer timml Leaky Line Doublet:sheet_pile and start the editing mode by clicking the Toggle Editing button ().

Two important remarks before drawing the sheet pile:

  • Start drawing at the most southern corner of the building pit to make a variation for an extra purpose later on in this tutorial.
  • In order to make a closed area surrounded with the sheet piles, the first and last point must have the same coordinate and intersect. For this we must activate the Snapping option in QGIS.
  1. In the main menu go to Project and select Snapping Options...
  2. From the new toolbar (see Figure 6), click on the Enable Snapping button ().
Figure 6: Switching on snapping options
  1. Then select the Add line feature icon () to start the drawing of the sheet pile.
  2. Draw the sheet pile around the building pit (start south corner) with your left mouse button and close the feature on the first point.
  3. In the pop-up window fill in the resistance (1000 d, see * below) and layer number (layer=0, see ** below) as in Figure 7 and click OK to close the window.
Figure 7: The characteristics of the Leaky line doublet

Remarks to the provided values:
* Because of a known issue in the TimML software, until this issue is resolved, the value of the resistance must be multiplied by the permeability. For normal sheet pile walls a resistance value of 100-250 d for interlock leakage is a good choice. In this case the chosen value needs to be increased by multiplying with the applicable aquifer permeability in the effected layer (R=100 d becomes R=10*100 d).
** In the real-world counting starts with 1. However, Tim is programmed in Python and in Python counting starts with 0. You will get used to it.

  1. Close the editing mode with a click on () and you are asked to save your changes.
  2. Open the attribute table for Leaky Line Doublet (click or press F6) to check hydraulic resistance values for the sheet pile walls.

Adding a Well

  1. In GIS-Tim go to the tab Elements and select the Well element.
  2. A name for the element can be given in the pop-up panel, e.g. pumping_wells.
  3. In the Layers panel right click on the layer timml Well:pumping_wells and start the editing mode by clicking the Toggle Editing button ().
  4. Next click on Add point feature ().
  5. With you left mouse button add the first well location similar to point 1 in Figure 8.
  6. In the pop-up window, fill in the well parameters discharge (50 m3/d), radius (0.1 m), resistance (1 d) and layer (0).
  7. Now add the other 7 wells locations in a fast way: do not import parameters with every single point. We will do that later. Just click OK on each Feature Attribute window.
Figure 8: Schematization of the pumping wells, sheet piles and observation locations
  1. In the Layers panel right select layer timml Well:pumping_wells and open the Attribute Table (F6).
  2. Start the editing mode and fill in the values shown in Figure 9. Don’t forget the column ‘Label’.
    The discharge value in the table is a first guess to be adjusted later to desired level of -3.5 m groundwater lowering in the building pit.
  3. Stop the editing mode, save your work and close the window.
  4. Select layer timml Well:pumping_wells again, click right and from the menu select Show labels.
Figure 9: Well properties

When you finished the input of the wells and their parameters, you have to look back to the combination of wells and sheet pile wall. The reason is that near a well the flow is strong and therefore the flow through the wall is increased. This may lead to extremities in the calculation if wall segments are chosen too large. Tim uses control points but they are set at regular distances on the doublet (the number of control points is 1+order). It is necessary to divide the wall in smaller sections near well locations. The length of each section should be chosen equal to approximately the distance of the wall to the nearest well.

  1. Select the Leaky Line Doublet in the Layers panel and start the editing mode.
  2. Activate the Vertex Tool () and then hoover along the sheet pile line in your drawing. The line element near your pointer lights up.
  3. With the combination of Shift + double left click, you can add a new vertex to the line. Repeat that for the number of points you want to add.

Figure 10 shows the result of positioned extra points in the line.

Figure 10: Extra points for segmentation of the sheet pile element

Computing the groundwater head drawdown

  1. Zoom in or out to desired domain for which you want to see the model results.
  2. In the QGIS-Tim panel select the tab Results.
  3. Select the button Set to current extent to define the Domain.
  4. Grid spacing will follow automatically but for now make the results mesh more dense by changing “Grid spacing” to 3.00 m.
  5. In the “Output” section give the name of the file where you want to store the results.
  6. For contouring, select the check box “Auto-generate contours”.
  7. Set the increment for contouring to a proper value: 0.5, 0.25 or 0.1 as applicable for your study.
  8. Press the Compute button to have the program perform the calculations.

A black Python.exe window pops up indicating that the TIM calculation started on the background. You can ignore this window but keep it open. Of course you van minimize it. If the calculation was completed successful, you will see this echo in QGIS. .

Studying output results

After the calculation you see that the result is automatically added to the Layers panel, probably called “case-TheHague output”. Results are presented as Mesh, Raster and Vector data. Contours are saved as vector.
Although these layers / groups are checked, the data is not visible. That is because the geopackage was added last, and QGIS adds layers at the end of the list.
Let’s move the layer “pastel” to the background.

  1. Select the layer “pastel” and drag it with your left mouse button to the bottom of the list of layers.
  2. Uncheck mesh and raster to only visualize contours on the base map.
  3. Uncheck contour lines of layer 1 and 2 to get only the result for layer 0 on screen as in Figure 11.
Figure 11: Updated results of groundwater head lowering contours due to 8 wells 50 m3/d and a sheet pile wall with 100 d wall resistance.

The lowering in the building pit is shown to be around -6.20 m. To get it around -3.50 m as demanded, the well discharge should be decreased accordingly. An extra calculation with half the flow per well will be sufficient.

  1. Go to the well attribute table, adjust the flow to 60% (30 m3/d per well) and compute the model again.

Adding observations wells

The city authorities that perform quality checks on the effect of construction projects insisted on the installation of some piezometers to assure reduction of risks for surrounding old monumental structures.

  1. In the QGIS-Tim panel go to the Elements tab and select the element Observation.
  2. Give a name in the pop-up panel, e.g. “observations”.
  3. Then go to the Layers panel and select the layer timml Observation:observations.
  4. Go to the toolbar with drawing options and activate Toggle editing ().
  5. Go right to the Add point feature () and drop some piezometers in the drawing.
    We propose to use 5 points (see Figure 8): 1 point in the centre of the building pit, 1 near the lower corner outside, 1 outside near the middle of the eastern sheet pile section and just 2 near street corners.
  6. Compute the model again.

Results of calculations at the observation locations are presented as Vector data in the Layers panel.

  1. Uncheck/check layers in order to display observation results only.

By default, the label at each observation location is the calculated head for layer 0. Perhaps you see a second number near the location. For your information: this is the “location number” label belonging to the model input in layer timml Observation:observations.

We can observe that the lowering of groundwater around the building pit is quite high due to leakage of the sheet piles or leakage through the bottom clay layer with small resistance. Therefore, it is needed to improve the wall quality or the length. First we check the effect of the clay layer by studying a cross section.

Creating a cross section

  1. To make a cross section, use the iMOD plugin toolbar and click on the Cross Section widget ()
  2. From the dropdown menu on the left of this panel, select the mesh () with *_layer_0.
  3. Press Add.
  4. Start to define your cross section with a click on the button Select location.
  5. Draw your cross section by left clicking on the map. Stop drawing with a right click. You can redraw this line anytime you like.
  6. Satisfied with your line? Click the button Plot to draw this layer in the cross section. Your screen might look like Figure 12.
  7. By using the Export button you can store results from the cross section in a CSV file.
Figure 12: Cross section of groundwater heads per layer, sheet pile wall R=100d and 8 wells at 30 m3/d each in layer0

Making calculations with parameter variations or checking bandwidth

The authorities demand a drawdown effect of dewatering at a maximum of 0.10 m at surrounding buildings. This means that improvements for leakage control are needed but first we need to discover what parameter to focus on.

To check whether the wall resistance or the bottom resistance of the layer 1 (below the building pit) is more important we can make 2 variations; one with C-clay=200 d and one with R-wall=500 d. 
Of course you can change your model input, rerun the model and overwrite your model results. The next steps show you how to change the model input and save the results in separate .gpkg and .nc files.

  1. In the input group select layer timml Aquifer:…
  2. Open the Attrribute Table (F6) and change the value for “aquitard_c” in layer 1 into 200 d.
  3. In the QGIS-Tim panel go to the tab Results and change the name of the output, e.g. case-TheHague_v1.
  4. Click Compute to run variant 1.

Check in the Layers panel and see that the results are not overwritten but added to the groups, e.g. layer case-TheHague_v1-timml Observation:observations is added to the group Vector.

  1. Fill in your calculated heads at the observation locations in Table 4 or use Excel.
Table 3: Table: calculated heads [m] in observation points for 2 variants compared to the initial situations. Your values will differ.
Observation
location		 Default:
Cc=40d
Rw=100d*
your value:		 Variant 1:
Cc=200d
Rw=100d*		 your value: Variant 2:
Cc=40d
Rw=500d*		 your value:
Pb1 centre pit -3.95 … -10.44 … -4.00 …
Pb2 south corner -0.43 … -0.54 … -0.35 …
Pb3 street corner -0.35 … -0.43 … -0.31 …
Pb4 south point -0.29 … -0.32 … -0.25 …
Pb5 east wall -0.51 … -0.70 … -0.41 …

* Known issue in TimML: you have to multiply Rw by layer permeability (Rw*10).

Let’s now run Variant 2.

  1. In layer timml Aquifer:… reset the value for “aquitard_c” in layer 1 to the default of 40 d.
  2. In layer timml Leaky Line Doublet:… change the value for “resistance” into 5000 d (500x10).
  3. In the QGIS-Tim panel go to the tab Results and change the name of the output, e.g. case-TheHague_v2.
  4. Click Compute to run variant 2.
  5. Fill in your calculated heads at the observation locations in the table above.

To get the same lowering in the building pit, in the second variation the well flow might be reduced to 33% (10m3/d per well). For the sheet pile wall, increasing the wall quality or decreasing interlock leakage doesn’t make a big difference.
We can conclude that the best investment during the phase of design would be to perform extra hydrogeological research, e.g. by making more cpt’s, borings or performing a pumping test.

Sheet piles with extra depth

Suppose the best guess value of the clay layer resistance was right, then a mitigation measure for the effect of dewatering in the construction phase could be the installation of the wall to a deeper level where additional hydraulic resistance of 100d can be found at a depth of -13 m to -15 m NAP.

If we want to create extra depth of the sheet pile we will have to introduce it in a deeper layer. There are 2 options to implement it in your model:

  • Add a copy of the geometry of the sheet pile wall to the existing Leaky Line Doublet shape and assign it to layer 1.
  • We recommend to create an extra Leaky Line Doublet element. In this case it is more easy to switch on/off this additional element in your sensitivity analysis.

How to copy the sheet pile wall to an extra Leaky Line Doublet element?

  1. In the QGIS-Tim panel go to the tab Elements and add a second Leaky Line Doublet and give it a name, e.g. “sheet_pile_L1”
  2. Go to the tab Model Manager and see that the element separately is added to the list. Here is can switch this element on / off for a calculation.
  3. In the Layers panel select the new layer timml Leaky Line Doublet:sheet_pile_l1.
  4. Open its Attribute Table (F6) and start the editing mode. The table is empty.
  5. Also open the Attribute Table of the first Leaky Line Doublet and select the existing element.
  6. Click on the Copy button () in the source table to copy the selected row to the clipboard.

  1. In the target table, paste it with the Paste button () as a new layer.
  2. Assign this new sheet pile to layer=1.
  3. Stop editing and save the new element.
  4. Click Compute to start the computation again.

The results are directly visible in the contours and cross-section again.
We can conclude that the drawdown in the building pit increases with a factor of almost 2 (-7.48 m at the centre of the building pit). Therefore, we adjust the well flow to 3.95/7.48*30=15.54 m3/d per well.

  1. Implement this change in the timm Well element and recalculate the model.
Figure 13: Cross section of groundwater heads per layer, sheet pile wall R=500d in layers 0 and 1 and 8 wells at 15.25 m3/d each in layer0

We conclude that the drawdown of groundwater level around the building pit with deep sheet piles decreased significantly.

Next also alternatives with a shallow concrete cut-off wall will be calculated for a shallow and a deep wall. In that case we have R=1000 d, but note that the input in the attribute than becomes k*R=10000 due to the error in Tim.

  1. After changing the value in attribute tables of wall elements, we can compute again, switching off and on the element for the deep wall section (tab “Model Manager” on the QGIS-Tim panel).

Again extra calculation is needed to adjust well extractions for drawdown in the building pit.
Results of calculations are gathered in the following table, showing extractions and head outside the wall at South East monitoring position.

Table 4: Effect of 4 Wall alternatives on extraction rate and drawdown South East.
Wall alternative
at c1=40 d		 Wall resistance
[d]		 Groundwater extraction
[m3/d]		 Head SE monitoring
[dh in m]
Shallow sheet pile wall 100 240 -0.39
Deep sheet pile wall 100 122 -0.16
Shallow cut-off wall 1000 208 -0.30
Deep cut-off wall 1000 80 -0.04

Installation of 15 m deep sheet pile wall or cut-off wall can be elaborated in a geotechnical design.
Still, probably some decrease of interlock leakage is needed when sheet piles are chosen. Interlock sealing or maybe irrigation of water in a shallow drain pipe around the building pit could lead to approval by authorities.

Effect of a not closed wall (about 20 cm)

Authorities are afraid that leakage incidents could be harmful for surrounding monuments. The effect of leakage can be studied by creating a fictitious little opening in the wall. This can be done e.g. by changing the values of first and last coordinate of the leaky line doublet.

  1. Select the Leaky Line Doublet in the Layers panel
  2. In the editing toolbar go to the toggle editing option () and check the vertex option ().
  3. On the drawing panel select a point in the line element of the wall and right click.

The list with coordinates appears in the vertex editor at the lower left corner of the screen. There you can edit the coordinates of all points in the line. Another option is to zoom in and select a point in the line element and drag it to a new position.

When the coordinates differ a gap results, in the studied case we created a 0.2 m wide gap. It might be elaborate to perform but the result is shown in the next figure. The result is calculated by using a small grid spacing. At this created gap a large flow results in the corner of the building pit, leading to an insufficient drawdown in the building pit and a 0.1 – 0.2 m larger lowering outside the building pit at that location.

Figure 14: Calculated effect of leakage at one corner in the building pit

When the contractor wants to restore the dewatering level inside the building pit, the drawdown outside the gap will even grow further.

Using Python for sensitivity analyses with a QGIS-Tim model

QGIS-Tim offers the opportunity to export the geopackage of the created model to a Python script. This makes it possible to use the script for other ways of calculation, e.g.: - Calculation of model results in other Python environments, like Anaconda or Spyder or in a notebook. - Use in other Python oriented programs, like the Probabilistic Toolkit.

  1. If input of all elements is ready and the model has proved to run properly, go to the QGIS-Tim panel and the tab Model Manager.
  2. At the bottom press the button Convert GeoPackage to Python script.
  3. After a short period for translation in Python the explorer panel appears where you can enter the name you want to give for the python file, e.g. “case-TheHague.py” and store it in a directory you choose to save your work. The file looks like:
Figure 15: Python script

As can be seen the converted Python file start with calls (e.g. import timml) to main necessary Python packages. After that, for each timm element in the model, all data coordinates and parameter values follow. At the end of the Python file the “model.solve” command is stated after which all head values in the domain with the desired mesh density are determined and also at demanded observation points.

As can be seen the Python script for TimML is written in a very dedicated and condensed manner.

To get the Python file running in a platform like Anaconda or Spyder, extra lines should be added at wish, to get the output that the user needs.

Next this file can be used for geostatistical and scenario-analysis. There are several ways of handling this kind of study, like writing an additional Python program to perform repeated calculations and statistical analysis on results. But an easy way is to use the Probabilistic Toolkit (PTK), a platform for statistical analysis to be used together with geotechnical design programs, developed by Deltares. The PTK can be used for study of model sensitivity for variation of parameter values or reliability analysis.

The PTK can be downloaded free of charge at Probabilistic Toolkit - download.

  1. Open the Probabilistic Toolkit from your desktop () or Windows menu (Deltares folder).

The Toolkit opens at the first of 5 tabs: Model.

  1. In the section Model Type check if the dropdown menu Type is set to Internal.
  2. In the section Model Type check if the dropdown menu Language is set to Python and a the field Version appears.
  3. Select the tree points (‘…’) in the field Version and give the path where PTK can find the Python interpreter, e.g. Spyder at “C:\roelofs_fs\_software\deltaforge\pythonw.exe”.
  4. Than you copy the content of the Python file we converted from QGIS-Tim into the window “Model code”.

We handle the process in this way because we want to change some lines in the source code to get the program running in PTK.

Next the specific parameters must be selected that are expected to be probably most relevant to variations in results. In the source code used in the PTK, those parameters will not have input on a value but need to be mentioned with a name that the PTK can use for input selection in the calculations.
The parameters that seem to be important are:

  • The hydraulic resistance of the sheet pile wall Rwall.
  • The resistance c1 of the first clay layer.
  • The permeability of the sand layer k01.
  • The resistance c2 of the second clay layer at -14 m NAP.
Figure 16: Declaration of Variables to control the analysis and first guess for values.

In PTK we have to set these parameters:

  1. In the PTK panel Input click 4 times on the Add button (), to add 4 variables and name them Rshp, czba, khol and cbasis (see Figure 16)
  2. In the tab Variables give these 4 variables the values 10000, 40, 10 and 100.
  3. Return to the tab Model and assign these variable values to the parameters in the source code by copying this code block to the top of the Python source code (see Figure 17).
    Rwall = Rshp
    c0 = czba
    k0 = khol
    c2 = cbasis
  4. Finaly, replace the fixed value (NB! 2 times “Rwall”, for each sheet pile definition) in the Python code for the parameter name (see red elements in Figure 17).
Figure 17: Replace parameters in the Python code.
  1. At the end of the source code, eliminate the following lines from the converted file:

head = model.headgrid(xg=np.arange(80940.0, 81166.0, 4), yg=np.arange(455554.0, 455360.0, -4) )

  1. Also eliminate the lines:

observation_peilbuis_2 = model.head(x=80970.16497991727, y=455495.54316352855) observation_peilbuis_3 = model.head(x=81004.86605597858, y=455389.68291883526) observation_peilbuis_4 = model.head(x=81049.06450220713, y=455485.0645022071)

  1. Because we only are interested in 2 points, use the first two observation points and add these lines at the end of the Python script:

pb0 = observation_observations_0[0]
pb1 = observation_observations_1[0]

  1. In the panel Output of the PTK use the Add button () the new defined variables “pb0” and “pb1”.
  2. Check if the program works properly. Field Analysis must contain “Run model”, field Results must contain “Single run”
  3. Press the RUN button () to calculate the model.

In the tab Run model we find the results of our calculation.

  1. Check if it complies with your earlier calculations in QGIS-Tim.
Figure 18: Result of test run in PTK to check model outcome.

If results are as expected, we can step to a sensitivity analyses.

  1. In the tab Field, set the Analysis option to “Sensitivity*.
  2. Go to the tab Variables.

For each selected parameter a distribution is defined with certain limits or characteristic values. Distribution formulas can be chosen based on knowhow of the user. Best guess of the parameter value distributions is given in the next example window.

  1. Change the distributon types for each Variable in the column “Distrubution” and fill in the other values.
Figure 19: Best guess of the parameter value distributions
  1. Select the row with sheetpile resistance (Rshp) and in the panel below the distribution is shown (see Figure 20).

Parameter value distributions of silt layer resistance (czba), basic peat layer resistance (cbasis) and permeability of Holocene sand (khol) are shown in graphs below.

Figure 20: Parameter value distributions of sheetpile resistance Rshp
Figure 21: Parameter value distributions of silt layer resistance czba
Figure 22: Parameter value distributions of basic peat layer resistance cbasis
Figure 23: Parameter value distributions of permeability of Holocene sand khol

Low and high value of resistivity are set in the tab Calculation to 10% and 90%.

  1. Press RUN () to recalculate the model in the Sensitivity mode.
  2. After the run is finished, a new tab Sensitivity is visible. Go to the tab.

Here we can find what extent parameter variations contribute to model results. In Figure 24 it is shown what parameter variation means for drawdown in the building pit. We conclude that for a situation with a sheet pile wall in only the first sand layer the variation of the resistance of the loamy layer (czba) determines the drawdown, and with that this factor determines the amount of extraction in that situation for the largest part. Translated to practical considerations, it is worthwhile to spend extra budget on determining the homogeneity of that layer and the vertical permeability of the loamy layer in more detail.

Figure 24: Sensitivity of drawdown in building pit to parameter value distributions of sheetpile resistance (Rshp), silt layer resistance (czba), basic peat layer resistance (cbasis) and permeability of Holocene sand (khol).

However, with this model we get a varying outcome for the level in the building pit and this is not realistic for a practical situation. In real projects, the dewatering contractor would like to know what the effect on dewatering demand and interaction with environment would be if the aim is to reach the designed groundwater head in the building pit while adjusting the rate of extraction in the wells.

We can do this quite simply with the aid of the Python file. To perform this analysis we follow the assumption that the relation between drawdown and well extraction is linear.

  1. Go to the tab Model and in the panel Model code remove from the end of the script the lines for pb0 and the line for pb1.
  2. On the top of the script, below the line k01 = khol add a value for the fixed total extraction flux and the fixed head within the building pit:
    Qws = 100
    pb0d = -3.5
  3. Now add the following lines to the code. NB! Use the same observation variable name as is used in your script!
    pb0 = observation_observations_0[0]
    fact = pb0d/pb0
    Qtot = 8*Qws*fact
    pb1 = observation_observations_1[0]*fact
  4. Within the section of each of the 8 wells, go to the line defining the pumping rate (Qw) and make it variable: Qw=Qws
  5. Recalculate the sensitivity.
  6. After the calculation is finisched, go to the tab Sensitivity and in the top left field select “Qtot” and or “pb1”.

Again, we see that the parameter value for the hydraulic resistance of the loamy clay layer (czba) is most relevant to the variation of the outcome, not only for the drawdown effect in the area surrounding the building pit (pb1, see Figure 25) but also for the total flow to the extraction wells (Qtot, see Figure 26).

Figure 25: Sensitivity of Pb1 (drawdown outside building pit) due to parameter value distributions of sheetpile resistance Rshp, silt layer resistance czba, basic peat layer resistance cbasis and permeability of Holocene sand khol with constant lowering inside the building pit.
Figure 26: Sensitivity of Qtot (groundwater extraction from the building pit) due to parameter value distributions of sheetpile resistance Rshp, silt layer resistance czba, basic peat layer resistance cbasis and permeability of Holocene sand khol, with constant lowering inside the building pit.

Based on the same statistic distributions for the parameter values we now are also able to calculate the uncertainty in output for this case, using a sheet pile wall in the first aquifer.
The distribution of the results for drawdown just outside the building pit is given in Figure 27, which in the lower part shows the distribution for total flow rate from the dewatering in the building pit.

Figure 27: Uncertainty in drawdown outside building pit due to parameter value distributions of sheetpile resistance Rshp, silt layer resistance czba, basic peat layer resistance cbasis and permeability of Holocene sand khol with constant lowering inside the building pit
Figure 28: Uncertainty in groundwater extraction from the building pit due to parameter value distributions of sheetpile resistance Rshp, silt layer resistance czba, basic peat layer resistance cbasis and permeability of Holocene sand khol with constant lowering inside the building pit

The probabilistic model gives a 50% probability on a flow of 297 m3/d, even with a 10% chance of an overrun to 659 m3/d. The groundwater drawdown at the standpipe outside the building pit amounts -0.44 m with 50% probability but with 10% value of -1.10 m.

Moreover, it is possible to perform a calculation of reliability with the PTK.

  1. Go to the tab Analysis and select the option Reliability.
  2. Then we check the variable distributions again.
  3. Next, we open the tab Calculation.

This gives the possibility to introduce a failure criterium for certain variable, like in our case pb1. We choose a less strict criterium because with a shallow sheetpile wall it is not possible to meet the requirements anyhow.

  1. Choose a comparison with undershooting a critical value of REF -0.5 m for the groundwater head in pb1 (see Figure 29).

At the right hand side in the settings window a calculation method needs to be selected. Many statistical methods are available, like Monte Carlo method etc. They are described in the help of the PTK in the main menu. We chose FORM (First Order Reliability Method) as a fast method because it is an automated intelligent iteration process to find a design point for reliability by using first calculated realisations (p.97 in the Help manual, see PTK menu).

Figure 29: Setting up a reliability calculation by choosing a failure criterium and selecting a calculation method.
  1. Recalculate the model and go to the tab Reliability.

From the result (see Figure 30) we conclude that the FORM method is able to find a solution in just a few iterations.

Figure 30: Result of the reliability analysis for the selected failure criterium.

It follows that the probability of failure for even a not very discriminating criterium as a drawdown outside the building pit of REF -0.5 m is 42%. The calculation method also returns the contributions of all relevant parameters to the outcome. It is obvious that the loamy layer at shallow depth below the building pit does not deliver enough hydraulic resistance to isolate the building activities from surrounding monuments.

As is clear, the studied layout of the building pit plus dewatering will not meet the demands of the local authorities regarding the limitations on drawdown. By waterlocking the sheet pile wall we can hardly expect a satisfactory decrease of drawdown effects around the building site, therefore, the best solution is to install the sheet pile wall to a deeper level.

Repeating the calculations with sheetpile wall in 2 layers the results will prove the possibility to meet the demands on groundwater effects with that building pit design.

This description of the possibilities like exporting a Python script from Qgis-Tim and using this in statistical analysis is just a short peek into all the options. To find more information we must refer to the manual of the Probabilistic toolkit PTK.

Creating a transient model

At this moment, it is not possible to create a transient (time dependent) model with correct answers to case studies with sheetpile walls. The TTim solution needs to be updated.
Nevertheless in this section we want to guide you on setting up a basic transient model.

  1. Return to QGIS and in the QGIS-Tim panel go to the tab Model Manager.
  2. In the section Model Setup turn on the option “Transient”.
  3. In the panel Layers select timml Aquifer:Aquifer.

As can be seen in Figure 31, the matrix of properties is extended in the transient mode with cells for transient parameters, like specific storage coefficients (aquitard_s, and aquifer_s) and porosities (aquitard_npor and aquifer_npor). Take care that specific storage coefficient is a storage per meter layer thickness. In older MWell course material a speadsheet is presented for estimation of coefficient values but there we calculated storage coefficients per whole layer. Those values should be divided by layer thickness.

Figure 31: Attribute Table for Aquifer in transient condition.
  1. Fill in the Table with the corresponding values.
  2. In the TTim section in the Layers panel go to the list of elements for transient calculations.
  3. Select ttim Temporal Settings:Aquifer and open its Attribute Table (F6). Here the length of the time series is specified.
  4. Set the starting time to 0.

But calculations take place for the range of time steps between “time_min” and “time_max”. The time unit is according to central settings in the project, mostly time in days. “Time_min” should not be zero but given a slight offset. Here we also find a “Stehfest_M” number, related to the number of terms in the calculation method performing transformation of formulas in solving differential equations. A reference date can be set to relate the calculation to real time data.

  1. Fill in the table according to Figure 32.
Figure 32: Table for ttim Temporal Settings:Aquifer

NB!: In QGIS it is an obligation to use a “reference_date” to get the model running and present outcomes!

  1. Select ttim Computation Times:Domain and set the time values for time steps where we want to get calculations of spatially distributed heads in mesh or raster points and contours.
Figure 33: Table for ttim Computation Times:Domain

For this case we will leave this table blank, because calculating several head distributions at different time steps is rather time consuming during this course.

  1. Select ttim Observation:observations and here we can set the time steps for a time series at certain points we want to monitor the calculation results. We set points at intervals that follow a (approximately) logarithmic increase in time intervals. Proposed values are shown in the next table.
Figure 34: Table for ttim Observation:observations

Don’t choose exactly the same time as when changing the flow. So if 60 is an end value, choose 59.9.

In the Attribute Table of layer ttim Well:pumping_wells we can set the discharge amount. The value might change during time when excavation is considered or construction is performed in phases. For each change we have to add a line of input. In case there are a lot of switch on/off moments, this is a lot of work.
We just leave this option without any change because there is an other option in case each well has only a singel start and end time.

  1. Therefor go to the layer timml Well:pumping_wells layer. This layer contains simple elements to make it behave as a simple Timm element.
  2. In the Attribute Table make start_time=5, end_time=30 and discharge_transient=15.
  3. But fill in a value of 0 m3/d in column “Discharge” because this variable is part of TimML well element. We use the “discharge_transient” variable to define flow in TTim. Remember: Solution of TTim is superposed on solution of TimML.