Hydrogeology Using numerical modelling to test the geological and groundwater conceptual understanding of a complex, layered aquifer: A case study from the Fell Sandstone, Northumbria

.


Introduction
Sandstone aquifers form globally important groundwater systems. Their relatively high porosity, and therefore storage, provides useable groundwater resources in many parts of the world. Examples include the Nubian sandstone aquifer in north Africa (MacDonald et al., 2012), the Great Artesian Basin in Australia (Hayes et al., 2019), the Guarani of South America (Hirata & Foster, 2020), and parts of the US (e.g. High Plains Aquifer; (Korus Jesse et al., 2022) ). In the UK, Permo-Triassic sandstone aquifers are of significant importance for water supply purposes, providing the second most exploited groundwater resource after the Cretaceous Chalk aquifer. Whilst Permo-Triassic sandstones are the most widely studiede.g. (Colyer et al., 2021;Lafare et al., 2021), with both Permian age (MacDonald et al., 2003) and Triassic-age sandstones (Medici et al., 2019) receiving attentionother, minor sandstone aquifers also contain important water resources for local communities. One such A C C E P T E D M A N U S C R I P T example is the Carboniferous age Fell Sandstone located in north east England, which has provided important water supply for the Berwick-upon-Tweed area since the early 1900's (Hodgson et al., 1971), as well as being exploited for geothermal purposes under Newcastleupon-Tyne (Younger et al., 2016). Belying its relatively small outcrop area, it is geologically and hydrogeologically complex (Turner et al. ,1992) and has not had the commensurate attention that other sandstone systems have received to understand it fully.
This paper provides a summary of the geological and modelling activities conducted to develop an updated conceptual groundwater flow model of the Fell Sandstone in the Berwick-on-Tweed area. We also present the first attempt of applying distributed numerical modelling to simulate groundwater dynamics of the Fell sandstone, a complex layered aquifer, which has been used to test the updated conceptualisation, and to assess the sustainability of the abstractions. The main findings resulting from these recent studies demonstrate the complexity of this groundwater system and offer lessons for similar layered sandstones both in the UK and worldwide.

Summary of previous investigations of the Fell Sandstone in Northumberland
The first geological assessments of the Fell Sandstone were published in the middle of the twentieth century by Robson (1956) and Smith (1967) focussing on the sedimentological nature. However, the earliest comprehensive hydrogeological study was provided by Hodgson et al. (1971). They examined meteorological, hydrogeological and geochemical data to develop a conceptual understanding of groundwater flow in the Fell Sandstone in the Berwick area, including Berwick itself, which has been subject to groundwater abstraction ( Figure 1) since the late 1930s. A recharge value of 392 mm/a was estimated for the whole Fell Sandstone outcrop, and it was concluded that the overall water balance was showing a surplus. Hodgson et al. (1971) also report results of pumping tests in two abstractions boreholes ("Thornton" and "Borehole A") showing transmissivities in the order of 50 m 2 /d.
Groundwater flow direction was assumed generally towards Berwick to the north-east, while chemical analysis of groundwaters samples indicated mainly neutral pH and chemically moderately hard waters with limited evidence of saline intrusion in Berwick. An area of reduction in the main chemical determinants (e.g. TDS) was observed around Borehole A thought to be related to pumping and induced recharge (Hodgson et al., 1971). Based on the results of these investigations, Hodgson et al. (1971) concluded that there was potential for further exploitation of the Fell Sandstone and in particular, the suitability for the provision of a groundwater supply for the Berwick area. Turner et al. (1992) then built on previous geological and hydrogeological studies to provide conceptualisation of the aquifer system consisting of a sequence of seven distinct sandstone A C C E P T E D M A N U S C R I P T units separated by mudstone units. From the top to the bottom of the stratigraphic sequence, the sandstone units were named "South Ord", "Murton Crags", "Murton Deane", "Peel Knowe", "Middle Ord", "Royalty" and "Thornton Park". This conceptualisation has been the mainstay of any hydrogeological investigations since.
There is some evidence suggesting karst behaviour in the Fell Sandstone, which was first identified and reported by Mullan (1989). Some indications of karst behaviour in boreholes as a result of fracturing and removal of weakly cemented sandstone are also reported in the BGS/EA Minor Aquifer Properties Manual (Jones et al., 2000).The theme of karstification was further developed by Self et al. (2005), who made the case for cave development and provide proof of the existence of sinking streams.
Over recent decades there has been a significant amount of work undertaken by environmental consultants on the conceptualisation of the aquifer system. Two main studies, commissioned by the Environment Agency (EA) and undertaken by Entec (now WSP, previously Wood and AMEC), presented a conceptual model of groundwater flow in the Fell Sandstone (Entec, 2009) and a revised recharge model and water budget using new data (AMEC, 2012). These studies were incorporated into one undertaken by Northumbrian Water Limited (NWL) as part of the National Environmental Programme (NEP) work. This study (NWL, 2018) builds on the AMEC work along with an MSc thesis completed at Imperial College, London (Colyer, 2018) and provides a more refined geological and hydrogeological understanding of the area. In terms of hydrochemistry, the studies of (Entec, 2009), (NWL, 2018, (Dearlove & Schwartz, 2022) and the more detailed EA nitrate report (2016) presented a more detailed discussion of the background groundwater geochemistry as well as the nitrate pollution in the aquifer.
Before this study, there have been few applications of groundwater modelling to study the Fell Sandstone aquifer. However, these only focussed on the estimation of the water balance without considering distributed models of groundwater flow dynamics. Examples of these previous modelling applications include a systems dynamic approach (Markou, 2013) and a lumped parameter model (Colyer (2018) based on AquiMOD (Mackay et al., 2014).
The latter showed that the Fell Sandstone response to recharge has significant spatial variability.

Location
The Fell Sandstone is situated in Northumberland, north east England, with Berwick-upon-Tweed to the north-east (Figure 1). The outcrop describes a thin curve adjacent to the Cheviot Hills with topography rising from sea level to a maximum of 440 m above Ordnance Datum (a OD). The outcrop forms a craggy escarpment with dip slopes to the south-east leading down to the River Tweed and eventually the North Sea. The land use is predominantly arable except the urban centre of Berwick. Long-term average annual rainfall for the period January 1987 to December 2018 is 695 mm/a, with variability in precipitation controlled by the topography. The main river in the area is the River Tweed, which is fed by a number of small tributaries running off the outcropping areas of the Fell Sandstone.

Geology
The region's bedrock succession and glacially dominated superficial geology is depicted in 1:50 000 scale map data (BGS, 1977) and summarised by Stone et al. (2010). The Carboniferous-age Fell Sandstone Formation has been the subject of extensive sedimentological studies (Smith, 1967;Turner et al., 1987;Turner et al., 1992;Turner et al., 1993). In the Berwick-upon-Tweed area, the 225 to 350 m thick Fell Sandstone mostly comprises interbedded sandstone, siltstone and mudstone. Seven main sandstone units and several lesser sandstone-dominated layers are recognised. Sandstone units vary in thickness (up to 20 m) and lateral extent. The sandstone units are commonly medium-to coarse-grained, moderately rounded, moderately sorted and ubiquitously cross-bedded. Mudstone units observed in drill core (e.g. Murton Craggy Bogs Observation borehole 1) are dark reddish brown, and lack internal structure (Dearlove et al., 2022).
The Fell Sandstone is underlain by the Ballagan Formation, a succession of thin and thick sandstones, mudstones and siltstones, with 'cementstones'. Due to lithological similarities between the formations (as reported in borehole records) and a lack of outcrops in the study area, the base of the Fell Sandstone is poorly constrained in the geological cross-sections.
The Fell Sandstone is overlain by the Scremerston Formation, comprising alternations of sandstone, siltstone, mudstone and coal, with occasional dolomite or limestone beds.
Descriptions of the Fell Sandstone that are based on its full geographic extent (BGS Lexicon) cite limestone, coal and seatearth (deposits underlying coal seams) as 'subsidiary' or 'trace lithologies'. However, a detailed review and correlation of borehole records indicates that these lithologies are absent from the formation in the study area. Hence, for this work the base of the Scremerston is taken below the lowest recorded interval of limestone, coal (or associated seatearth, fireclay etc.).
In the study area to the south of Berwick-upon-Tweed, the Fell Sandstone generally dips towards the southeast by 8° to 12°. An asymmetrical southeast-plunging anticline in the Berwick-upon-Tweed area is associated with dips of up to 25°. Two mapped faults exist in the study area to the north and south of the Bleak Ridge borehole, cutting west-east. These structures are inferred to be normal and have downthrows of up to 40 m to the north. Current structural mapping in the study area underrepresents the extent of faulting and associated deformation. It is highly likely that unmapped faults are present, and that these may limit the lateral continuity of the sandstone units or juxtapose stratigraphically different units.
The superficial geology of the study area ( Figure 3) broadly consists of: glacial till (semicontinuous sheets in low-lying parts of the study area and irregular patches on the eastfacing slopes); patches of sand and gravel (associated with the till areas); head (ribbons between till deposits and areas of no superficial cover, especially in the south of the study area); river alluvium and accumulations of peat (occupying modern-day drainage system and topographic depressions, including the Thornton Bog area).

Hydrogeology
The main groundwater system consists of sandstone bodies within the Fell Sandstone, which form different aquifer units separated by less permeable mudstone units. Flow in the sandstone is assumed to be predominantly intergranular, but there is some evidence of fractures in the core as well as karstic features, which have been identified further south towards Wooler, i.e., Routin Lynn Cave (Mullan, 1989).
Relatively few pumping tests have been undertaken in this area, e.g. see Jones et al. (2000). S C R I P T Ambient regional groundwater flow direction is difficult to discern from analysis and interpolation groundwater heads measured in boreholes given the relatively sparse monitoring network and the influence of abstractions. Historic analysis presented in Hodgson et al. (1971) seems to suggest the groundwater flow direction is north/north-east towards Berwick. However, this simplistic interpretation is unlikely to represent the current state of the flow field, given the layered nature of the aquifer and the cones of depression caused by the abstractions.
A conceptual groundwater balance capturing the main components can be defined as follows (e.g. (Fetter, 2018)): The main inflow is likely to be rainfall recharge, predominantly direct, i.e., occurring where precipitation directly infiltrates at outcrop. Secondary to this is the runoff from other adjacent areas (SW inflow) and any exchange with other groundwater units. The main outflow is likely to be NWL abstraction and other local abstractions with secondary outflows being springflow and limited baseflow to streams.
Infrequent spot gauging is carried out at the streams draining the area, namely Horncliffemill Burn and Newbiggin Dean. It has been reported that groundwater levels are below the base of the stream bed and the streams could be losing flow to the sandstone units (Environment Agency (2016).
The main features of the conceptual model are how "layered" the system is thought to be with sandstone units separated by mudstone. The main inflows are rainfall recharge and subsurface runoff from the mudstone. Groundwater flow is along the sandstone units, dipping down towards the centre of the basin. Faulting may play a part in restricting flow within the aquifer.

Geochemistry
Groundwaters are typically Ca-Mg-HCO3 type, with some exceptions including higher chloride, sulphate and sodium-potassium water in localized areas. The piper diagrams presented in Entec (2009) show relatively uniform chemical compositions of the borehole waters, while those presented by NWL (2018) seem to indicate a wider heterogeneity.
In terms of impact of anthropogenic activity on groundwater quality, nitrate pollution is an issue in the catchment (Dearlove et al. (2022)). Elevated concentrations of chloride are also observed in boreholes mainly in the urban area of Berwick, but values are very much lower than sea water concentrations and have been interpreted as not likely related to saline intrusion (Entec, 2009). The Environment Agency (2016)

Methods and materials
The following section describes the activities conducted to revise the geological understanding, to quantify aquifer recharge and estimate the water balance, and to implement a numerical groundwater flow model of Fell Sandstone in the area immediately west of Berwick.

Revision of the geological understanding
The first step in the revision of the geological understanding was to identify relevant data sources and collate them into a systematic dataset (see Table 1). An initial set of over 350 boreholes was identified and then a sub-set identified having depth greater than 10 m.
These >10 m deep boreholes were coded by an agreed internal BGS protocol. Coding of the boreholes recorded as much sedimentological information (grainsize, composition texture, structure etc.) as possible to aid subsequent correlation. Of the initial borehole records identified, 83 have been used in the construction of the schematic cross-sections and revised 1:50,000-scale geological linework. The understanding of the 3D geology was supplemented with the examination of cores where available as well as fieldwork involving visits to quarries and other surface outcrops. A number of field sedimentary logs were recorded at the available outcrops in the region (Murton High Crags and two localities in Shoreswood; Quarry and field cutting). These sedimentary logs were taken to identify the sedimentological architecture of the exposed sandstones. Building on the 83 coded boreholes as well as other information, schematic cross-sections were created using BGS in-house software 'Groundhog' 1 using a Bald Earth digital elevation model (DEM) as the capping surface.

Recharge modelling and water balance
Recharge was estimated with the distributed recharge model ZOODRM (Mansour et al., 2018). ZOODRM uses the UN Food and Agriculture Organisation (FAO) method (Allen et al., 1998) to estimate evapotranspiration and excess soil water, and the excess water is split into surface runoff and recharge according to run-off coefficients. Surface runoff is routed across the land surface and can re-infiltrate into non-saturated soil: for example, in this case, runoff generated on the mudstone may infiltrate into downslope sandstone outcrops.
The datasets used in the recharge modelling are shown in Table 2 incorporated into the groundwater model. The recharge was averaged over 6-month periods (January to June and July to December). The MODFLOW-2005 numerical model was implemented in two stages: a preliminary model based on previous geological understanding was initially implemented to test the interaction of the main sandstone units (Bianchi et al., 2019). This was followed by a revised second model (Bianchi et al., 2020), which took into account the novel geological understanding (Ford et al., 2019) and considered an expanded modelling domain to the south-west including the area surrounding the NWL Borehole F (Figures 1 and 5).

Groundwater flow model
The model domain was discretised with a block-centred finite difference grid consisting of square cells with sides equal to 25 m, organised in 156 rows and 336 columns. In the vertical direction, the domain was discretised into seven layers each representing a single sandstone unit. The correspondence between model layers and current geological interpretation and historical nomenclature of the sandstone units is as follows: • Layer 1: "South Ord" • Layer 2: "Murton Crags" • Layer 3: "Murton Dean" • Layer 4: "Peel Knowe" • Layer 5: "Middle Ord" • Layer 6: "Royalty" • Layer 7: "Thornton Park" The mudstone units separating the sandstone units ( Figure 5) were simulated as "quasi 3D confining beds". With this choice, groundwater flow is not explicitly simulated, but their and S with depth. The later was used to represent vertical variation in hydraulic parameters within a single layer of a groundwater model, e.g. (Rushton et al., 1989).

Geological understanding
Collectively, the units within the Fell Sandstone were deposited within a fluvial-deltaic setting, with the sandstone units representing the higher energy deposits (compared to the mudstone) generated within the river channels and on the delta that these channels feed into . The mudstones within the Fell Sandstone are either the product of floodplain deposition (associated with the channel sandstones) or were laid-down within shallow marine environments. However, it was not possible to discriminate between these two types of mudstone, though the ability to do so could allow for more accurate subdivision and correlation of the Fell Sandstone. At a larger scale, the separation of the succession into predominately mudstone-and sandstone-dominated layers (e.g. Table 3) is likely driven by cyclical changes in the sea level during the period of deposition of the Fell Sandstone (e.g. Turner et al. (1997)). The

A C C E P T E D M A
N U S C R I P T start of the deposition of the sandstones (their bases) likely represents periods of relatively low sea level, whereby the rivers of the Fell Sandstone flowed out over the continental shelf for greater distances to reach the ocean into which they drained. The fall in relative sea level would have also had the effect of causing the fluvial system to incise and erode underlying material. Some of the mudstones including the laterally more persistent interleaving mudstones between the sandstones (Table 3)

Net recharge and water balance
Long-term average distributed recharge generated by ZOODRM can be seen in Figure 7.
Recharge on the outcrop varies between 150 and 225 mm/a. This range is somewhat lower than previous assessments (e.g. Infiltration of 392 mm/a in Hodgson et al. 1971). Estimated values for the components of the water balance extracted from the recharge model are presented in Figure 8 where Details of the long-term groundwater balance can be found in Table 4. It should be borne in mind that these estimates are for the whole recharge model area, and, whereas all of the abstractions are from the Fell Sandstone, not all of the recharge will reach the Fell Sandstone. Therefore, the water balance is for an area larger than the groundwater abstractions obtain their inflow. The recharge produced by ZOODRM is potential recharge and doesn't take into account the distribution of the Fell sandstone outcrop, therefore that reaching the water table (actual recharge) is calculated within the groundwater flow model.
There is a large amount of uncertainty in the estimation of outflow as baseflow/springs. The estimate was derived from the two spot gauging sites. Owing to the limited amount of data and the fact that groundwater heads are generally dominated by large inter-annual trends rather than seasonal fluctuations, annual baseflow was estimated as the annual minimum gauged flow. The long-term average baseflow/spring flow was then defined as the median annual baseflow and assumed to be constant across the whole model area. The two subcatchments associated with the spot gauging sites comprise 33% of the total area and 45% of the Fell Sandstone outcrop area of the recharge model area. Inflow/outflow from/to other groundwater units was assumed zero because of a lack of information. However this assumption may lead to water lost or gained from the aquifer, particularly from flow outside the model boundary, but also from the Scremerston or the Ballagan.

A C C E P T E D M A
N U S C R I P T

Simulated groundwater heads and budget.
The optimal values for the input parameters estimated after calibration of the MODFLOW model are presented in Table 5 Whilst these are different from those derived from pumping test analyses, they are consistent with the lower end of the range of transmissivities reported in the literature. It is, however, to be noted that pumping tests analyses have a large degree of uncertainty associated with them related to measurement errors, design of test, availability of observation boreholes and its interpretation. Therefore, they are treated as guide values for groundwater modelling as they rarely produce transmissivity values that can be placed directly into a regional model.  Table 4). Other reasons include but are not limited to erroneous pumping data since the modelled hydrograph does not follow the behaviour of the observed hydrograph with respect to sudden changes in the pumping rate over time.
Overall, the mismatch between simulated and measured values measured in term of RMSE is equal to 8.8 m. However, when the most problematic boreholes (i.e. the bottom four in Table 6) were excluded from the analysis of the residuals, the average RMSE falls to a value around 4.5 m. This value is considered an acceptable margin of error considering the complexity of the aquifer system and the uncertainty in the interpretation. In particular, the main sources of uncertainty are the geometry and continuity of the sandstone and mudstone units in certain areas of the model domain, the lack of data for the definition of the boundary and initial conditions, particularly the lack of knowledge regarding the natural initial state of the flow field predating the abstractions.
Model estimates of the groundwater budget are presented in Figure 10. Over the simulated period (30 years), the total abstracted groundwater volume (62,634 Ml) is about 7% lower than the cumulative recharge (67,228 Ml). This discrepancy shows that to balance the model water balance, inflow is required from across the boundary via the GHB and that the size of the model may need to be increased. However, the recharge water balance (see The sustainability of supply from the groundwater system has been investigated.
The main findings from the work are as follows: (1) The system is complex, with variation in spatial interactions between the sandstone units . There is a need to conceptualise the system as seven discrete sandstone units, unusually for a British aquifer and its associated groundwater system, particularly of this size. The interaction of different sandstone units vertically can be important. Generally, they are separated by very low K mudstone units. However, this relationship changes spatially and certain sandstone units, namely Murton Dean and Royalty, have limited extents as well as direct interconnection in places.
(2) Differing time responses between variations in the hydraulic stresses (i.e. recharge or abstractions) and the resulting changes in groundwater heads exist. There are also differences in temporal response between units with a general slow accumulation of recharge (accumulating over six months). However, some hydrographs show very rapid response, e.g. Murton Craggy Bog which can see a response to storm events (Dearlove et al., 2022). Other boreholes such as The Kells demonstrate a very unusual groundwater response which requires further investigation.
(3) The geology and hydrogeology of the Fell Sandstone is complex, with significant vertical and lateral interactions in a small area. Whilst the Scremerston overlays the Fell Sandstone and is considered as a non-aquifer, there remain questions as to its behaviour. One of these is how permeable its lower formation is and, in particular, whether this horizon produces inflow to borehole A.
Despite the recent investigations, there are uncertainties that could be reduced through further study: including the 3D geometry of the units; the further delineation and investigation of the role of the faults; stream-aquifer interactions; and a water balance for the whole area.
The next planned phase of work is to widen the area of interest both toward Berwick and further south-west. This will take the form of combined geological-hydrogeological fieldwork alongside groundwater modelling. It is expected that, by continually iterating the geological and numerical models along with the conceptual understanding of groundwater flow to address particular questions in specific locations, an holistic understanding can be formed to assist both the water company in its exploitation of the aquifer and the regulator in its oversight.   Figure 1)    Outflow to other groundwater units