OR/14/029 Model construction and workflow
The construction of the GSI3D model was carried out on a tile by tile basis by the geologists listed in Table 6. A Metadata diary recorded the modelling process for each individual tile, with this overall metadata summary document prepared for the combined model.
A standard GSI3D workflow (Kessler & Mathers 2004; Kessler, Mathers & Sobisch, 2009) was followed for constructing the cross-sections and geological unit distributions (outcrop and/or subcrop).
Table 6 Allocation of model construction work
|Tile||Modeller||'Start' date||Completion' 'date|
|Area' 1'||H Burke/S Mathers||2006||2008|
|Area' 2'||H Burke||2007||2008|
|Area' 3'||J Ford/H Burke||2006||2008|
|Area' 4'||S Mathers||2007||2007|
|Area' 5'||S Mathers||2007||2008|
|Area' 6'||J Ford/H Burke||2007||2008|
|Area' 7'||H Burke||2010||2011|
|Area' 8'||H Burke||2009||2011|
|Area' 9'||S Thorpe||2009||2011|
|'Area' 10||S Mathers||2009||2009|
|'Area' 11||S Mathers||2009||2009|
|'Area' 12||H Burke||2009||2010|
|'Combined' model||R Terrington, H Burke||2012||2014|
A framework of 922 cross-sections was constructed in the modelled area, spaced up to 3 km apart (Figure 7). This includes shallow ‘helper sections’, added to aid the volume calculation of particular units. Helper sections are especially needed along the length of alluvium and through polygons that fall between sections to provide extra depth constraint during calculation. On completion of a tile, docking sections were constructed along all the bounding grid lines, these were iterated with the adjacent model tiles as described above.
For guidance, the 1:50 000 scale geological map polygons were rendered to the DTM and displayed during section construction. However, where borehole evidence contradicted the mapped linework, precedence was given to the borehole. During borehole coding for the project, the borehole start height was entered when recorded on the log, or taken from the DTM in the absence of a start height. Thus, true borehole start heights were honoured wherever possible during cross-section construction.
Figure 7 Framework of cross-sections used to construct the GSI3D model
When the cross-sections for a particular tile were completed, the envelopes (coverages) of each of the geological units were constructed. DiGMapGB-50 polygons were used and/or edited to delineate the outcrop extent of the geological units, and as necessary, these were combined with the subcrop portion defined by the cross-sections and boreholes.
To create the combined model, all envelopes from each tile were exported as polygon shape files using the ArcGIS tools for GSI3D. A single shape file was then created for each geological unit using the ArcGIS merge' tool on polygons with the same model code/stratigraphy (e.g. merge all ‘alv’ polygons). Next, the dissolve function was used to remove any overlaps or internal boundaries in each unit to make continuous polygons. All cross-sections from the modelled areas were then loaded into an empty GSI3D project and the newly produced coverages were imported into the corresponding geological unit.
Several checks were carried out at this stage, such as ensuring that only one version of each cross-section was loaded into the GSI3D project. Particular attention was paid to the original area boundaries, where duplicate versions of docking sections were removed if they had been loaded from more than one tile. The distribution of each geological unit was checked to ensure that real ‘holes’ within coverage polygons had been preserved following the GIS processing. The polygon data was also checked for inconsistencies, such as duplicates and mismatches across geological map sheets.
Once calculated, the surfaces generated for the combined model were checked for artefacts, especially along tile boundaries. These inconsistencies were addressed by iterating the cross- sections in the affected area.
The 25 m DTM files used in the original model tiles were far too large for the combined model to process. Therefore, a more generalised Bald Earth DTM with a 100 m cell size was deployed for the entire model area, the outcrops of superficial deposits were fitted to this dtm. Each cross- section in the model was examined and adjusted accordingly to ensure that artificial and superficial geological unit bases correspond at crossing points and to match their envelope boundaries to cross-sections. Whilst obeying the borehole data, river terraces and alluvium were re-shaped in the cross-sections and their coverages adjusted to give geologically sensible results.
A DiGMap mismatch at the north-south oriented boundary between 1:50 000 scale map sheets 255 (Beaconsfield) and 256 (North London) was also addressed in the combined model. On the
western side of this boundary, the most recent edition of DiGMapGB-50 at the time of writing (version 7.22) shows Winter Hill Gravel Member on the Beaconsfield sheet and Westmill Gravel Member to the east, on the North London sheet, with the dividing line running along the map sheet boundary. Taking into account the survey dates of these two geological maps, precedence was given to the Beaconsfield sheet, which was surveyed more recently, and the entire polygon was modelled as the Winter Hill Gravel Member.
Tidal deposits were modelled in the south-east corner of Tile 3 in accordance with the geological map of the Thames Estuary (Sheet 272, Chatham). However, these tidal deposits continue south onto Tile 6 and westwards on either bank of the River Thames as far as Silvertown, they are not differentiated from the alluvium elsewhere in the model. So to ensure consistency across the model, tidal deposits were reassigned to alluvium in Tile 3 (see also Section 6 below).
The distribution of Thanet Formation was also revised in the combined model. The modelled subcrop of Thanet Formation was based on its mapped distribution and thickness in the BGS publication Geology of London (Ellison et al, 2004, Figure 9). Following a reassessment of borehole data for the HS2 3D model, the Thanet Formation subcrop was revised in the HS2 model and incorporated into the London Basin model. To ensure that the re-interpreted boundary of Thanet Formation matched the wider London Basin model, borehole data used in tiles 4 and 12 was re-examined in a GIS and the Thanet Formation boundary was adjusted accordingly.
Because of the modifications outlined above, the combined geological model supersedes all the individual model tiles.
The final version of the combined model is LL50k_Area1_to_12_v204.gsipr.
A regional GVS (stratigraphic sequence file) was used for all individual model tiles and the combined model for continuity. The ‘London’ GVS was based on the pre-existing Thames Gateway GVS, but was adapted to generalise the level of detail, particularly in the alluvium. The code ‘Alv’ (an abbreviation of alluvium) was used in the London GVS to define the base of all alluvium deposits, whereas the Thames Gateway model separated out individual peat horizons and intervening silt and clay layers; these are not included in the combined model.
Progressive versions of the GVS were created when new geological units were added as they were encountered with the expansion of the modelled area. On completion of the combined model, a new GVS was created that lists only the artificial and superficial geological units in the model. This revised GVS is named Areas1_12_Quaternary_GVS. Each geological unit in the GVS file is attributed with stratigraphy, lithology and age, with stratigraphy used as the primary basis for modelling.
Similarly the London Basin GLEG (colour legend) file applies to the region as a whole, and is also based on the Thames Gateway file. To match the geological map sheets, the DiGMapGB-50 colours were used in the London legend file (Figure 2). A legend file specifically for the superficial model was created, named Areas1_12_Quaternary_GLEG.
For each bedrock unit in the GSI3D model, the interpretations in the cross-sections were exported en masse as single Curve objects to a GOCAD® ASCII file (one file per unit base); these files were then loaded into GOCAD®.
Unit envelopes (coverages) were also imported into GOCAD® as Curves with a z-coordinate of zero metres.
The DTM was supplied as a GOCAD® surface exported from GSI3D: this was loaded into GOCAD®.
The overall fault pattern is shown in Figure 8. The faults were initially digitised in GSI3D and then exported and adjusted in GOCAD® to model the bedrock units. Figures 9 and 10 show more detailed views of the modelled fault network, with individual faults labelled.
Figure 8 Overview of fault pattern, with the eastern (a) and western (b) faulted areas shown
Figure 9 The eastern fault area with faults in the model labelled
Figure 10 The western fault area with faults in the model labelled
- GOCAD'®' 'MODELLING WORKFLOW Derivation of 3D subcrop information'
In order to model each surface correctly using the supplied datasets, both correlation lines along sections the extent (subcrop) polygons must be taken into account. GSI3D can export directly the base of a geological unit across all sections, so these data were straightforward to obtain. However, because the unit envelopes, or subcrops, are defined only as 2D map polygons, a procedure must be defined in order to assign z-coordinates to the subcrop data so that they can be used for 3D modelling. This process is handled automatically within GSI3D but requires a
manual implementation in GOCAD®.
The general idea is that surfaces S'i, S'i+1, …, S'N , i = 1, …N are the N unit bases that comprise the model (in GSI3D these would comprise the GVS). S0 is defined as the model capping surface and S'N is the lowest surface in the stratigraphic sequence. In other words, for all points (x, y) in the district, S'i(x, y) > S'j(x, y) if i < j (i' ≥ 0 and j > 0).
We define a set of intermediate capping surfaces C'i , i = 1, …N, where each surface C'i is the minimum of surfaces S0, …, Si-1 where they exist. By definition all points on the subcrop line of unit i lie on C'i and hence z-coordinates can be assigned by querying C'i at all (x, y) points on the subcrop line.
The problem is therefore to compute the set of capping surfaces; an implementation in GOCAD® is as follows:
- On the initial model capping surface we create N new properties Z'i , i = 1, …N, one for each unit base in the model (note that the Z property of the surface (no subscript) is the one that defines the geometry of the surface). For the London model these properties were Z01LNM,
Z02CMBS, Z03STHP, Z04WIDS, Z05SAHP, Z06SWCL, Z07BGS, Z08CLGB, Z09LC, Z10HWH, Z11LMBE, Z12TAB
- Starting at the first unit (i=1):
- The outline curve for the i’th unit base is projected onto the capping surface. This interpolates the Z property of the surface onto the Z property of the nodes of the outline curve.
- The points for the unit base are assembled from the unit’s correlation lines and the 3D outline curve that was defined in step 3.
- The unit base is modelled within the geographical extents defined by its subcrop polygons.
- The z-coordinates of the modelled surface are transferred onto property Z'i of the capping surface, where i is the number of the unit base in the sequence (i = 1, …, N).
- Using a property script on the capping surface, we set property Z equal to Z'i (e.g. Z = Z01LNM for the first horizon)
- Repeat for the next unit base in the succession (i=i+1) – go to step 3
Applying the above procedure to the horizon extent polygons, a set of 12 3D polygons was generated (one per unit base in the model). Each extent polygon was turned into a set of points and combined with the corresponding GSI3D unit base to give a single set of points that defines the known unit base.
'Area' of interest
Because some unit bases have many disconnected parts, it is impractical to model each patch separately, as would normally be done. Instead, the model was constructed over the entire area of interest and then cut by the outline curves, with the unwanted parts being discarded
'Modelling' the unfaulted surfaces
For each unit base, a set of 3D points was obtained from the correlation lines along GSI3D sections and the 3D subcrop lines obtained as described in 7.1 above. Each unit base was then modelled across the entire area of interest using the GOCAD® Structural workflow.
'Faults' in the model
It was initially hoped that fault surfaces exported from the incomplete GSI3D London bedrock model could be used unchanged in the GOCAD® model. Unfortunately, the variable quality of fault meshes in the exported surfaces led to problems with computing correctly the contacts between faults; the decision was therefore taken to re-model the faults in GOCAD® in order to get clean fault meshes.
The remodelled faults were introduced to the Structural workflow and Fault-Fault contacts were established. After checking these, Horizon-Fault contacts were set up and the fault cuts computed. The first pass of this threw up many errors that were due to points lying on the wrong side of the fault surfaces (something that will be common in older versions of GSI3D). These were corrected by a combination of exclusion by distance from the fault surface and by manual inspection (both of these operations are in the Structural workflow).
The faults that are modelled are those which have been observed or inferred in the shallow subsurface; others may exist that have not yet been seen, or incorporated into the model. Faults have been extrapolated downwards into the underlying strata, where their form and extent is not known with certainty, although they are probably steeply-dipping structures.
'Further 'tidying up
A number of artefacts were also identified with respect to the subcrop polygons, where there were occurrences of section interpretations extending outside these polygons (this is obviously
physically impossible in the general case). These were again manually excluded using the tool in the Structural workflow.
The Chalk Group, Gault and Upper Greensand (combined), Lower Greensand and Base Jurassic bedrock surfaces were constructed using SKUA-GOCAD™ 2013.1. These were constructed as unfaulted surfaces and used the following data which included the rockhead model produced from the London Lithoframe 1:50 000 scale model which gave the relative elevation for the mapped outcrop, DigMapGB 1:50 000 mapped bedrock to constrain the surfaces at outcrop, and cross-section correlation points from National Fence Diagram (GB 3D 2014). All of the surfaces apart from the Base Jurassic were also fitted to deep boreholes from the BGS Stratigraphic Surfaces database.