Three-dimensional modelling of superficial geology is challenging, particularly for complex industrial sites such as the Sellafield nuclear site, Cumbria, UK. A lithostratigraphic approach is sufficient to show distribution of key units where only a general understanding of hydrogeology is required. For many sites, however, a much higher-resolution understanding is needed to better characterize heterogeneity and identify potential pathways for contaminant transport within safety case development and remediation planning. Models that focus on often thin, laterally discontinuous, heterogeneous sedimentary units are beset with issues including (1) a lack of exposure coupled with a dearth of data gathered by low-intrusive and intrusive means, (2) structural, stratigraphic and sedimentological complexities making correlation between data points difficult and (3) uncertainties associated with data quality. However, these problems may be overcome by combining land-system-based sedimentary and glaciotectonic architectural analysis of local field analogues with the creation of a portfolio of generic unit-bounding surface geometries. We test this through the development of a high-resolution 3D geological model of Sellafield nuclear site, where exposures are sparse, but shallow site-investigation boreholes are many. The results demonstrate that land-system-based sedimentary and glaciotectonic architectural element analysis can be applied to capture substantial subsurface heterogeneity resulting from dynamic glacial processes.

Research into the subsurface behaviour of radioactivity at complex industrial sites such as the Sellafield nuclear site in the UK has shown that coarse-grained unconsolidated deposits have the potential to form preferential pathways for the transport (within groundwater) of contaminants such as radionuclides, whereas finer grained deposits act as aquitards that baffle, hinder or redirect flow (Kuras et al. 2016; Smith et al. 2020). Understanding the location, extent and 3D geometry of these deposits therefore provides an important step towards improving the conceptual understanding of contaminant transport. It underpins contaminant ‘fate and transport’ modelling and subsequent predictions of risk to human health during safety case development in response to regulatory requirements.

High-resolution ‘real’ (as accurate or representative as possible) computer-developed 3D models, as opposed to more simplified ‘conceptual’ models (e.g. Fookes 1997) can better facilitate understanding of the distribution and geometry of these sediment (soil) pathways. However, geological investigation at industrial sites is challenging owing to a general lack of exposure. Suitable subsurface data for 3D modelling are typically confined to pre-construction geotechnical borehole logs, trial pits and sparse geophysical survey data. Although such data may allow the overall extent of a body of sediment to be defined, they represent only a very small percentage of that body. The building of a high-resolution 3D model that is geologically robust, and can communicate the main sources of heterogeneity, whilst understanding the degree of uncertainty (i.e. risk) therefore presents challenges, including: (1) achieving better computational representations of geological information whilst highlighting the differences between digital representation and the real phenomena; (2) responding to uneven spatial distribution of data, whilst allowing for variations in data quality; (3) creating simulated, digital versions of geological phenomena that are indistinguishable from their real counterparts (Wycisk et al. 2009).

Three-dimensional geological models where borehole logs are the primary data source can be created in several different ways, but most rely on connecting units that possess similar features. For example, the modelling of mineral ore bodies has for many years involved interpolation via a ‘kriging’ algorithm to create a ‘volume’ bounded by a 3D envelope that defines the geometry of the ore body (Krige 1951). More recent refinements have involved the use of geometric reconstruction and finite-difference modelling (Feltrin et al. 2009), and incorporation of geostatistical variogram analysis and fractal modelling (Wang et al. 2013). Inevitably, correlation involves some type of representation of ‘bulk lithology’, but the key to this is how it is interpreted and generalized. In many pseudo-lithostratigraphic models, it is often this generalization that removes the understanding of heterogeneity.

Models created using borehole data are invariably based on the correlation between boreholes of the tops and bottoms of lithologically or stratigraphically similar units, typically resulting in a grid of cross-sections (Kessler et al. 2009; Lelliott et al. 2009; Wycisk et al. 2009) and the creation of a ‘modelling stratigraphy’. Depending on the level of detail, this hierarchical sequence of units, which is used to construct the 3D model, may reciprocate known regional ‘lithostratigraphic’ sequences (e.g. Lee 2018). However, pure ‘layer-cake’ lithostratigraphy (Salvador 1994) is limited as a stratigraphic technique for glacial successions because units are commonly heterogeneous, disparate and fragmented (Rose and Menzies 1996; Lee 2018), resulting from a variety of tectonic, sedimentary, gravitational and other processes (e.g. Phillips et al. 2007; Lee et al. 2017; Lee 2018). In effect, lithostratigraphy relies upon lateral continuity of lithologies and simple rules relating to way-up, facies replication (or lack of it) and superposition, but Quaternary environments (and especially glacial and constrained fluvial settings) lack this owing to the restricted accommodation space within the landscape and presence of geological processes that were active at multiple levels within that space.

Units may also be spatially limited for purely depositional reasons: for example, channels within a braided river system change course, resulting in repeated, stacked, amalgamated fluvial sequences (Medici et al. 2015). Furthermore, in formerly glaciated regions such as west Cumbria, many units within glacigenic sequences have been affected by subglacial and glaciotectonic processes at various scales, resulting in localized thickening or removal, structural repetition, deformation, post-depositional subglacial injection of sediment by hydrofracturing, shrink-and-swell architectures and the development of widespread, gently undulating, sharp bounding surfaces (Nirex 1997a, b, c, d, e; Busby and Merritt 1999; Merritt and Auton 2000; Williams et al. 2001; Smith and Merritt 2008; Cross et al. 2018; Smith et al. 2020).

A different approach to pure lithostratigraphy is therefore required to satisfactorily model a higher level of heterogeneity, particularly regarding granular aquicludes. We describe a case study here in which these issues are addressed by utilizing and developing ‘architectural element analysis’ (Miall 1985; NACSN 2005; Slomka and Eyles 2013) and ‘glaciotectonic architectural analysis’ (sensu Pedersen 2014), together with the assessment of structural style related to deposition environments (sensu McCarroll and Rijsdijk 2003).

Our procedure first involves the definition of a sedimentary and glaciotectonic architectural element framework based on prior knowledge of the Quaternary sequence and placed within a land-system classification (see Cross et al. 2018; Evans and Benn 2021). A new portfolio of generic geometries for erosional bounding surfaces of architectural elements or lithofacies is then developed, which is used to aid confident correlation between boreholes. Our workflow comprises seven phases:

  1. literature review of site and regional geology (geological background);

  2. analysis of sedimentary and glaciotectonic architectural elements, facies associations, palaeo-environmental reconstruction and definition of generic 3D geometries of erosional-bounding surfaces;

  3. assessment of appropriate level of heterogeneity for high-resolution 3D modelling;

  4. manipulation (rationalization and homogenization) of borehole logs;

  5. cross-section construction using geological correlation of specific units within a modelling stratigraphy;

  6. refinement of modelling stratigraphy and redraft of sections according to new modelling stratigraphy;

  7. conversion to a 3D surface model and visualization and interrogation of the model itself.

The key to this approach is a thorough understanding of the geological and geomorphological processes that created the geometries and bounding surfaces of the units to be correlated. This study focuses on available field exposures of Pleistocene sediments within 5 km of Sellafield (Fig. 1), the log of borehole BH7 from the site investigation described by Kuras et al. (2016) and Smith et al. (2020), together with representative logs from the same study, and thin sections obtained from samples from boreholes (e.g. Drigg offsite borehole C) and field exposures examined during the Nirex investigations in the 1990s.

The study area forms part of the west Cumbrian coastal plain, between the Irish Sea and foothills of the Lake District National Park (Fig. 1). As shown in Figures 1 and 2 (a representative correlation panel using sections of exposures and borehole logs throughout the study area as shown in Fig. 1), superficial, recent sediments such as blown sand are underlain by glacigenic sediments that lie mainly on sandstones of the Permo-Triassic Sherwood Sandstone Group and partially on older Paleozoic rocks including limestones (Barnes et al. 1994; Jones and Ambrose 1994; Akhurst et al. 1997; Stone et al. 2010; Kuras et al. 2016; Smith et al. 2020). Nearby to the east, the foothills are formed of volcaniclastic rocks of the Borrowdale Volcanic Group and fine-grained, weakly metamorphosed, sedimentary rocks of the Skiddaw Group (Akhurst et al. 1997).

Early in the 20th century a ‘tripartite’ Pleistocene stratigraphy was recognized in western Cumbria (Eastwood et al. 1931; Trotter and Hollingworth 1932; Trotter et al. 1937). This comprised a ‘Lower Boulder Clay’, ‘Middle Sands and Gravels’ and an ‘Upper Boulder Clay’, although additional units were known to occur locally and the two tills (for a definition of till see p. 83 of Griffiths and Martin 2017) could not be distinguished where the intervening sand and gravel was absent. The lower till was linked to the last major glaciation of the region (Late Devensian), whereas tills higher in the sequence were linked to readvances of Scottish ice into the Irish Sea basin, either the Gosforth Oscillation, or the succeeding Scottish Readvance (Evans et al. 2005). Subsequent work along the coast by Huddart (1971), Huddart and Tooley (1972) and Huddart et al. (1977) broadly reconfirmed the tripartite stratigraphical model, but dismissed evidence for the earlier Gosforth Oscillation.

Following extensive investigations in western Cumbria undertaken in the 1990s on behalf of UK Nirex Ltd (Nirex 1995, 1996, 1997a, b, c, d, e, f, g, h, i, j, k, l, m; Michie and Bowden 1994; Bowden et al. 1998), a more comprehensive lithostratigraphical framework was established (Nirex 1997l; Merritt and Auton 2000). This scheme (Fig. 3), which has been further developed by Smith and Merritt (2008) and McMillan et al. (2011), includes three subgroups, 13 formations and 46 members. The Central Cumbria Glacigenic Subgroup includes erratics derived predominantly from the Lake District and Shap Fell, whereas the Irish Sea Coast Glacigenic Subgroup (ISCD) includes material from southern Scotland, the Solway Lowlands, the west Cumbrian Coalfield and the northern Cumbrian mountains. The deposits of the latter subgroup were laid down by ice flowing through the Solway Lowlands, swinging around the northwestern side of the Lake District, and feeding into the Irish Sea Ice Stream (Livingstone et al. 2012). Formations and members within the two subgroups interdigitate locally, notably beneath Lower Wasdale, where they reveal former interactions between locally sourced and far-travelled ice. An increase in the proportions of Scottish erratics in tills occurs towards both the coast and the top of several multi-till sequences in the district, indicating that Scottish ice flowing into the Irish Sea basin gradually became dominant over ice emanating from the mountains to the west (Eastwood et al. 1931; Trotter et al. 1937; Nirex 1997g).

Merritt and Auton (2000) concluded that the products of at least one major pre-Devensian glaciation are represented in cored boreholes around Drigg and Lower Wasdale. Ice possibly flowed out of Wasdale during Marine Isotope Stage (MIS) 4 to lay down the Maudsyke Till. Varved glaciolacustrine deposits (Carleton Silt) were subsequently deposited in proglacial lakes when the glacier retreated, possibly on two occasions. The Carleton Silt passes upwards into cold-water marine deposits (Glannoventia Formation) recording a marine transgression and contemporary sea-level of c. 20 m below Ordnance Datum (OD), probably during MIS 3. Glaciers emerged once again from the western valleys of the Lake District, particularly Wasdale, to lay down the very stony Holmrook Till. This occurred during the build-up of ice towards the Last Glacial Maximum (LGM), which probably reached its maximum extent between 27 and 24 ka (Clark et al. 2012; Livingstone et al. 2012; Chiverrell et al. 2013, 2018; Scourse et al. 2021), when the whole region was ice-covered. Scottish ice eventually deflected Lake District ice southwards along the coast, flowing across the valleys of the Ehen, Calder and Lower Wasdale, laying down the Low Wath and Ravenglass tills (Table 1).

The Irish Sea ice stream became dominant during several readvances recorded by thinly interbedded tills, sands and gravels within the Gosforth Glacigenic Formation. The first major readvance (Gosforth Oscillation) followed significant deglaciation. During this first readvance the Irish Sea ice stream thickened and encroached inland, and thick sequences of fine-grained, laminated glaciolacustrine sediment accumulated within Lower Wasdale (Whinneyhill Coppice Clay) and other major valleys of the district (e.g. Ehen Valley Silt). In addition, ice over-rode most of the coastal plain, much as proposed by Trotter et al. (1937). Local ice possibly expanded across Lower Wasdale to lay down the Green Croft Till. During the subsequent retreat of the Irish Sea Ice Stream, an ice marginal lake re-formed in Lower Wasdale, within which the glaciolacustrine Holmeside Clay and deltaic Mainsgate Wood Sand and Gravel were deposited. The margin of the ice stream oscillated several times, laying down the clay-dominated Drigg Beach, Fishgarth Wood and other thin tills along the coast. These minor oscillations were possibly contemporaneous with the creation of the push moraine at St Bees, widely attributed to the Scottish Readvance (Huddart and Glasser 2002, and references therein; Williams et al. 2001).

Subsequent research in the region has generally developed and established this account of events (Thomas et al. 2004; Evans et al. 2005; Roberts et al. 2007; Thomas and Chiverrell 2007; Ballantyne et al. 2009; Stone et al. 2010; Clark et al. 2012; Livingstone et al. 2012; Chiverrell et al. 2013; Cross et al. 2018; Smith et al. 2020). A radically different interpretation has been offered by Eyles and McCabe (1989), who attempted to reconcile much of the evidence described here into their controversial model of glaciomarine ice sheet retreat caused by a rapid marine ingress into a glacio-isostatically depressed Irish Sea basin. For example, they reinterpreted the widespread ‘upper till’ as a glaciomarine mud ‘drape’, inferring former relative sea-levels of up to 152 m OD. However, their local evidence has been dismissed by Huddart (1991, 1994), Huddart and Clark (1994), Nirex (1995, 1996, 1997a, b, c, d, e, f, g, h, i, j, k, l, m), and Merritt and Auton (2000). A key difference between submarine and subglacial depositional systems is the types of deformation associated with each. Large- and small-scale soft-sediment deformation features within stratified glaciomarine and associated valley infill complexes are the result of syndepositional gravity loading in response to high deposition rates and/or post-depositional slumping, whereas the compressional and extensional features in subglacial material result from glaciotectonism (McCarroll 2001). Several researchers (e.g. Harris et al. 1997; Williams et al. 2001; Phillips et al. 2007; Thomas and Chiverrell 2007) have demonstrated the unequivocal and common presence of the latter at sites around the Irish Sea Basin such as Abermawr, Dinas Dinlle, St Bees, Aberdaron, Porth Neigwl, Gretna and Killiney Bay. Glaciomarine sediments such as mud drapes have a high water content (Powell 1984) and are simply not coherent enough to have been deformed by such brittle deformation features. Eyles and McCabe (1989) did not identify evidence of glaciotectonic deformation at these sites, meaning that their sedimentary interpretation was severely compromised. In fact, their glaciomarine paradigm pertaining to the Irish Sea basin has since been categorically rejected (McCarroll 2001).

The idea of mud drapes, however, has been resurrected recently by Coleman et al. (2020), who concluded that that the massive, clay-dominated ‘upper tills’ exposed in the cliff at Drigg Beach and elsewhere in the district are glaciolacustrine clay, rather than till. Together with the intervening glaciofluvial units, they identified three lake fill–lake drainage cycles, rather than the generally accepted glacial advance and retreat scenario. They adopted this alternative idea to develop a new simplified regional lithostratigraphic framework and 3D geological model for the area around the Low-Level Waste Repository (LLWR) at Drigg. This reinterpretation has implications for future 3D modelling and site investigations in the region and is addressed further in the Discussion below.

Method for architectural element and facies analyses

Architectural element analysis allows the characterization of sedimentary geometry and facies relative to a hierarchy of erosion bounding surfaces that define the architecture of the sediments (Miall 1985, 1988, 1996, 2006a, b; Fielding 2006; Hughes 2010; Slomka and Eyles 2013, 2015). It has been shown to be very effective in the characterization of clastic fluvial sequences (Mountney and Thompson 2002; Cain and Mountney 2009; Ghazi and Mountney 2009; Banham and Mountney 2014), and has been applied to a range of glacigenic settings (Eyles et al. 1998; Boyce and Eyles 2000; Duller et al. 2008; Slomka and Eyles 2013, 2015), improving attempts to extrapolate the understanding of sedimentary geometries and properties in the subsurface, particularly in connection to hydrogeological modelling (Slomka and Eyles 2013; Medici et al. 2018). It is therefore appropriate for the high-resolution 3D geological modelling described here, as well as for modelling lower-resolution, lithostratigraphic bounding surfaces, although the effects of glaciotectonic deformation cannot be ignored (see below), for these can generate new elements (e.g. tills) or may alter pre-existing sediments in such a way as to affect hydraulic transmissibility (e.g. pervasive glaciotectonites).

Sedimentary architectural element analysis is essentially a hybrid of lithofacies analysis and allostratigraphy. The former is the interpretation of the genetic origin of sediments utilizing observations of lithology and sedimentary structures (Cain and Mountney 2009; Slomka and Eyles 2013; Medici et al. 2018), whereas allostratigraphy is the hierarchical arrangement of discontinuities or bounding surfaces that compartmentalize packages of rock or sediment (Chang 1975; NACSN 1983; Lee 2018). Allostratigraphy focuses on the geometrical arrangement of units rather than simply their lithological properties (Salvador 1994; Rawson et al. 2002) and provides a suitable stratigraphic framework to investigate detailed architectural relationships that a purely lithostratigraphic approach cannot (Lee 2018). Although it shares principles of sequence stratigraphy (Emery and Myers 1996, 2009; Miall 1997), it is not controlled by cyclical changes in sea-level or sedimentation (Slomka and Eyles 2013, 2015). The lithostratigraphy adopted here, as outlined by Nirex (1997l) and Merritt and Auton (2000) and incorporated into a larger scheme proposed by McMillan et al. (2011), is fundamentally a hybrid of litho- and allostratigraphy (McMillan and Merritt 2012).

Glaciotectonic architectural analysis is based on the identification of packages of glaciotectonically generated or deformed sediment bounded by structural features such as décollement surfaces (Pedersen and Gravesen 2009; Pedersen 2014; Pedersen and Boldreel 2015; Lee et al. 2017; Gehrmann et al. 2019). It recognizes that certain deformation styles comprising combinations of structural types are associated with specific depositional environments (McCarroll and Rijsdijk 2003), providing additional constraints on the correlation between subsurface data points in 3D geological modelling. Local and regional field analogues of glaciotectonic deformation can be used to better predict the nature and geometry of likely deformation features within glacigenic and non-glacigenic sediments underlying poorly exposed industrial sites.

The first stage of architectural element analysis is typically an assessment of all the lithofacies present within a study area (Medici et al. 2015). Considering the complex glacial sequence in west Cumbria, particularly with respect to till units that locally include other interbedded deposits, Cross et al. (2018) proposed a land-system-based, process–form model for the Sellafield area that relates glaciation style and dynamics. This approach provides a broad template for interpreting glacigenic landform–sediment or ‘facies’ associations (Evans 2017). Although, from observations made in the field, we have separated architectural elements possessing ‘ice contact’ from ‘ice marginal’ architectures (sensu Slomka and Eyles 2013; although we use the term ‘proglacial’ rather than ‘ice marginal’, as the latter must logically include ice contact architectures). The base of the ice sheet was probably underlain by a deformable bed of water-saturated sediment during the glacial readvances that affected the district (Cross et al. 2018), which would have aided rapid ice flow and the deposition of subglacial traction till and glaciotectonite (sensu Benn and Evans 1996; Evans et al. 2006; Evans 2018).

The next stage involves defining contact relationships. Miall (1985, 1988, 2006a, b) defined several orders of unit-bounding erosional discontinuities for clastic sedimentary sequences, but here we adopt a modified version of the hierarchies suggested by Boyce and Eyles (2000) and Slomka and Eyles (2013, 2015), as shown in Figure 4. For instance, we define facies associations by laterally continuous fifth-order erosional bounding surfaces that can be traced within a defined lithostratigraphic unit, thus aiding unit correlation during 3D geological modelling. Typically, as also noted by Slomka and Eyles (2013), facies associations record changes in allocyclic-controlled conditions within a depositional environment, whilst also systematically organizing complex, heterogeneous and otherwise unmappable units. Our scheme also embraces variations in surface geometry caused by the effects of glaciotectonic processes. At sites such as Sellafield, where there is little or no exposure, any closely available natural and anthropogenic exposures provide important information regarding the geometry and architecture of units to aid correlation between boreholes. Such information, as presented below, allows a more informed 3D model to be created than relying on computer algorithms alone. This has been a focus in recent years of geological modelling packages such as GSI3D (Kessler et al. 2009; Smith 2010) and Groundhog, both developed by the British Geological Survey, wherein it is up to the geologist to draft correlation lines (themselves representing surfaces) with geometries that are appropriate to the units they bound.

Architectural element analysis: description and interpretation of lithofacies and architectural elements

In the following section we present the results of architectural element analysis, providing a description and interpretation of each element, with associated field examples. Where possible, on each graphic log of a field exposure we note lithofacies, architectural element, relevant bounding surfaces according to the approach described above and palaeoflow indicator rose diagrams (data comprise azimuth and dip of clinoforms, and we use GeoRose software; Yong Technology Inc. 2014). We also identify relevant lithostratigraphic units that have been formally established for the region by McMillan et al. (2011), the abbreviations of which are shown in Table 1. Particular attention is paid to structural evidence observed at exposures, as it is crucial to fully understanding the 3D relationships of units, for accurate correlation and aiding critical examination of borehole core.

Following the results of comprehensive investigations in west Cumbria reported by Huddart (1991, 1994), Huddart and Clark (1994), Nirex (1995, 1996, 1997a, b, c, d, e, f, g, h, i, j, k), Merritt and Auton (2000), Smith and Merritt (2008) and Smith et al. (2020), most clay-rich diamicton units are judged here to have been formed subglacially; gravel, sand and silt deposits record glaciofluvial deposition, whereas laminated clay deposits formed mostly in a glaciolacustrine setting. The term ‘diamicton’ has been used following the definition of Martin et al. (2017), which states that ‘a diamicton is a granular deposit with limited fine-grained material that engineering geologists would describe as either gap-graded or part of a soil sequence containing lenses/layers of different grain-size materials’.

In total, 19 lithofacies are recognized in the study area as shown in Figure 5 and Table 2 (see also figs 5 and 7 of Smith et al. 2020). These occur in five facies associations that accommodate 11 sedimentary architectural elements (AE) and two mélange elements (Fig. 6). Deposition of sheet-like massive diamicton (facies association FA1, Fig. 6) is recorded by the occurrence of massive diamicton (matrix- and clast-supported) and laminated diamicton on bedrock sedimentary architectural elements (MD, MDc and DLbr), and pervasive and non-pervasively deformed mélange architectural elements (Mp and Mnp). Gravelly confined flow to unconfined flow deposition (FA2, Fig. 6) is recorded by presence of channel bar confined flow sedimentary architectural element (F1), unconfined flow gravel sheet element (F2), confined flow gravel clinoform element (F3), solitary sand/silt concave fill element (F4) and stacked, amalgamated sand/silt unconfined flow element (F5). Deposition of the transition from confined flow to lacustrine deposition (FA3, Fig. 6) is characterized by stacked amalgamated gravel, sand/silt element (F6) and rhythmite element (Rh). Sand-dominated, unconfined flow deposition (FA4, Fig. 6) is recorded by the occurrence of channel bar unconfined flow element (F1), unconfined flow gravel sheet element (F2) and stacked, amalgamated gravel, sand sheet element (F6), whereas high-energy channelized flow (FA5, Fig. 6) is characterized by the presence of large-scale high-discharge channel element (HDC). As shown in Figure 6, groups of facies associations record deposition in two distinct ice sheet related land systems: those displaying ice-contact architectures (FA1, FA2 and FA3), and those displaying proglacial architectures (FA1 (Mnp mélange element only), FA2, FA4 and FA5).

Channel bar concave element (F1)

Element F1 is a fine- to coarse-grained gravel with some cobbles and boulders. Clasts are subrounded to rounded and up to 70 mm in diameter. Imbrication together with planar and trough cross-bedding (Gp and Gt; Table 2 and Fig. 5) is displayed within the Peel Place Sand and Gravel Member (PPG) at Peel Place Quarry (section A–A', Fig. 7; section B–B', Fig. 8; section D–D', Fig. 9), whereas a much smaller development comprising massive, unstructured gravel (Gm) was observed within the Drigg Holme Sand Member (DHS) at Drigg Beach (location g, Fig. 10f). Element F1 typically rests on, or is interbedded with, fluvial sheet sands of element F5 as observed in BH7 (Smith et al. 2020), and in the exposures at Peel Place and Drigg Beach. Hence the upper boundary between F1 and overlying sediments is sometimes gradational (e.g. where it underlies element F5), but usually forms an abrupt second-order bounding surface. The lower bounding surface of the element is typically concave-upward and erosional, and is commonly lined by a coarse-grained basal lag as shown in Figure 9 (section D–D') and Figure 10g (location g). The surface itself typically inherits the upper boundary of the unit below, which is commonly a fourth-order surface. Units of F1 are not typically laterally extensive and may be interbedded with or isolated within other elements such as F2 and F5.

Interpretation. The stacked sets of massive and planar or trough cross-bedded pebble gravel of element F1 are interpreted to form the infill of a variety of channel types within a fluvial environment generated by repeated meltwater flood cycles (see Miall 1985). The finer part of the element records cyclic sequences of reactivation of ephemeral channels during high stage flow conditions, depositing coarse-grained bedload during waning flow and finer grained material during abandonment (Williams and Rust 1969; Slomka and Eyles 2013). Coarser material may record the formation of scours owing to channel avulsion or confluence during higher energy flow events, with deposition of coarse-grained facies and boulders within the scours (Marren 2005; Kostic and Aigner 2007; Slomka and Eyles 2013), as shown in Figures 79. F1 elements probably represent incision of stable channels owing to base-level change and/or avulsion and deposition of subaqueous bedforms within the channel (Miall 1985; Olariu and Bhattacharya 2006; Slomka and Eyles 2013).

Unconfined flow gravel sheet element (F2)

Element F2 comprises a medium- to very coarse-grained sand and gravel. It is locally massive (Gm) or horizontally laminated (Gh), but trough and planar cross-bedding (Gt and Gp) are most typical, with the geometry of foresets recording a southeastward palaeo-flow direction. The element consists of thick, laterally extensive, fine- to coarse-grained gravel sheets with a tabular geometry within the PPG in Peel Place Quarry (section A–A', Fig. 7) consisting of clast-supported pebbly gravel with some sand. The element is generally massive and structureless (lithofacies Gm as in Table 2 and Fig. 5) as displayed in the PPG at Peel Place Quarry where F2 forms several laterally extensive and laterally continuous units overlying a thick sequence of silts and sands of element F5 (section A–A', Fig. 7; section C–C', Fig. 9), and at Drigg Beach where element F2 is represented by the lower gravelly unit of the PPG (Fig. 10f–h), together with smaller developments at the base of the DHS, and stratigraphically higher in the sequence (Fig. 10b and c). Element F2 also can be seen in the cliffs at Seascale beach, where it forms a laterally continuous sheet (Fig. 11), and at Ravenglass (Fig. 12). It locally displays planar cross-bedding (lithofacies Gp) as within the PPG at Drigg Beach, where clinoforms indicate a southeastward palaeoflow direction.

Both upper and lower surfaces typically display a horizontal to subhorizontal topography and are covered unconformably by till sheets or other fluvial sheet-like elements, as at Sellafield (Fig. 13) and Drigg Beach (Fig. 10), although both surfaces may also display hummocky or non-horizontal relief.

Interpretation. Element F2 is interpreted as fluvial sheet gravel or channel floor sediment (sensu Thomas and Chiverrell 2007). Observed clinoforms are interpreted to be foresets, recording downstream aggradation and migration of gravel bars in a braided fluvial environment (see Miall 1985; Reesink and Bridge 2009).

Confined flow gravel clinoform element (F3)

Element F3 comprises stacked sets of large planar cross-bedded, fine- to coarse-grained gravel (Gp). As the exposures shown in Figures 8 and 13 illustrate, the unit is several metres thick and mostly conformable with under- and overlying units. Multi-sensor core logging (MSCL) data from borehole BH7 at Sellafield suggest that the element contains clay (Smith et al. 2020), although it is not clear whether this is within the matrix or as clay lenses. The element is sometimes continuous with F2, but with a visible lateral transition from Gp to massive, structureless Gm lithofacies typical of F2 (e.g. below location h, Fig. 10f). F3 possesses an upper surface that is typically flat and mostly horizontal, and mainly represented by an abrupt transition from coarse-grained gravel to overlying finer-grained sediment (compare the boundary between F3 and overlying F5 within the PPG at Drigg Beach; Fig. 10f–h), or an unconformable surface where it is overlain by subglacial traction till sheets. The lower bounding surface is also typically flat or horizontal, although units of F3 display onlap onto existing topography (e.g. location d, Fig. 10f), with clinoforms downlapping onto the surface of the underlying unit.

Units of F3 thus possess a similar sheet-like geometry to F2, with which they are often laterally continuous, ranging in thickness from <0.5 m to tens of metres, but with a much smaller lateral extent than F2. Both upper and lower surfaces display a topography that ranges from relatively flat and subhorizontal, when deposited on a similar topography and covered unconformably by (for example) till sheets, to hummocky. Contrary to element F2, however, there is a prevalence for the lower surface to possess a concave appearance, whereas the upper surface sometimes displays a hummocky relief, as demonstrated at Sellafield (Fig. 13) and Drigg Beach (Fig. 10).

Interpretation. The element is interpreted to represent migration of large-scale bedforms during periods of very high-energy flows (see Miall 1985, 1988; Marren 2005; Ghazi and Mountney 2009; Slomka and Eyles 2013). The presence of clay may be the result of reworking of the underlying clay unit, as at Sellafield (Smith et al. 2020). Alternatively, if it is a component of the matrix, it may indicate sudden drops in energy regime with widespread dumping of sediment. The presence of clay lenses, however, implies more quiescent phases of sedimentation within relatively still water conditions. Like element F2, sheet-like geometries of lower surfaces reflect deposition on existing low-relief topography formed by sheet-style erosion during lower energy flooding events, or during emplacement of subglacial traction till. Concave lower surfaces reflect incision during high-energy flow, whereas hummocky relief upper surfaces typically represent erosion during deposition of overlying sediments or the results of proglacial deformation. As with element F2, clinoforms are interpreted to record downstream aggradation and migration of gravel bars in a braided fluvial environment (see Miall 1985; Reesink and Bridge 2009).

Solitary sand/silt concave fill element (F4)

This element comprises mostly medium- to fine-grained sands exhibiting horizontal lamination (Sh). At Peel Place Quarry it is 1–3 m thick, several tens of metres in lateral extent, overlies elements F1, F5 and F6, and is overlain by element MD (Figs 7 and 9). Geometrically, F4 forms laterally thinning, onlapping units that infill depressions in the surface of underlying sediment. The upper surface is typically a subhorizontal erosion surface forming the base of overlying coarser material.

Interpretation. This element records a train of superimposed bedforms that represents a decrease in energy of the depositional system in response to declining discharge, leading to short-lived sand-dominated deposition in a braided river or mid-fan environment (see Miall 1977, 1985; Thomas and Chiverrell 2007; Reesink and Bridge 2009; Slomka and Eyles 2013).

Stacked, amalgamated sand/silt unconfined flow element (F5)

This element is similar in style to F4, but rather than a solitary development, it comprises laterally extensive sheets up to 20 m thick, comprising stacked horizons of sand and silt. Horizontal laminations (Sh) are common, as are planar cross-laminations (Sp). Exposures of F5 have been observed at several locations, including within the Sellafield site (Fig. 13), where it forms horizontally laminated sand units up to 5 m thick. At Peel Place Quarry it forms the lower 3–7 m of a thick (<14 m) coarsening-upward sequence of interbedded silts, sands and gravels (PPG), as shown section B–B', Figure 8, and section D–D', Figure 9. Thinner (<1 m) exposures of element F5 have been observed in the cliffs at Drigg beach within the DHS, where a thick sequence of silts and fine to medium sands are interbedded with clays of element Rh, sheet gravels (F2) and lag gravels (F1). Glaciofluvial palaeoflow indicators such as clast orientations and cross-bedding indicate SSE-directed palaeoflow (Merritt and Auton 2000). Element F5 is clearly visible in exposures observed at Seascale (Fig. 11a and b) and Ravenglass (Fig. 12a and b), where it forms thin sheets of sand overlying and interbedded with elements F1, F2 and F3, with which it forms fining-upward cycles. F5 typically underlies till sheet elements MD and MDc, and at several locations has been glaciotectonically modified to form new mélange elements Mp and Mnp (described below).

The observed thickness of <1 m to over 10 m is comparable with that observed in MSCL data from Sellafield (Smith et al. 2020). This element has a conformable, commonly gradational fourth-order bounding upper surface, which is commonly difficult to define as it has typically been created by erosion during deposition of the overlying unit. In such circumstances it is represented by either a fourth- or a sixth-order surface. The lower bounding surface of the element is typically conformable on underlying elements.

Interpretation. Element F5 is interpreted to represent fields of superimposed fluvial bedforms, deposited in an unconfined, sand-dominated braided river or mid- to lower-fan depositional environment, also as sheet flood deposits during periods of less energetic flow (see Miall 1977, 1985; Reesink and Bridge 2009; Slomka and Eyles 2013).

Stacked, amalgamated gravel, sand and silt confined flow element (F6)

This element comprises planar and sigmoidal cross-bedded (Sp, Ss) and horizontally laminated (Sh) sand, together with units of trough and planar cross-stratified pebble gravel (Gt and Gp) as shown in the exposures in Peel Place Quarry (Figs 7 and 9) and at Drigg Beach (Fig. 10c). Clay is sometimes present in the matrix and as clay stringers, leading to a variable density profile on MSCL data (see Smith et al. 2020, fig. 7). The element is a few metres thick at Sellafield, but >6 m at Peel Place Quarry, where couplets of gravel and sand are clear. It is highly likely to occur throughout the succession at Sellafield. The upper boundary is undulating, reflecting either ridge development or erosion during deposition of the overlying unit. The lower bounding surface is not definable.

Interpretation. The cross-stratification visible at outcrop and on X-ray computed tomography (XCT) radiographs suggests that element F6 was deposited in laterally migrating, incised channels. Typically located at the base of gravel–sand assemblages, where it overlies sand, the element probably represents fresh incision into older depositional surfaces after a shift in distributary location.

Rhythmite element (Rh)

Element Rh comprises thinly interlaminated clay and silt (Fl), typically forming graded couplets (Flv). Similar deposits form two thick (<30 m) sequences beneath Lower Wasdale (Merritt and Auton 2000). Laminated clays are recorded on numerous borehole logs in the region, but three thin (<1 m) units exposed in the cliffs at Drigg Beach are much more restricted in their lateral extent, infilling and onlapping onto pre-existing topography of underlying units (Fig. 10c, f and h). We recorded no evidence of lonestones within the exposed units.

Interpretation. This element represents subaqueous deposition in a low-energy, glaciolacustrine environment. A thin section of a sample of the lower sequence (Carleton Silt Formation, Fig. 3) taken from borehole QBH 20 near Hall Carleton, in Lower Wasdale (Fig. 1), revealed rhythmic, upward-fining silt–clay couplets, well-developed fine internal lamination and dropstones, confirming that the graded laminae are true, annually deposited varves (Nirex 1997k). The upper sequence (Whinneyhill Clay Member, Fig. 3) includes similar graded couplets (Merritt and Auton 2000), but other occurrences may simply represent episodic input of slightly coarser material during fluvial or slump activity into otherwise still water. Deposition mostly occurred within ice-dammed lakes that formed at several localities and times (Trotter et al. 1937; Huddart 1991; Merritt and Auton 2000; Coleman et al. 2020).

Large-scale high-discharge channel element (HDC)

Element HDC has been observed at Nethertown (Fig. 14), where it comprises a thick, chaotic body of cobbles and boulders of up to 2 m in diameter (the Townhead Boulder Gravel Member of Merritt and Auton 2000). The lower boundary of the unit is highly erosive, channelized and locally undercuts the underlying F5 element. The whole body exposed at Nethertown is over 25 m wide and a little over 6 m in height.

Interpretation. Element HDC records very high flow energies, transporting huge boulders that infilled channels created by the flow itself. Recent jökulhlaup events in Iceland recorded flows of around 4500 m3 s−1 in similar ice-proximal areas (Maizels 1997; Russell et al. 2000, 2010; Duller et al. 2008; Slomka and Eyles 2015).

Laminated diamicton on bedrock element (DLbr)

The most striking exposure of element DLbr can be seen at St Bees where it is represented by the Lowca Till Member, which can be observed directly overlying fractured bedrock of the Triassic-aged Chester Formation (St Bees Member). In the Sellafield–LLWR area the element also rests on bedrock and is represented in Blengdale and beneath Lower Wasdale (for example) by the Holmrook Till Member (HRT) of the Blengdale Glacigenic Formation, which is exposed downstream of Bleng Bridge, near Gosforth (Fig. 1). At Ravenglass (Fig. 12) and Seascale (Fig. 11c) this element is represented by the Ravenglass Till Member (RVT), which comprises a very stiff, matrix-supported, massive, stony, silty clayey diamicton (Dmm or Dml containing moderately well-dispersed, angular to well-rounded clasts up to boulder size of granophyre, granite, Borrowdale Volcanic Group lithologies and siltstones; some red sandstone, ironstone concretions and shell fragments are also present (Merritt and Auton 2000). At Ravenglass this element possesses bed-parallel and near vertical joints that strike 074–096° (Fig. 15a–c). DLbr possesses a well-developed clast microfabric, which gives it a laminated appearance.

The element is <5 m thick at Ravenglass, where it is deformed into roughly east–west-trending folds with a wavelength of 20 m and an amplitude of 3 m. It has a sixth-order upper bounding surface that has a gently undulating or hummocky relief, also observed at Seascale. Folds are bisected by shallow-angle to steep thrusts that are sometimes lined by wispy laminae of sand. At Ravenglass the lower surface is represented by a seventh-order lithostratigraphic boundary and unconformity where the element rests on bedrock. In other locations, such as in the vicinity of the LLWR site, it rests unconformably on an older sequence of tills and outwash sediments and is therefore underlain by a sixth-order surface. At both Ravenglass and Seascale the element is overlain by a conspicuous thin (5–8 cm) seam of yellow–brown, thinly laminated silt or very fine-grained sand that locally merges with wispy veins within the diamicton below (Fig. 12b and c) (Merritt and Auton 2000). At St Bees clear deformation of underlying bedrock can be observed under the Lowca Till Member, in the form of fractured and displaced blocks of sandstone. Clay can be observed infilling the fractures.

A thin section of a sample obtained from the Holmrook Till Member at Bleng Bridge is of a coarse- to very coarse-grained, very poorly sorted, immature, lithic-rich, open-packed, matrix-supported, massive diamicton, illustrating microfeatures of element DLbr (Fig. 16a), but containing no obvious sedimentary structures (Nirex 1997k). Clasts are angular to subangular in shape with low to moderate sphericity, and comprise andesite, dacite, felsite, microcrystalline granodiorite and granite. Also present are several diffuse, poorly defined or preserved pebbles of diamicton (possibly till). Clasts are locally enclosed within a rim or coating of slightly darker clay or silt, and a well-developed concentric to weakly radial, circular arrangement of finer grained clasts can be observed locally adjacent to larger lithic clasts. This arrangement is commonly accompanied by a skelsepic plasmic fabric within the clay matrix or coating. Also noted is the localized development of a weak circular arrangement of medium to coarse sand-grade clasts without any ‘core stones’ (Fig. 16a). The red–brown to orange–brown matrix is composed of fine silt and clay. At higher magnification a weakly developed lattisepic plasmic fabric, defined by short highly birefringent clay domains, has been observed. Towards the bottom of the thin section is a set of irregular fractures and voids filled with very fine-grained sand and silt. Some are filled by a microbreccia formed of small, randomly oriented fragments of the fracture walls.

Interpretation. Element DLbr possesses many of the attributes of a subglacial traction till (sensu Evans 2018). This interpretation is reinforced by the brittle deformation of underlying bedrock observed at St Bees. At Ravenglass and Seascale, initially, intra-unit deformation and intermixing would have occurred through subglacial entrainment and soft-sediment deformation as the unit was originally accreted, but the truncated, intensely sheared top of the element observed at Ravenglass and Seascale was possibly formed during a subsequent glacial advance (Merritt and Auton 2000).

The observations made of the thin section described above are fully consistent with a subglacially formed diamicton containing abundant Borrowdale Volcanic Group and Ennerdale Intrusion pebbles derived from the Lake District. The microstructures and larger structures both indicate a subglacial traction till. The microbreccia-filled fractures appear to have formed after the till was partially lithified and probably resulted from subglacial hydrofracturing (see Phillips et al. 2013a, b).

Massive diamicton element (MD)

Element MD is generally a massive, poorly sorted, clay-rich, sometimes silty, fine- to medium-grained diamicton (Dmm) (Fig. 5 and Table 2), with blocky jointing and well-dispersed, angular to subrounded granules and pebbles. Units range in thickness from c. 50 cm to >5 m, as observed at Drigg Beach (Fig. 10h), Ravenglass (Fig. 12), Seascale (Fig. 11) and Nethertown (Fig. 14). Two units of MD have also been recorded in a historical (1970s) excavation at Sellafield (Fig. 13) (Smith et al. 2020), and a probable correlative of the element is exposed at Cookson Wood, Holmrook (Fig. 1), where it includes striated blocks of sandstone (Nirex 1997g, k). It has also been penetrated by numerous boreholes, including Drigg offsite borehole C (Fig. 1), as described here.

The lowest exposed unit of MD (the Drigg Beach Till Member; DBT) at Drigg incorporates intraclasts of gravel, sand, silt and clay derived from underlying units (Merritt and Auton 2000). Few examples of macro- or meso-scale deformation were observed at any of the coastal exposures of the DBT, apart from large-wavelength folds and both compressional and extensional faults. These structures also deform the overlying sequence (Fig. 10a and b) and are thus post-depositional in nature. However, three thin sections obtained from the DBT and RVT in the 1990s provide useful microstructural information.

Thin section PY432 (Fig. 16b) was obtained from a sample of moderately brown, sandy clay diamicton (DBT) some 17 m down Drigg offsite borehole C (Nirex 1997k) (Fig. 1). The thin section is of a poorly sorted, matrix-supported, fine- to medium-grained massive diamicton containing several well-dispersed, pebble-sized (>7 mm) lithic clasts. It reveals a broadly similar clast assemblage to that observed in the Holmrook Till thin section described above, including lithic clasts with a Borrowdale Volcanic Group and an Ennerdale Granite provenance. The clay to locally fine silt-grade matrix has no obvious planar fabric(s) or discrete shears, although a weak, sporadic skelsepic plasmic fabric and clay coatings are locally developed upon or adjacent to larger lithic clasts. The sample is cut by several irregular veinlets, which are infilled with opaque material and brown coloured clay. This poorly sorted, sandy diamicton is massive with no obvious sedimentary lamination or fissility, and no obvious deformation fabrics or textures. Apart from minor brittle fractures the sample is relatively undeformed.

Thin section PY434 (Fig. 16c) was obtained from exposures of the DBT at Drigg Beach and comprises a thinly interlaminated, very compact, fine- to medium-grained sand and pebbly silty clay diamicton containing very well-dispersed fine pebbles (Nirex 1997k). The thin section shows a fine-grained, matrix-rich and weakly laminated diamicton, where the lamination is defined by a slight colour variation within the clay matrix. No obvious grading has been recognized. The sample possesses a similar clast assemblage to PY432. The matrix is very fine-grained (clay- to silt-grade) and is locally stained dark by a secondary opaque or organic phase. A variably developed lattisepic to omnisepic plasmic fabric, defined by aligned, high-birefringent (in cross-polarized light) clays, is locally developed within the matrix, although the later staining overprints these fabrics. Also visible are circular arrangements of quartz clasts, and alignment of fine silt-grade clasts adjacent to the relatively large clasts. There are hints of a weakly developed skelsepic plasmic fabric next to the latter.

X-ray computed tomography data of core from a borehole drilled within the Sellafield site exhibits additional evidence of two anastomosing microfabric crenulation fabric domains in the form of clast axis alignment, within the uppermost diamicton of a ‘tripartite’ sequence near ground level (Smith et al. 2020). The unit is matrix-supported with well-dispersed clasts of at least two distinct mineralogies. It probably correlates with the Fishgarth Wood Till Member (FWT) at Drigg Beach, which reveals more associated macro- and meso- structures than the lower diamicton (DBT) at that site. There is little apparent macro- or meso-scale intra-unit deformation, with only sporadic shallow-angle fractures (Fig. 17a–d) and imbricate stacks of shallow-angle thrusts resulting in some intra-unit thickening (locations m and t, Fig. 10g).

A thin section (TMF1) of a sample collected from a superficial, probable correlative of the FWT at Cookson Wood (Figs 1 and 16d) shows no preferred shape or length alignment of clasts but does reveal irregular clusters and a locally developed circular arrangement of sand-grade clasts without obvious ‘core stones’ (Nirex 1997g, k). A thin, opaque oxide coats lithic clasts, lines the margins of fractures, fills veinlets and stains the reddish-brown matrix, which is composed of fine silt and clay. Clasts within a higher unit of element MD (FWT) at Drigg Beach are similar in distribution and provenance to those within the lower (DBT) unit of the element. A weak lattisepic plasmic fabric has developed within the matrix of the thin section (Fig. 16d), and a variably developed skelsepic plasmic fabric occurs on, or adjacent to, larger clasts, locally associated with a concentric, circular arrangement of fine-grained sand grains around a central ‘core stone’ formed by the larger clast. Two sets of fractures may be observed: (1) oblique, irregular fractures (with limited minor displacement of material either side) that are locally offset by short, more steeply inclined fractures; (2) complex steeply inclined fractures (in this plane of section) that, in detail, are composed of an anastomosing network of finer, irregular features. A thin, finely laminated coating of fine silt or clay clearly rims voids. The clay fillings and coatings are characterized by a concentric to planar, highly birefringent fabric that is locally deformed by very fine kinks that are recognized by their sharp changes in extinction, indicating a sequence of till accretion followed by fracture, fracture infill and subsequent shearing.

Throughout all exposures, element MD is penetrated and deformed by variably common, shallow-angle to relatively steep thrusts associated with ramp anticlines and fault propagation folds, together with conjugate sets of normal faults, many of which are truncated at the base of the overlying sedimentary unit (locations b, d and f, Fig. 10f). This is especially apparent at location f in Figure 10f, where the sequence has been deformed into an eroded, NW-verging thrust-duplex sensu McClay (1992), preserved in the footwall of the large normal fault zone. Parts of the sequence have been completely detached and thrust over themselves, preserving the sequence of glaciotectonic deformation features described above, within detached and overthrust units as shown in Figure 10e (and inset 1) and Figure 17c. The duplex is not extensive, with overall individual fault displacements of 50–200 mm and an overall shortening of c. 0.5 m for the whole duplex, which, allowing for the post-depositional erosion of its top, contributes to an overall structural thickening of <0.5 m. The deformation has resulted in local structural repetition of the sequence and a vertical crenulation fabric (S2) (Fig. 17c), which appears to deform an older S1 fabric within the detached portion of DHS.

Element MD forms thin (<1 to 5 m), sheet-like units that generally become thinner higher up in the succession. They locally have an erosive base (Fig. 11c) truncating underlying strata, and commonly drape underlying topography. Although typically laterally extensive over tens, if not hundreds of metres, units may be absent or thin considerably, especially where the upper sixth-order surface is erosive. These features typically create an undulating, ‘shrink-and-swell’ architecture, whereas the topography of both upper and lower surfaces may itself have a hummocky appearance. The upper surface of exposures of element MD are locally convoluted, incorporating sediments from the overlying element into the convolutions.

Interpretation. It is important to critically examine the genesis of this element and the underlying one described in the next section, as, together, they are the subject of some debate (Trotter et al. 1937; Eyles and McCabe 1989; Merritt and Auton 2000; Smith and Merritt 2008; Cross et al. 2018; Coleman et al. 2020). The clay-rich, pebbly lithology, with an overall absence of sedimentary structures (for example, grading or lamination), suggests that subaqueous reworking or deposition is unlikely, and that deposition directly from glacier ice is more plausible. For instance, the sparse, well-dispersed clasts within a massive structure, together with an overall lack of macro- or meso-deformation structures, yet a demonstrable presence of ductile and brittle micro-structures, together with an observed sharp or abrupt décollement-like base, are indicative of subglacial deformation (Benn and Evans 1996; Phillips et al. 2002; Evans et al. 2006). Also diagnostic, are the horizonal to subhorizontal fractures, which form when bed-parallel or bed-subparallel fissility is reactivated during unloading of the till as the glacier retreats (Evans et al. 2006).

The element may represent the dilated ‘A’ horizon of a subglacial traction till (see Evans et al. 2006, 2016). However, as such horizons are typically much thinner and rarely preserved (Evans 2018), it is more likely that it represents pre-existing sediment (e.g. glaciolacustrine clay) that has been homogenized subglacially into a subglacial traction till. This is supported by the extensive investigations summarized by Merritt and Auton (2000), which revealed that the lithology of tills in the region, notably those associated with glacial readvances, are strongly influenced by the lithology of underlying sediments. For example, tills thin, become more clay-rich and include far fewer, more dispersed clasts shortly beyond where the ice has crossed over fine-grained glaciolacustrine deposits.

The sedimentological features observed in the two thin sections of DBT, particularly the brittle minor fractures and associated veinlets in Figure 16b and the aligned fabrics illustrated in Figure 16c, are typical of a subglacial traction till that has undergone pervasive deformation. The Cookson Wood thin section (Fig. 16d) also displays sedimentary features and brittle glaciotectonic features typical of a subglacial traction till.

The NNW-verging thrust duplex and associated brittle structures that deform both the FWT and underlying DHS at Drigg Beach post-date the deposition of the FWT but predate deposition of the overlying PPG. The duplex records re-deformation of subglacially deposited and deformed sediments along thrusts that appear to have propagated from pre-existing subglacial thrusts but are steeper than the underlying décollement planes they propagate from. This suggests that the sediments have drained and consolidated a little between deformation events. The structural style corresponds to submarginal, subglacial simple shear (see Hart and Boulton 1991; Hart 1999). The vergence is opposite to the overall ice palaeoflow direction, but this may reflect a northward component of easterly ice movement.

Massive, clast-supported diamicton element (MDc)

Element MDc appears to form laterally extensive sheets (<2 m) of stiff to very stiff, mostly clast-supported diamict with bed-parallel to subparallel shears giving the impression of laminations (Dcm-s). Never more than 0.3 m (and often no more than 0.05 m) in thickness, the element incorporates less compacted, dilated zones of matrix-supported, non-sheared diamicton (Dmm), resulting in a clear platy fissility. This is illustrated in exposures at Ravenglass (Fig. 15a–c) and Drigg Beach (Fig. 17a–e).

The element mostly underlies element MD in the coastal exposures, where the boundary between them typically is a definable fourth-order glaciotectonic surface. Variations in this were observed during this study. For instance, high in the sequence penetrated by BH7 (Fig. 13, this study; figs 8 and 11 of Smith et al. 2020) this surface seems to be represented by a relatively thick transition zone between the two elements, with a non-definable top or base.

At Lowside Quarter (Fig. 18a), the element is represented by a c. 30 cm thick, sheet-like seam of stiff, thinly interlaminated clay, silt and sandy silt that dips 6°/180 S (Nirex 1997g, k), generally overlain by 3 m of gravel (Gm), probably representing element F2. The seam sharply truncates climbing-ripple cross-laminae in the underlying 2 m of silty, fine-grained sand. It is formed of fissile, alternating plates of stiff, clay-rich diamict (Dcm) and less stiff, silty, matrix-supported diamict (Dmm), separated by laminae of silt. Meso-scale brittle deformation occurs in the form of shallow and steep-angle thrust faults and associated folding, but there is little evidence of deformation in the underlying sediments except for minor soft-sediment deformation and possible water escape or drainage structures.

A thin section (Fig. 18b and c; Nirex 1997g, k), shows that the lamination is highly disturbed by soft-sediment deformation, including micro-scale flame structures, crenulation-style folds and rounded clay-rich till ‘pebbles’ that locally occur within the silt laminae. Thin sand laminae within the seam are highly deformed, displaying lobate, irregular and wispy flame structures that deform upper and lower boundaries and cause local attenuation. Detachment of a prominent sand lobe or ball has caused intense, localized deformation within underlying sand, silt and clay, which itself has been injected upwards to accommodate lobe formation. Open, synformal, convolute fold structures are preserved within the lobe, deforming fine lamination. A slight convex-upward warping of clay and silt laminae is visible in the upper part of the thin section (Fig. 18b and c), whereas the laminae are offset by several small microfaults that have both normal and reverse senses of movement.

In general, the base of element MDc is sharp, forming an abrupt, laterally continuous and extensive sixth-order erosion bounding surface (fifth-order glaciotectonic surface) that commonly overlies subglacially deformed sediment (i.e. glaciotectonite sensu Benn and Evans 1996; Evans et al. 2006; Evans 2018), as shown in Figures 1012. Underlying sediment, however, is not invariably deformed, as at location x in Figure 10h.

Interpretation. The observed fissility displayed by interfingering of layers of Dmm and Dcm lithofacies and their stiff nature together with observed ‘solid-state’ meso- and micro-deformation structures such as shears, suggest that the occurrences of element MDc probably represent a component of subglacial traction till where dilation was ceasing, or had ceased, and a coherent till matrix was forming, or had formed. Such features may be a response to progressive formation and collapse of the till-matrix framework within a developing ‘B’ horizon (Evans et al. 2006, 2016; Evans 2018). With each successive framework collapse the intensity of bed-parallel foliation increases, resulting in each plate forming a localized subzone of ‘B’ horizon. This whole process is typically driven by repeated fluctuations in porewater pressure, resulting in the accretion of layers of variably deformed, strain-hardened, clast-supported ‘B’ horizon till separated by planar detachment, or erosion surfaces underlying individual, thin layers of massive, homogeneous, matrix-supported ‘A’ horizon till (Evans et al. 2006), as observed at Drigg, Ravenglass and Lowside Quarter. The observations of intra-stratification and deformation indicate that the seam at Lowside Quarter is a glaciotectonite sensu Evans (2018) (glaciotectonic architectures are described in the next section). However, Phillips et al. (2007) recorded the focusing of subglacial deformation within the water-saturated, stratified base of a till unit near Gretna, Scotland. Indeed, the zone of disturbed laminated silt and fine sand in the central part of the thin section appears to have undergone partial liquefaction and remobilization with the fine sands having been injected such that they locally cross-cut the lamination. This suggests that the zone here may have acted as a fluid pathway during deformation with water or sediment escaping from the underlying sequences via the water escape conduits flowing along the actively deforming zone. The seam also shares many characteristics of hydrofractures, some of which may stretch laterally (subhorizontally) and conformably for several metres within a sediment sequence, before turning abruptly up or down discordantly (Phillips and Merritt 2008; Phillips et al. 2013a, b).

The convolute or disharmonic nature of the folding is consistent with deformation of water-rich sediment resulting in localized liquefaction, whereas the presence of water escape or drainage conduits within sediment underlying element MDc at Drigg Beach (Figs 10f and 17d, f) and Lowside Quarter (Fig. 18a) probably records also liquefaction and water escape in the structurally lower parts of the sequence resulting from loading caused by the tectonic thickening of the sequence during high water saturation. The latter probably contributed a water-rich surface on which the ice moved, thus significantly reducing the amount of shear translated into underlying sediments, confining deformation to within element MDc. This resulted in little or no ductile or brittle deformation of the underlying sediment (see Evans et al. 2006; Phillips et al. 2007). Indeed, if there is a higher than critical porewater pressure, the development of a till matrix framework may be limited.

Although several, not necessarily mutually exclusive interpretations of this element have been provided, in the context of this paper, the importance of this deliberation is to demonstrate that great care is required in interpreting thin, horizontally laminated units within glacigenic sequences. Even if they are formed of stacked graded couplets, they may not be glaciolacustrine in origin.

Pervasive and non-pervasive mélange architectural elements Mp and Mnp

Although sediment underlying massive diamicton elements MD and MDc is locally undeformed (e.g. Figs 10h and 18a), these elements commonly rest on zones of deformed sediment (Fig. 6). For example, pervasive mélange element Mp directly underlies MDc at Sellafield (Fig. 13), Drigg Beach (Fig. 10f, location h; Fig. 10g, location t; Fig. 10h, location w; Fig. 17), Seascale (Fig. 9) and Ravenglass (Figs 12 and 15). At these locations the mélange element comprises sediment that has been affected by brittle deformation, resulting in both shallow-angle and steep thrusts together with conjugate sets of extensional faults. Ductile deformation has also occurred, resulting in rootless and sheath folds, boudins, overthrust material and commonly an S1 fabric that has obliterated the original sedimentary structure and locally is itself homogenized into a mélange of material. The fabric is formed of rotated and realigned clasts.

Sediments of element F5 appear to have undergone the most intense deformation. However, the exposures of element Mp in the vicinity of location h at Drigg Beach (Fig. 10f) also incorporate less intensely sheared clay of element Rh, as shown in detail in inset 3 (Fig. 17b). Additionally, as shown in insets 1 (Fig. 10f) and 4 (Fig. 17c), the northerly verging thrust duplex described above also incorporates thrust and displaced units of element Mp. Within that duplex, the S1 fabric has been itself deformed, resulting in an S2 crenulation fabric. Some exposures of Mp show distinct water escape fractures and conduits (Fig. 17d). Element Mp was observed to invariably underlie massive diamict elements MD and MDc, and, where present, to overlie non-pervasive mélange element Mnp (Fig. 6).

Element Mnp is widespread throughout most exposures, incorporating brittle deformation features that include shallow- and steep-angled thrusts with associated folds and ramp anticlines and conjugate extensional faults (Figs 7, 8, 1013, 15 and 17).

The original sedimentary structure was observed at all locations where element Mnp was identified (hence the ‘non-penetrative–non-pervasive’ sobriquet). Deformation features could be roughly divided into two sets: those that affect only sediment underlying massive diamicton elements (and are thus truncated on the underside of those units), and those that also affect the diamictons. For instance, shallow-angle thrust faults and associated folds were observed to have deformed sediment of element F5 at Seascale (Fig. 5c and d) and Drigg Beach (Fig. 10h, locations s and u), but are truncated by the sixth-order basal surface of overlying units of MD and MDc. Some structures have propagated up through that surface (e.g. location t, Fig. 10g) and deform the overlying MD and MDc elements, but are truncated at the base of overlying outwash sediments (e.g. locations b–d at the previously described thrust duplex, t and u; Fig. 10f–h), whereas other structures were observed to have deformed the whole sequence up to the seventh-order surface marking the base of post-glacial sediments. For instance, steep-angled, southerly verging thrusts and back thrusts deform the exposed sequence at Drigg Beach up to the base of the uppermost diamicton (Sandy Acre Till Member; SAT) (Fig. 10h, location u), forming a southerly verging, positive half-flower structure with tight ramp anticlines deforming a fining-upward sequence of gravel and sand unit of PPG (representing element F5), and the underlying diamicton (FWT) representing element MD. Some of the steep thrusts appear to have propagated up through the sixth-order surface into the overlying diamicton (SAT), which thus represents a younger occurrence of element MD.

Further south along Drigg Beach, beyond the large normal fault, a large set of low-displacement (<50 cm shortening) thrust faults deform a unit comprising rhythmites of clay and silt (we refer to this as SACS, representing architectural element Rh) and overlying sequence of sands and gravels of elements F5 and F6 (Fig. 10a and b). These themselves have been folded and possess a structural dip towards the SW. Compressional structures and folds were also observed throughout the sequence, where they are bisected by extensional faults, including a large strand with a displacement of c. 30 cm and a northerly dip (Fig. 10b), which is probably related to the larger extensional fault. Between the two faults, grey clay and overlying peat infills a depression, onlapping underlying sediment and clearly overstepping linked conjugate, extensional fault strands.

Interpretation. Zones of glacial mélange element Mp result from the brittle and ductile pervasive modification of existing sediment (such that original sedimentary structure is obliterated) during the subglacial development of an overlying traction till. The observations above are consistent with a Type A (pervasive–penetrative) glaciotectonite (sensu Benn and Evans 1996). The S2 crenulation foliation observed within the thrust duplex described above was probably emplaced during the subglacial formation of the duplex and records the deformation of the S1 fabric.

The thin and laterally discontinuous occurrence of element Mp, and its local dislocation from the underside of the overlying subglacial traction till, probably results from lateral fluctuations in subglacial water pressures. Translation of shear into sediments underlying till sheets occurred during lower water pressure, and none during higher pressures (see Phillips et al. 2007). The much thicker zones of purely brittle, non-pervasive deformation (here termed ‘mélange element Mnp’), which commonly underlie element Mp, record the formation of a Type B (non-penetrative) glaciotectonite (sensu Benn and Evans 1996). The variable presence of brittle deformation also probably results from fluctuations in subglacial water pressure. Observations of structural relationships commonly help in the identification of Type B glaciotectonite and thus provide additional supporting evidence for the identification of overlying subglacial traction tills.

Discussion 1: palaeoenvironmental reconstruction, glaciolacustrine v. glacigenic?

We agree with Coleman et al. (2020) that their new ‘glaciolacustrine’ interpretation presented by Coleman et al. (2020) has significant implications for future 3D modelling and site investigations, both in the district and in other areas of comparable glacial history. For example, by adopting their glaciolacustrine model, these researchers have projected the enigmatic thin seams of pebbly clay (diamicton) exposed at the coast near horizontally inland. In contrast, previous studies have demonstrated that they gently rise inland towards former glacial limits and take into account the former direction of ice flow (Nirex 1997l). For example, in early phases of the last glaciation ice flowed through, and directly out of Lower Wasdale, but during subsequent readvances the Irish Sea Ice Stream became progressively dominant and glided more or less horizontally across deposits buried in the valley (Merritt and Auton 2000).

Distinguishing between glaciolacustrine and glacigenic sediments presents similar challenges to those encountered in correctly identifying glaciomarine deposits, as opposed to non-marine subglacial or fluvial outwash. For example, McCarroll (2001) concluded that sedimentological evidence alone is not enough to identify glaciomarine sediments, particularly if they were deposited in a glaciotectonic setting. The issue is not unique to the Cumbria area and has been discussed in relation to other areas around the UK and Ireland with similar glacial history (e.g. east Yorkshire, north Norfolk, north Wales, Irish Sea coast; Phillips et al. 2007; Thomas and Chiverrell 2007; Evans et al. 2013; Lee et al. 2017; Merritt et al. 2018). Rather than seeking unequivocal sedimentological evidence the key to successful discrimination lies in the correct identification of glaciotectonic as opposed to gravity-driven deformational features and understanding the role of subglacial processes in sediment deposition, accommodation, accumulation and preservation (Hart and Roberts 1994; Roberts and Hart 2005).

The evidence that Coleman et al. (2020) cited against subglacial deposition of pebbly clay units in the Drigg Beach area includes sedimentological, glaciotectonic and accumulation style. We examine each facet in turn.

Sedimentology. Coleman et al. (2020) interpreted particle size distribution curves of the controversial units to indicate gap-graded clays of glaciolacustrine origin, especially as they assumed that the well-dispersed lonestone pebbles are dropstones. Coleman et al. (2020) invoke the presence at both macro- and micro-level of finely draped structures and gradational transitions at most of the sand–clay boundaries, like those reported as diagnostic features of glaciomarine sediments by Eyles and Clark (1985) to support their interpretation. In addition, they interpreted pillow and ball structures such as described by Murton and French (1993) to exclusively indicate lacustrine sedimentation. Furthermore, they suggested that the water-escape structures they observed in thin section are the result of glaciolacustrine sedimentation and proposed that thin wispy stringers of silt and sand within the pebbly clay units are similar to those interpreted by Eyles and Eyles (1983) to have resulted from water current re-sedimentation.

In reply, it is generally accepted that tills have the most diverse range of particle size distribution (psd) of any soil (Clarke 2018). The psd curves of Coleman et al. (2020) fall partially within the published ranges for many subglacial clayey diamictons, including, for example, the Brewood Till Formation in Staffordshire (Culshaw et al. 2017), till and sheared clay from thrust planes collected at Dinas Dinlle, on the Llyn Peninsula (Harris et al. 1997), and tills within the Vale of York Formation (Culshaw et al. 2017).

Pillow and ball structures such as those described by Murton and French (1993) can be created during lacustrine deposition, but the features described by those researchers are typically meso-scale (i.e. centimetres to metres) in size, and usually form when sand with a higher bulk density sinks into mud with a lower bulk density, rather than the other way round as Coleman et al. (2020) described. Furthermore, such features are more comparable with the complex, lobate, soft-sediment deformation structures described in a broadly similar and correlative readvance till near Gretna (British Geological Survey 2006; Phillips et al. 2008) and other locations (e.g. Lee et al. 2017).

Mass flows may be introduced locally into a glaciolacustrine environment and others may be generated during periods of exposure (Brodzikowski and Van Loon 1987). However, ultimately the extent and thickness of the mass flow-body depends on the extent of the glaciolacustrine basin and on a range of factors that act to limit the flow-body's preservation. As Merritt and Auton (2000) demonstrated, the units of pebbly clay diamicton in question form units with subplanar boundaries that are tens to hundreds of metres in length yet commonly less than 1 m thick. They can be lithostratigraphically correlated over tens of kilometres via outcrop and borehole data (Fig. 2). Units that are unequivocally lacustrine clay are more geographically confined and commonly much thicker. For instance, the Whinneyhill Coppice and Holmside Clay members located beneath Lower Wasdale are 2–3 km in extent and up to 30 m thick (Merritt and Auton 2000).

Lonestones are indeed very common within the pebbly clay units, but clast-poor subglacially deformed materials often appear this way and in the absence of cited diagnostic evidence of impact and draped lamination (see Benn and Evans 2010) most are more likely to have resulted from subglacial processes (see Piotrowski 1992; Hart and Roberts 1994; Roberts and Hart 2005).

Glaciotectonics. Coleman et al. (2020) cited that an apparent absence of glaciotectonic deformation structures (micro-scale folding and faulting) both within the pebbly clay units and in underlying sediment supports subaqueous deposition (see Eyles and Eyles 1983). Furthermore, they suggested that observed downward-dipping clasts and poorly oriented fabrics in the pebbly clays indicate deposition as mass flows.

In contrast, we have observed numerous meso- and macro-scale deformation structures at Drigg that affect the pebbly clay units together with underlying sediment, and understand that structures such as these, and poorly developed tectonic fabrics, are a common feature in the wider study area (Nirex 1997e, g, k). We have explained how such deformation features are generally associated with the passage of grounded ice over wet sediment and that they are commonly partitioned into discrete structural zones of deformation (see Phillips et al. 2007), as shown in the palaeoenvironmental reconstruction in Figure 19. However, although the presence of these structures is a useful criterion for defining a subglacial traction till (Evans et al. 2006), their local absence is not a useful indicator to distinguish subglacial and lacustrine deposition.

Accumulation style. Coleman et al. (2020) argued that preferential accumulation within topographic lows indicates gravity-driven sedimentation. However, although some fine-grained units towards the top of the sequence (e.g. Sandy Acre Clay and Silt Member; SACS) undoubtedly occupy basins and exhibit thin, graded, draped, onlapping lamination (Fig. 10b–d) characteristic of lacustrine sedimentation, others do not. We maintain that three of the units are much more likely to represent subglacial traction till. At Drigg Beach, the boundaries of the pebbly clay beds are generally non-gradational (location y, Figs 10f and 17) and thin stringers of silt and sand were not observed within them.

Summary. Our study demonstrates that it is crucial to ensure distinct glaciotectonic evidence (including published information) is not overlooked, where the interpretation of ice-marginal/ice-contact successions is used as the basis of robust 3D geological models, particularly where it can support unconvincing sedimentological evidence. They then used the apparent absence of glaciotectonic deformation to support their glaciolacustrine depositional interpretation without further investigation. In contrast, the comprehensive structural and micropetrological evidence we present strongly indicates the presence of at least three subglacial traction till units interbedded with clastic deposits, which fully supports the generally accepted glacial readvance scenario based on considerable evidence reported by Nirex and summarized by Merritt and Auton (2000).

The controversy regarding the ponding of ice-marginal lakes and associated glacial readvances during latter stages of the last glaciation in west Cumbria is not unique in the UK. For example, analogous large ice-marginal lakes also occupied the Durham lowlands, particularly in valleys of the Tyne and Wear, and Teesside (Smith and Francis 1967; Smith 1994; Stone et al. 2010; Davies et al. 2012). The largest was Glacial Lake Wear, which stood at several levels up to 132 m OD, governed by the elevations of spillways that became available sequentially and ponded along the coast behind ice streaming southeastwards in the North Sea basin. The lakes were filled with up to 60 m of mainly thinly laminated clay and silt with dropstones, intercalated with pebbly clay diamicton and gravel (Tyne–Wear Glaciolacustrine Formation of McMillan et al. 2011). Similar to sediments beneath Lower Wasdale, complex interdigitation of units is common, especially towards the coast, and understanding the composition and 3D distribution of the fine-grained glaciolacustrine deposits is also important, particularly because they pose major geotechnical problems. The relatively high water content of the laminated silts and clays renders them plastic and potentially very weak with low shear strength (Smith 1994).

As in west Cumbria, these glaciolacustrine deposits are widely capped by thin (<2 m), laterally persistent beds of soft, plastic silty clay with well-dispersed pebbles (‘Pelaw Clay’, ‘Prismatic Clay’) (Smith 1994). Elsewhere they are capped by more clast-rich diamicton (‘Butterby Till’), but with low shear strength, inverted strength profile and local absence of a preferred fabric (Smith and Francis 1967). As with the glaciolacustrine model of Coleman et al. (2020), these units were generally believed not to be subglacial tills, but to have formed as cohesive debris flows into water from ice margins, as solifluxion flows from surrounding deglaciated slopes, or by redistribution of unconsolidated lake sediments following periodic lake drainage (Smith and Francis 1967; Hughes et al. 1998; Smith 1994). However, more recent research in the region, particularly at the coast, has firmly established that the ice-marginal lakes existed and that they were affected by several surges of the North Sea Ice Stream (Davies et al. 2009, 2012).

Discussion 2: heterogeneity and its application to high-resolution subsurface 3D geological modelling

Three levels of heterogeneity have been recognized in the study area (sensu Slomka and Eyles 2013). Level 1 heterogeneity is recorded on a vertical scale of centimetres to metres, whereas Level 2 is on a lateral resolution of <10 m to several hundreds of metres and a vertical resolution of <1 m to several metres. The lowest resolution of heterogeneity, Level 3, is recorded on a vertical scale of <1 m to several metres and a horizontal scale of up to several hundred to thousands of metres (i.e. the scale of Figs 1 and 2). Level 3 generally encompasses lithostratigraphic units and is identified at outcrop through integrating the topology, geometric characteristics and genetic relationship of component architectural elements constrained within fifth-, sixth-and seventh-order bounding surfaces. Correlation between borehole logs based on Level 3 heterogeneity and using fifth-order bounding surfaces will consequently result in the lowest resolution of 3D geological model.

Finer grained facies within the sheet-like FA1 and FA4 facies associations form lower permeability confining layers that have in places been eroded out beneath younger channelized units of higher permeability FA2, and FA5. This results in vertical hydraulic connectivity between facies associations of different lithostratigraphic units, and, for example, results in ‘perched’ aquifers above the main groundwater aquifer at Sellafield (Kuras et al. 2016). The geometry of the sixth-order surfaces that bound facies associations ranges from virtually flat to low-relief, long-wavelength, long-amplitude undulations. Modifications to this geometry result from subglacial or glaciotectonic deformation, which may lead to localized ‘shrink-and-swell’ relationships, removal or stacking. However, understanding the complex arrangements of coherent lithological units, their extent and their geometries, is the key to ascertaining where potentially contaminated groundwater may flow. Of paramount importance is the location and geometry of fine-grained, impermeable lithologies with respect to coarser material. Guidance such as BS5930:2015, together with corresponding lithology codes (e.g. the c. 120 AGS codes of the Association of Geotechnical and Geoenvironmental Specialists), provides some basis for correlation, but the possibility of up to 10 depositional sequences could result in a modelling stratigraphy of over 1000 units, which would be challenging for even the most enthusiastic 3D correlation-based modeller.

We address this issue here by mapping the permeability ranges, predominant flow mechanisms and ranges in hydraulic conductivity for various lithological sequences onto lithofacies and/or relevant architectural elements and facies associations that are of interest to the hydrogeological modeller (as per Lewis 1989; Lewis et al. 2006) (see Fig. 20). This demonstrates that most of the end-members of permeability that are of interest can be correlated and modelled at Level 2 (Slomka and Eyles 2013), apart from units of gravel, sand and silt within architectural elements MD (convoluted), F5 and F6, which ought to be correlated at Level 3 heterogeneity. However, field evidence including graded third-order bounding surfaces between silt and sand lithofacies units within element F5, the complex arrangement of fining-upward, small (metre-scale) sand and gravel units within element F6, and the chaotic array of various, similarly small lithofacies units (resulting from cryoturbation) with a huge range of upper and lower bounding surface geometries within element MD (convoluted), mean it is sensible to base correlation on defining the Level 2 bounding surface geometries. For elements F5 and MD (convoluted), we base the permeability range and hydraulic conductivity on the dominant lithology denoted by the borehole lithology log. Neither gravel nor sand is dominant within element F6, so we have assigned the sand and gravel hydraulic conductivity and permeability ranges of Lewis (1989) and Lewis et al. (2006), as in Figure 20.

The parameters of Lewis (1989) and Lewis et al. (2006) suggest that for the definition of architectural elements, borehole lithology logs may be rationalized down to bulk lithology, whilst utilizing full log descriptions to inform unit definition and correlation of architectural elements and their tops and bases. In Figure 21, we define 21 generic architectural element geometries including undeformed units and those that have undergone secondary subglacial or proglacial deformation. These have been used along with the architectural element analysis to facilitate the creation of several small-scale (hundreds of metres) 3D geological models of the Quaternary sequence at the Sellafield site using the workflow shown in Figure 22a-l. They are based on the gradual development of a 3D fence model of correlation cross-sections (a)–(e) constrained by a topography upper surface (in this case we have used the most recent Environment Agency aerial LiDAR digital terrain model) and a bedrock (rockhead) surface created in this study from interpolation of borehole and outcrop data. The architectural element analysis within this study illustrates that a typical glacial advance–retreat–readvance cycle comprises a fine-grained till at the base followed by a fining-upward sequence of glaciofluvial material (Fig. 22g) representing high-energy, ice-proximal deposition with upward-decreasing flow. However, this cycle may be repeated, because prograding and retrograding glacial advance and retreat episodes may result in apparently coarsening-upward sequences, such as at Peel Place Quarry. The gradual creation of a 3D fence model of sections within this study illustrates some of these complexities, which are represented by the development of an appropriate 3D geological modelling stratigraphy encompassing the bounding surfaces of over 400 units (Fig. 22h). Although this may seem to be a large number, it considers both vertical and lateral heterogeneity, and captures features such as braided river systems, glacially dammed and oxbow lakes, and till deposits, all of which may have undergone syn- and post-depositional deformation. In this study, the resulting 3D fence diagram (Fig. 22i) subsequently allowed the creation of a high-resolution 3D surface model (Fig. 22j) that could be explored in various ways, including looking at individual architectural elements such as till units (Fig. 22k) on their own, or in the context of other units such as sands (Fig. 22l).

To increase the effectiveness of high-resolution 3D modelling of geologically complex areas, such as around Sellafield, Lelliott et al. (2009) recommended utilizing a high density of site investigation data points within the dataset, based on the premise that a model is only as geologically valid as the data used to create it and that ‘something is usually better than nothing’; missing out data points just increases uncertainty, unless the geology is simple. Kessler et al. (2009) and Lelliott et al. (2009) noted that each cross-section should be constructed considering data points that may not lie along the correlation section line by projecting nearby boreholes onto the section line. Furthermore, such data points can be used to inform the creation of new sections that bisect them. This technique has been successfully deployed by the British Geological Survey within their GSi3D (Geological Surveying and Interpretation in 3D) workflows (Kessler et al. 2009; Lelliott et al. 2009; Smith 2010) and Groundhog software. Quantifying uncertainty is important where the data density is low, but calculations for this type of exercise are relatively easy to undertake (Lelliott et al. 2009). For example, uncertainty estimates calculated for each geological surface generally match the intuitive interpretation of the geologist building the model, who decides where the greatest uncertainty occurs. Lelliott et al. (2009) noted, however, that the calculated estimates required validation through the drilling of additional boreholes, meaning that an uncertainty estimate exercise can effectively form the basis of a data gap analysis, alongside other statistical approaches such as the Data Quality Objectives (DQO) process (US Environmental Protection Agency 2006).

In this paper we recommend that a reliable, high-resolution 3D geological model can be built following thorough data review and utilizing all available data points, and indeed have constructed several ourselves, one of which is shown in Figure 22 (another is described by Kuras et al. 2016). Reducing data decreases uncertainty in lithostratigraphic correlation, but removal of such data takes away possibly the only source of information on the lithology of a unit, its distribution and geometry. The case example published recently by Coleman et al. (2020) took a more simplified, less time-consuming approach to 3D geological modelling at the UK's LLWR, incorporating fewer data. These researchers did not include, for example, data such as historical records not completed to current standards such as BS5930:2015, data not supported by geotechnical records such as cone penetration tests (CPT), or data collected by potentially less reliable investigation techniques relative to rotary cored boreholes. We agree that ‘understanding and modelling the complex glacial geology in West Cumbria is challenging’ (Coleman et al. 2020, p. 1), but it is essential to build upon all available knowledge for developing robust hydrogeological models. Coleman et al. did not cite, for example, any of the detailed and comprehensive research work undertaken for Nirex in the 1990s (Nirex 1995, 1996, 1997a, b, c, d, e, f, g, h, i, j, k, l, m).

The new wholly glaciolacustrine interpretation proposed by Coleman et al. (2020) for the succession of clays and outwash sands and gravels exposed in the cliffs of the study area, and concealed beneath Lower Wasdale, is an interesting supposition that those researchers then extrapolated to the LLWR site, some 500 m distant. They integrated it with lithological and geotechnical information obtained from boreholes drilled within and in the vicinity of the LLWR site, geophysical data and regional geomorphological information, and then used it in the construction of a 3D geological model of the site. For this they deployed a prioritized borehole review approach resulting in the initial utilization of data of only 64% of 746 boreholes (Coleman et al. 2020), coupled with a revision of the site ‘event stratigraphy’ that aligned with their new interpretation and that therefore precluded the possibility of subglacial traction till being present.

This paper, which describes the application of architectural element analysis of local outcrop analogues to high-resolution 3D geological modelling for the nearby Sellafield nuclear site, has, by necessity, included a critical review of the evidence on which the glaciolacustrine hypothesis rests. We find it inappropriate to dismiss the generally accepted glacial readvance scenario. In contrast, by collecting detailed sedimentary, structural and micropetrological evidence, we demonstrate the unequivocal presence of up to three well-exposed subglacial traction tills, in a succession that includes glaciofluvial and glaciolacustrine sediments. We suggest that this provides a much more reliable basis for the 3D geological modelling of the nuclear sites within the study area.

This paper presents a detailed analysis of the sedimentary and glaciotectonic architectural characteristics of fine-grained glacigenic and fine- to coarse-grained glaciofluvial deposits exposed at field locations in and around the Sellafield nuclear works, Cumbria, England. In total, 19 lithofacies units have been recognized within 11 sedimentary architectural elements and including two glaciotectonic–mélange elements that result mainly from modification of existing sediment during emplacement of overlying subglacial diamicton.

We show how the high-resolution lithological and geometrical information required for contaminant transport modelling at complex industrial sites such as nuclear plants can be obtained through 3D geological modelling via correlation of erosion bounding surfaces at an architectural element level of heterogeneity. We present a new workflow showing how the architectural elements recognized in the Sellafield area have been used to develop a suitable 3D modelling stratigraphy and a high-resolution 3D geological model. The latter includes the definition of 21 generic unit or bounding surface geometries at local outcrops to help inform accurate correlation between data points (mainly boreholes) at the Sellafield site.

This study confirms the generally accepted evidence presented by Merritt and Auton (2000) for the presence of at least three laterally extensive, subglacial traction tills along the coast between Sellafield and Ravenglass. These units are closely associated with two glaciotectonic mélange architectural elements, which represent penetrative (Type A) and non-penetrative (Type B) glaciotectonites (sensu Benn and Evans 1996; Evans 2018) respectively.

These elements are not laterally continuous, and whereas the first element (Mp) always occurs in tandem with an underlying zone of the second (Mnp), both may be absent whilst locally only Mnp may be present.

This is consistent with published literature, but importantly it shows that whereas the presence of a sub-till glaciotectonite is indicative of subglacial traction, its absence does not negate the unit above being subglacial traction till.

At Drigg Beach, a remarkable NW-verging thrust duplex is preserved, resulting in overthrusting of two members and the creation of a second S2 crenulation fabric within overthrust Type A glaciotectonite. This indicates a local lowering of subglacial water pressure during emplacement that resulted in coupling of the ice mass to its bed (sensu ‘sticky spots’ of Stokes et al. 2007). Elsewhere along the coast, clay-lined hydrofractures have been identified that were formed by pressurized water, which is likely to have played a significant role in controlling subglacial processes. The local uncoupling of ice from its base (sensu Lesemann et al. 2010) may have provided accommodation for the subglacial deposition of thin, lensoidal units of glaciofluvial and rhythmically laminated sediment that commonly lie in the former ‘down-ice’ lee of anticlinal structures affecting subglacial till units, creating a shrink-and-swell architecture. At least two episodes of proglacial compressional and extensional deformation have been identified.

The depositional and structural findings presented in this paper support the lithostratigraphy and history of events outlined by Merritt and Auton (2000) and Smith and Merritt (2008), including deposition and deformation associated with several glacial readvances and not the revised, wholly glaciolacustrine model proposed by Coleman et al. (2020).

This paper is published with the permission of the National Nuclear Laboratory, UK and the Executive Director of the British Geological Survey (UKRI). The authors thank Sellafield Ltd for provision of borehole data, and for permission to reproduce the photograph used in Figure 13. Tendley Quarries Ltd are thanked for permission to access Peel Place Quarry during several field visits over the last 15 years and for use of images of the quarry. The paper benefited from reviews by an anonymous reviewer and J. Lee (British Geological Survey), who are both thanked for their comments and suggestions. N.T.S. was supported by the National Nuclear Laboratory (NNL) Core Science Research and Development (R&D) Programme (Decontamination Science Theme). J.W.M. and E.R.P. publish with permission of the Director of the British Geological Survey. J. Hyde, J. Griffiths, M. Smith, A. Bull, A. Wareing, C. Lennon, S. Kwong, L. Edmiston, J. Graham and D. Connor (all NNL), and J. Taylor-Bray (ONR) are thanked for assistance in the field and for informal manuscript review.

NTS: conceptualization (lead), data curation (lead), formal analysis (lead), funding acquisition (lead), investigation (lead), methodology (lead), project administration (lead), resources (lead), software (lead), validation (lead), visualization (lead), writing – original draft (lead), writing – review & editing (lead); JWM: data curation (equal), formal analysis (equal), funding acquisition (equal), investigation (equal), validation (equal), visualization (equal), writing – original draft (supporting), writing – review & editing (equal); ERP: formal analysis (equal), funding acquisition (equal), investigation (equal), validation (equal), visualization (equal), writing – original draft (supporting), writing – review & editing (supporting)

This research received no specific grant from any funding agency in the public, commercial or not-for-profit sectors.

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Some borehole logs analysed during the current study are not publicly available owing to nuclear site restrictions but may be made available from the corresponding author on reasonable request. Please note some historical field exposures are now lost owing to erosion, quarrying or construction.

This is an Open Access article distributed under the terms of the Creative Commons Attribution 4.0 License (