Groundwater flow velocities in a fractured carbonate aquifer-type: Implications for contaminant transport

Contaminants that are highly soluble in groundwater are rapidly transported via fractures in mechanically re- sistant sedimentary rock aquifers. Hence, a rigorous methodology is needed to estimate groundwater flow velocities in such fractured aquifers. Here, we propose an approach using borehole hydraulic testing to compute flow velocities in an un-faulted area of a fractured carbonate aquifer by applying the cubic law to a parallel plate model. The Cadeby Formation (Yorkshire, NE England) – a Permian dolostone aquifer present beneath the University of Leeds Farm - is the fractured aquifer selected for this hydraulic experiment. The bedding plane fractures of this dolostone aquifer, which are sub-horizontal, sub-parallel and laterally persistent, largely dominate the flow at shallow (<~40 mBGL) depths. These flowing bedding plane discontinuities are separated by a rock matrix which is relatively impermeable ( K well-test / K core-plug ~10 4 ) as is common in fractured carbonate aquifers. In the workflow reported here, the number of flowing fractures - mainly bedding plane fractures - intersecting three open monitoring wells are found from temperature/fluid conductivity and acoustic/optical televiewer logging. Following well installation, average fracture hydraulic apertures for screened intervals are found from analysis of slug tests. For the case study aquifer, this workflow predicts hydraulic apertures ranging from 0.10 up to 0.54mm. However, groundwater flow velocities range within two order of magnitude from 13 up to 242m/ day. Notably, fracture apertures and flow velocities rapidly reduce with increasing depth below the water table; the upper ~10m shows relatively high values of hydraulic conductivity (0.30–2.85m/day) and corresponding flow velocity (33–242m/day). Permeability development around the water table in carbonate aquifer-types is common, and arises where high pCO 2 recharge water from the soil zone causes calcite/dolomite dissolution. Hence, agricultural contaminants entering the aquifer with recharge water are laterally transported rapidly within this upper part. Computation of groundwater flow velocities allows determination of the Reynolds number. Values of up ~1, indicating the lower limit of the transition from laminar to turbulent flow, are found at the studied site, which is situated away from major fault traces. Hence, turbulent flow is likely to arise in proximity to tectonic structures, such as normal faults, which localize flow and enhance karstification. The occurrence of turbulent flow in correspondence of such tectonic structures should be represented in regional groundwater flow simulations.

Fractured carbonate aquifers, which represent the focus of this work, underlie a land area covering~15% of the earth's surface, sup-plying~25% of world population need of drinkable water (Ford and Williams, 1989). Rigorous coupling of hydraulic monitoring (hydraulic head) and field tests (televiewer and/or fluid logs and hydraulic tests) is required to determine parameters (e.g., fracture hydraulic aperture and flow velocity) which control solute transport velocities in their saturated zones.
A methodology for computation of groundwater velocities within fractures has been developed by several authors (Tsang et al., 1990;McKay et al., 1993;Lo et al., 2014;Medici et al., 2016Medici et al., , 2018aMaldaner et al., 2018;Ren et al., 2018a). Initially, optical and acoustic imaging and fluid logging, or borehole dilution testing are used to determine the number of discontinuities that effectively flow. Next, hydraulic tests on screened intervals are undertaken for computation of fracture permeability and hydraulic aperture using the cubic law. At this stage of the workflow, a large number of authors have derived fracture flow velocities for characterization of fractured igneous and sedimentary aquifers, assuming a parallel plate model of fracture flow (Witherspoon et al., 1980;Ji et al., 2008;Quinn et al., 2011). Indeed, application of the cubic law to a parallel plate model is considered reasonable in absence of karstic cavities for fractured aquifers with flow dominated by sub-parallel discontinuities which are separated by an impermeable matrix (Witherspoon et al., 1980;Berkowitz, 2002). This is the approach of our research for determination of groundwater flow velocities and flow regime (laminar or turbulent as indicated by the Reynolds number) within the fault blocks of a fractured aquifer.
In our investigation, groundwater flow velocities within fractures in a carbonate aquifer under natural conditions have been found using optical and acoustic televiewer / fluid logging and slug tests, with focus on variation with respect to depth and seasonality. Our research contrasts with the previous approaches to determine velocities and type of flow regime for groundwater under natural conditions (Cappa et al., 2008;Quinn et al., 2011;Maldaner et al., 2018). Indeed, the focus is on the variation with both depth and time in a specific case of a nonfaulted area. We selected a non-faulted block at specific field site, the University of Leeds Farm (Spen Farm) in NE England between the cities of Leeds and York (see Fig. 1A-D), to orient future research in terms of hydraulic testing and modelling strategies for fractured aquifers.
The University of Leeds Farm is an experimental site dedicated to the Critical Zone, which is defined as the outer layer of earth's surface extending from the top of the vegetation cover down to the depth of fresh water circulation in the subsurface (Anderson et al., 2008). Research at this experimental site is focused on the migration of the impact of human activities, specifically agriculture in this case, on this outer layer which is of paramount importance for both human health and economic development (Banwart et al., 2011(Banwart et al., , 2012. In this paper, we report a hydro-geophysical characterization of the Critical Zone of the carbonate aquifer of the Permian Magnesian Limestone in Yorkshire (NE England, UK) at the University of Leeds Farm. This experimental farm site (see Fig. 1A-B) is underlain by this dolomitic aquifer and provides an opportunity to assess risk of contaminant transport in an area of the UK which is characterized by intense agricultural and farming activity. In fact, 80% of NE England is arable and livestock farmland; as a result the concentration of nitrates in groundwater has steady increased since 1940 due to introduction of nitrogen fertilizers (Forster et al., 1982). Previous research describes the Magnesian Limestone Group as a carbonate aquifer where flow is largely dominated by fractures (Cairney, 1972;Allen et al., 1997). Hence, the site at the University of Leeds Farm is an optimum laboratory to investigate fracture flow in an area where hazards from contaminant dispersal is relatively high (Foster and Crease, 1974;Forster et al., 1982).
Previous published studies on the hydraulic properties of the Magnesian Limestone aquifer in NE England have focused on core plug and pumping test analyses and hydro-geochemistry (Cairney, 1972;Aldrick, 1978;Younger, 1995;Allen et al., 1997;Mayes et al., 2005). More recently, a steady state groundwater flow model was published for an area further north in County Durham (Fig. 1B;Neymeyer et al., 2007). Here, we propose an investigation using a combination of experimental techniques including fracture analysis on optical and televiewer logs, monitoring of the water table, slug tests and temperature and electrical conductivity fluid logs, new to this aquifer.
Specific research objectives are as follows: (i) estimate hydraulic conductivities and hence apertures of flowing discontinuities, (ii) estimate groundwater velocities and Reynolds number and their variation with depth and time, and (iii) address future research needs for carbonate aquifer-types.

Geological setting
The Magnesian Limestone aquifer is characterized by carbonate rocks derived from shallow water sedimentation at the margins of the Zechstein Basin (Aldrick, 1978;Lee, 1993). In NE England, the Magnesian Limestone aquifer attains a typical thickness of 120 m ( Fig. 1A-D) and is formally part of the Zechstein Group, which is characterized by a succession of marine dolomite, limestone, evaporites, mudstone and siltstone of exclusively Permian age (Smith et al., 1986;Tucker, 1991). The Magnesian Limestone Group is formally divided into three different formations: the Cadeby, Edlington and Brotherton formations (Fig. 1B;Smith et al., 1986).
The Cadeby and Brotherton formations represent the mainly carbonate units; they are separated by the marls, mudstones and evaporites of the Edlington Formation. The base of the Cadeby Formation is characterized by 5 m of marls which separate the aquifer from the sandstone, mudstone and clays of the Pennines Coal Measures Group (Fig. 1B, C). The Brotherton Formation passes upwards into the Roxby Formation which is entirely characterized by anhydrite and gypsum (Allen et al., 1997).
The Cadeby Formation is characterized, with the exception of the basal 5 m of marls, by thinly bedded dolostone showing ooids, corals and bivalves (Tucker, 1991). This formation, which represents the focus of this study, was deposited during the Lower Permian (272-251 Ma) in a shallow marine environment under sub-tropical climatic conditions. Then, secondary dolomitization occurred at early (Permian) and middle (Triassic) diagenetic stages (Peryt and Scholle, 1996). Dolomitization fluids were supplied from the gypsum and evaporites of the overlying Edlington Formation (Smith et al., 1986). The Brotherton Formation above represents the uppermost part of the Magnesian Limestone aquifer in Yorkshire and is characterized by thinly bedded dolomitic limestone with ooids, algae and bivalves. The Brotherton Formation was also deposited in a shallow water environment under tropical conditions; but it differs from the Cadeby Formation due to less intense dolomitization (Peryt and Scholle, 1996).
In NE England around the University of Leeds Farm (Fig. 1B, D), the tectonic structures which characterize the UK Magnesian Limestone aquifer are represented by normal faults and non-stratabound joints (sensu Hartmann et al., 2007). Outcrop studies and seismic lines carried out near the field site show how normal faults of Mesozoic age are mainly oriented ENE-SSW ( Fig. 1B; Aitkenhead, 2002, Cooper andLawley, 2007;Medici et al., 2015). The Cenezoic Alpine orogenesis results in the gentle dip (< 5°towards E) of the Permo-Mesozoic deposits in the Yorkshire area (Fig. 1B, C;Bottrell et al., 2006;Hartmann et al., 2007;Bricker et al., 2012). Outcrop observations of the Cadeby and Botherton formations show how normal faults represent areas of high fracture density (Aldrick, 1978). However, non-faulted sections in Highmoor Quarry at the University of Leeds Farm site (Fig. 1D) and elsewhere show high angle (dip 50°-80°) joints of~1-2 m persistence cross-cutting bedding parallel fractures (Fig. 2;Walker, 2006;Boddy, 2018). Such high angle joints have been previously mapped in the Magnesian Limestone aquifer in Yorkshire (Kortas and Younger, 2013). These fractures arise from the Cenozoic uplift of NW Europe and are related to unloading due to removal of overlying strata (Odling et al., 1999;Kortas and Younger, 2013). Other discontinuities that characterize the Magnesian Limestone include sub-horizontal bedding plane fractures and stylolites. Bedding plane fractures in this aquifer are considered beds that have been reactivated by the Cenozoic uplift (Kortas and Younger, 2013). However, stylolites ( Fig. 2A, B) represent pressure solution structures that formed during burial diagenesis. Bedding plane fractures exceed the outcrop scale (~20 m) and are much more persistent than stylolites, which in the Cadeby Formation show persistence ranging from 0.5 up to 1.5 m ( Fig. 2A; Walker, 2006;Gillham, 2017).

Hydrogeological background
The Magnesian Limestone aquifer represents the fifth and the third most important aquifer in terms of abstraction volume in the entire UK, and NE England, respectively (Cairney, 1972 Abesser and Lewis, 2015). Contamination in this carbonate aquifer has arisen from the use of fertilizers in agriculture; nitrate concentrations in groundwater rose from average values of 0.28 mg/L to 75 mg/L between 1948 and 2006 in the study area (Aldrick, 1978;Walker, 2006). At the University of Leeds Farm (Fig. 1A-D) the Magnesian Limestone aquifer is unconfined due to absence or presence of only~1 m thick Quaternary cover (Gaunt et al., 1970;Aldrick, 1978). The Cadeby and Brotherton formations represent aquifers, separated by the Edlington Formation, which is considered a "leaky" aquitard. Hydraulic connectivity between the two separate formations is either related to faulting or gypsum dissolution that has resulted in the collapse of strata in the Edlington Formation (Farrant and Cooper, 2008).
The Magnesian Limestone Group is typically highly mechanically resistant (UCS of 48-75 MPa); it has been used as a building stone in NE England since Roman times (Lott and Richardson, 1997;Lott, 2013). In fact, Highmoor Quarry near the study site is one of the few main sources of dimension stone of this lithology, which is elsewhere typically heavy fractured in outcrop (Walker, 2006).
Groundwater flow is directed towards the east driven by topography which is characterized by higher elevations in the adjacent area of Coal Measures strata outcropping to the west (Fig. 1C). The Cadeby Formation interquartile interval ranges for intergranular hydraulic conductivity and porosity are 2.9 × 10 −4 to 0.9 × 10 −3 m/day and 8.5 to 18.7%, respectively. These ranges substantially overlap those of the Brotherton Formation, which range from 4.0 × 10 −4 to 1.5 × 10 −3 m/ day and 9.9 to 19.0% for intergranular hydraulic conductivity and porosity, respectively (Allen et al., 1997).
Pumping test values from the UK Magnesian Limestone aquifer are reported by Allen et al. (1997); transmissivities range from 6 to 4300 m 2 /day with median of 299 m 2 /day. These pumping tests highlight a lower transmissivity of the Cadeby Formation (T median = 163 m 2 /day; n = 16) with respect to the Brotherton Formation (T median = 420 m 2 /day; n = 15). This has been interpreted as resulting from the higher dolomitisation which characterizes the Cadeby Formation (54% CaCO 3 ; 46% MgCO 3 ); this reduces the rate of dissolution in correspondence of fractures. By contrast, the dolomitic limestone of the Brotherton Formation is more abundant in calcite (65% CaCO 3 ; 35% MgCO 3 ), hence more prone to dissolution and permeability development (Allen et al., 1997;Lott and Cooper, 2005).
However, the highest transmissivities in the Magnesian Limestone (T median = 2000 m 2 /day; n = 7) are associated with fault zones due to the high degree of fracturing (Allen et al., 1997;Cooper and Lawley, 2007). The transmissivities within fault blocks are generally lower but extremely variable depending on fracturing density and proximity to rivers (Allen et al., 1997). Flow occurs essentially in correspondence of fractures, i.e. evidenced by the large difference between permeability from pumping and core plug tests (K well-test /K core-plug~1 0 4 ; Allen et al., 1997). Scanline surveys which have been carried out in the quarries nearby the field site indicate development of fissures of a few centimetres aperture exclusively in correspondence of extensional faults (Walker, 2006). Water geochemistry in the Cadeby Formation in the area of Durham (Fig. 1A) shows relatively high alkalinity (380-450 mg/ L), and relatively high concentrations of sulphate and chloride from dissolution of halite and gypsum from the overlying Edlington Formation (Younger, 1995;Mayes et al., 2005).

Experimental methods
Three wells (BH1, BH2, BH3) were drilled in the lower part of the Magnesian Limestone aquifer, which is represented by the Cadeby Formation (NE England, UK) at the University of Leeds Farm (Fig. 1A, B). The three wells (30 to 40 m depth) are located within fault blocks avoiding known fault traces. The layout design approximates an equilateral triangle configuration as illustrated in Fig. 1D. The boreholes were drilled with diameters ranging from 0.23 and 0.28 m to install a maximum of three piezometer tubes, each with a 0.08 m diameter with screen length ranging from 2.3 to 3.0 m (see construction details in Table. 1).

Optical and acoustic televiewer logging
Optical and acoustic televiewer (Mount Sopris ALT QL40 mk5) logs were recorded in open boreholes prior to piezometer installation, and the data were analysed by structure picking (sensu Williams and Johnson, 2004) of discontinuities. This structure picking allowed sterographic plotting of discontinuities (n = 264) as well as the computation of mean vector statistics for each group of fractures using the Stereonet 9 software package (Allimendinger et al., 2012). Discontinuities recorded in logs in the Cadeby Formation have been divided into four groups: bedding plane fractures (D1), stylolites (D2) and SE (D3) and NW (D4) dipping high angle joints.

Fluid temperature and conductivity logs
Fluid temperature and electrical conductivity logs have been recorded in the three open boreholes using a Geovista mk2 multi-parameter probe. Measurements of fluid temperature and electrical conductivity were undertaken under ambient conditions. Flow logging was realised prior to installation of piezometer tubes in order to define sections dominated by flowing fractures for location of screened intervals (construction details in Table 1). Boreholes were logged both downwards and upwards but the downward fluid logs are presented here, as these have better data quality.

Slug tests
Three slug tests have been performed in each of the screened intervals of BH1, BH2 and BH3 boreholes (Table 1) using 5, 2 and 5 l as withdrawal volume following Butler (1998). A Solinst Levelogger Edge 3001 data logger was used to measure hydraulic head variations. A preliminary slug test was carried out in each screened interval to estimate the return time of the initial displacement to 5% of rest water level. Based on this preliminary test, time intervals between the slugs were fixed to 10 min for BH1 and BH2. However, a longer period of 30 min was used for BH3.
Slug tests analyses in BH1, BH2 and BH3 have been realised using the AquiferWin32 software package. Horizontal hydraulic conductivity (K h ) has been computed for each slug (n = 24) using the methods of Bouwer and Rice (1976) and KGS (1985) applicable to partially penetrating wells in unconfined aquifers.

Hydraulic head monitoring and determination of groundwater flow vector
Hydraulic gradient and direction of the groundwater flow have been estimated by hydraulic head monitoring in each piezometer from 1st October 2017 to 30th June 2018 using a sampling frequency of 1 h, correcting for barometric pressure variations. Position and elevation of the top of the wells (see Table 1) were surveyed with a Leica GPS TS15 device for accurate determination of hydraulic gradient. The accuracy of this GPS device is ± 10 mm and ± 20 mm for lateral and vertical position, respectively.
The water table was assumed as planar and the three point estimation method has been used for determination of the groundwater flow vector. The arithmetic averages of the piezometers heads in each borehole were used as variations between piezometers within the same borehole, i.e. differences in head are ≤0.5 mm. Hourly data of hydraulic gradient magnitude and azimuth of groundwater flow have been calculated with the 3PE Spreadsheet (Beljin et al., 2014).

Potential recharge
Potential recharge (R) is defined as the difference between total precipitation (P) and the potential evapotranspiration (ET 0 ). Daily totals of rainfall (P) are available at the field site, i.e. a rain gauge is installed at the National Centre for Atmospheric Science (NCAS) weather station (see location of the station with respect to BH1, BH2 and BH3 wells in Fig. 1D). Air temperature is recorded every 30 min at this weather station. This allowed computation of potential evapotranspiration (ET 0 ) using the Hargreaves Eq. (1), where R a is extraterrestrial radiation, T m is the daily mean air temperature, computed as an average of the maximum and minimum air temperatures, TR is the temperature range between minimum and maximum values. Extra-terrestrial radiation, R a , is computed using HYDRUS-1D from the latitude and altitude of the field site.

Hydraulic aperture determination
Average fracture hydraulic aperture has been found in order to determine the average groundwater flow velocities in fractures intersecting each screened interval. The hydraulic aperture (b) represents the aperture of a parallel walled fracture of equivalent permeability. This hydraulic aperture is expressed by Eq. (2) based on the implementation of the cubic law assuming laminar flow (Romm, 1966) where T s is the screened interval transmissivity, g is the gravitational acceleration, v the kinematic viscosity of water (1.307 × 10 −6 m 2 /s based on the 9.5°fluid temperature at the three boreholes) and N the number of flowing fractures intersecting the screened interval. Flow in fractures intersecting each screened interval is considered bed-parallel and sub-horizontal as around 90% of flowing fractures are bedding plane fractures with < 4°dip (see section 4.2 below). Screened interval transmissivity (T s ) is the product of the hydraulic conductivity (K h ) from slug tests and the piezometer screen length (w). Eq. (2) is commonly used for computation of the hydraulic aperture in both fractured carbonate and crystalline aquifers at depths <~50 mBGL (Quinn et al., 2011;Ren et al., 2018a).

Computation of groundwater fracture velocity and flow regime
Average horizontal flow velocity (V) in fractures intersecting each of the studied screened intervals can be computed using Eq. (3) which is derived from the cubic law assuming fractures as parallel plates (Snow, 1965;Romm, 1966;Qian et al., 2011).
where b is the average hydraulic aperture, g is the acceleration due to gravity, v is the kinematic viscosity of water and i is the hydraulic gradient along the fracture direction. Here, the average velocity (V) in fractures was computed for each screened interval after determination of the hydraulic aperture (b) using Eq. (2) from slug tests, and the hydraulic gradient from the groundwater flow vector (see Section 3.4). After determination of flow velocity at each screened interval the Reynolds (Re) number was found using Eq. (4) to assess whether laminar or turbulent flow occurs under natural conditions (Ji et al., 2008;Ji and Koh, 2015;Ren et al., 2018a). Reynolds numbers < 1 indicate laminar flow, 1-10 transitional and > 10 turbulent flow, i.e. laboratory experiments with single fracture models indicate that non-Darcian flow is significant with Re > 1-10 ( Witherspoon et al., 1980;Zimmerman et al., 2004;Ranjith and Darlington, 2007;Ji et al., 2008),

Aquifer fractures
Outcrop observations at Highmoor Quarry (see Fig. 2) show how the dolomite aquifer of the Cadeby Formation is thinly bedded, i.e. the typical vertical spacing of bedding planes is 0.4 m in quarry faces. High angle joints (50°-80°) form a non-stratabound fracturing system, i.e. they typically cross-cut the bedding planes ( Fig. 2A). The fracturing network in the three boreholes is representative of a similar non-faulted area. Optical and acoustic televiewer logs do no show displacement of bedding planes or presence of cataclasites confirming the validity of previous geological surveys at the field site (see Fig. 1D).
Good quality optical and acoustic televiewer logs (see Fig. 3) in the three studied wells allowed the structure picking of discontinuities, which are represented by bedding plane fractures (D1), stylolites (D2) and sub-vertical joints (D3, D4), see Fig. 4. Low values of amplitude and travel time in acoustic televiewer logs show where bedding plane fractures (D1) and joints (D3, D4) are open, hence representing potential flow pathways (Figs. 3A, B, D). Here, optical and acoustic televiewer images show no large (> 1.5 mm) karstic features in correspondence of bedding plane fractures and joints (Figs. 3A, B, D). Stylolites (D2) are typically sub-horizontal like the bedding plane fractures ( Fig. 2A, B). Despite this, they can be recognized in outcrop and optical and acoustic televiewer logs due to their typical irregular shape (Figs. 2B, 3C; Park and Schot, 1968;Rolland et al., 2014). Acoustic televiewer logs also show how stylolites are characterized by lack of travel time anomalies. Thus, such discontinuities are substantially closed and likely unconducive for flow.
The optical and acoustic televiewer logs acquired from the three wells confirm that the characteristics of the Cadeby Formation observed in outcrop are also present in the subsurface. The three vertical boreholes show a layered aquifer; 92% of discontinuities are sub-horizontal (< 30°). As a consequence, relatively high angle (50°-80°) fractures characterize only 8% of total discontinuities. However, note that moderate and high angle joints are under-sampled in vertical boreholes. The Terzaghi (1965) correction has been applied in contouring to correct for bias due to the vertical orientation of the three studied boreholes.

Fluid logs
Fluid temperature and electrical conductivity have been logged under ambient conditions in all three wells BH1, BH2, BH3 (Figs. 1D, 5, and Fig. S1). Fluid temperature ranges from 9.2°up to 10.0°; the vertical profiles do not increase with depth showing lack of correlation with the geothermal gradient. Fluid electrical conductivity shows values up to 110 mS/m ( Fig. S1 and Fig. 5); relatively high for potable waters. These values arise from natural alkalinity from carbonate dissolution but also gypsum and halite dissolution in the overlying Edlington Formation; ions flow into the underlying Cadeby Formation though normal faults which juxtapose the two geological formations in the field site area ( Fig. 1C; Aldrick, 1978;Mayes et al., 2005).
The three logged boreholes show several inflow horizons identified by changes in gradient in fluid temperature and electrical conductivity in correspondence with visible traces in the optical and acoustic logs G. Medici, et al. Journal of Contaminant Hydrology 222 (2019) 1-16 (Figs. 5, S1). The number of inflows identified are 5, 7 and 6 for BH1, BH2 and BH3 boreholes, respectively. Each inflow horizons is either characterized by single bedding plane fracture or by a cluster of 2-3 discontinuities which are vertically spaced < 0.10 m. Effective flowing fractures represent between 30% and 100% of the visible features on the optical and acoustic logs sections corresponding to each screened interval (Table 3). The principal flow pathways in the Magnesian Limestone aquifer are represented by the bedding plane fractures (D1) such as those represented in Fig. 3. Indeed, sixteen of the eighteen inflows identified in the three boreholes (BH1 and BH2, BH3) are characterized by bedding plane (D1) discontinuities; nearly all such inflow points are either represented by a single major bedding plane discontinuity (Fig. 3B) or a cluster of minor bedding plane fractures (Fig. 3D). Another, single inflow in BH2 corresponds with two high angle joints (D4); one inflow occurring in BH3 corresponds with both a low angle bedding plane fracture (D1) and a high angle joint (D2). Thus, 90% of the total flowing fractures (n = 28; Fig. 3) are bedding plane fractures with < 4°dip. Although high angle joints (Fig. 3A) do not represent the principal flowing structures in the aquifer of the Cadeby Formation, our data indicate they are capable of flow, hence establishing vertical hydraulic connections between bedding plane fractures. No inflows are localised in correspondence of stylolites (D2). To summarise: flow in the Cadeby Formation is largely dominated by bedding plane fractures (D1), sub-vertical joints (D3, D4) secondarily contribute and stylolites (D2) seem to be non-conductive.

Hydraulic conductivity
Horizontal hydraulic conductivity (K h ) values from slug tests are summarized in Table S1 and range from 0.07 to 2.85 m/day showing medians of 0.49 and 0.57 m/day computed using the Bouwer and Rice (1976) and KGS (1985) methodologies, respectively. Comparison of hydraulic conductivities reported in Table S1 indicates that values analysed with Bouwer and Rice (1976) are 2% to 23% lower than those obtained using the KGS (1985). 1 The latter method has been selected for computation of hydraulic aperture, i.e. specifically developed for fractured anisotropic aquifers (Hyder and Butler, 1995). This fits our hydraulic characterization of the Cadeby Formation, i.e. temperature and electrical conductivity logs show how sub-horizontal bedding plane fractures largely dominate flow.
Hydraulic conductivity computed both using Bouwer and Rice (1976) and KGS (1985) techniques progressively decrease with depth moving from shallower to the deeper screened intervals in all three boreholes. This hydraulic conductivity decrease has been plotted in Fig. 6 which reports the KGS (1985) values as arithmetic means of the tests in each screened interval.
These screened intervals (BH1 P1, P2; BH2 P1-3; BH3 P1-3) show length ranging from 2.3 up to 3.0 m and are characterized by a similar number (n = 2-5) of sub-horizontal bedding plane fractures (D1). Such sub-horizontal discontinuities represent the principal flow pathways in the Cadeby Formation and are characterized by a relatively regular vertical spacing (0.38 m; σ = 0.57 m; Table 2) in the three boreholes. The number of fractures that localize inflows does not decrease with depth (Table 3). Hence, the observed hydraulic conductivity decrease (Fig. 6) cannot be related to vertical variation of fracturing density.

Hydraulic monitoring 4.4.1. Hydraulic head fluctuation and potential recharge
The hydraulic monitoring realised at the University of Leeds Farm shows a similar pattern of head variation in the three boreholes, with head range of 2.5-4.0 m between low in December 2017 and high in  April 2018 (Figs. 7, S2), with almost identical results from loggers installed in piezometers in the same borehole (max difference 5 mm), suggesting good hydraulic connectivity between bedding fractures. Temporal changes in hydraulic head are consistent with the potential recharge estimated from the NCAS Weather Station data, i.e. a gentle fall in head which reaches the minimum at the end of December 2017 corresponds to a period of negative potential recharge (R), see Fig. 7A and B. A time lag of three weeks is observed between potential recharge becoming mainly positive in early December and the minimum hydraulic head occurring at the end of that month. Positive potential recharge (R) from 26th December up to 5th April (Fig. 7A), arises from low values of potential evapotranspiration (ET 0 ) occurring during this period. This matches a steep rise in hydraulic head during the winter-spring period which is~2.5 and~4 m in BH2 and BH1/BH3

Table 3
Fractures detected by the optical and acoustic televiewer logs vs. flowing fractures picked using fluid and temperature fluid logs within each screen interval.  respectively (Fig. 7B, Fig. S2). The hydraulic head reached a maximum in the middle of April (Fig. 7B); the time lag between the switch in potential recharge towards negative values in early April and the peak hydraulic head is~10 days. Finally, a drop of the water table characterizes the period of late April, May and June when high values of potential evapotranspiration (ET 0 ) exceeds rainfall (P) in this period of the hydraulic year (Fig. 7A, B).

Groundwater flow vector
Hourly data of hydraulic gradient and azimuth (found using the averages of hydraulic heads from the piezometers in each borehole, see Figs. 7C, D, 8A, B) show that the hydraulic gradient ranges from 0.0210 to 0.0247 and the arithmetic mean is 0.0236 (Fig. 8A), i.e. relatively minor seasonal variation. The azimuth varies by only 5.8°with a Fisher Mean of 68.4° (Fig. 8B).
The hydraulic gradient is characterized by a very gentle rise during the months of October, November and December (Fig. 7C). This plateau during the autumn produces a cluster of hourly data in the range 0.0240-0.0247 (where 45% of hourly data are grouped). The hydraulic gradient drops during the months of January, February, March and early April (Fig. 7C); corresponding with the rise in hydraulic head at the field site. Then, it rises contemporaneously to decreasing hydraulic head in late April, May and June (Fig. 7B, and Fig. S2).
The direction of the groundwater flow vector is also characterized by minor seasonal variation (Fig. 7D); the azimuth of the vector rotates 5.8°eastwards during the steep rise in hydraulic head during the winter-spring period, and shows a bimodal distribution (Fig. 8B), with 48% and 23% of readings clustered between 66.5°-67.5°and 71.0°-72.2°, respectively. The 66.5°-67.5°cluster is representative of the months of October, November and December; the 71.0°-72.2°cluster is representative of end of April, May and June which are characterized by both a gentle drop of the water table and a little eastward rotation of the groundwater flow vector (Fig. 7B, D). This contrasts the hydrologic scenario of the winter-spring months which are characterized by both a steep rise of the water table (~3 m; Fig. 7B) and a sharp eastward rotation of the flow vector (Fig. 7D).

Hydraulic aperture, fracture flow velocities and regime
Average fracture hydraulic aperture (b) at each of the screened intervals is computed using Eq. (2); b ranges from 0.10 up to 0.43 mm (Table 4), with smaller values found for increasing depths. Such values are in the range of shallow (<~100 mBGL) fractured aquifers for bedding plane discontinuities and joints (~0.01-0.8 mm; Cappa et al., 2008;Quinn et al., 2011;Maldaner et al., 2018;Ren et al., 2018a).
Determination of hydraulic aperture allows computation of the average flow velocity within fractures using Eq. (3), based on the measured horizontal hydraulic gradients. Groundwater flow velocities computed using this approach range from 13 up to 242 m/day (Fig. 9). Velocities rapidly decrease with increasing depth in all three boreholes; as a consequence of the depth sensitivity of hydraulic conductivity and hydraulic aperture (see Fig. 6). The first~10 m below the water table represents the part of the Cadeby Formation which is characterized by relatively high horizontal hydraulic conductivity ( Fig. 6; K h = 0.30-2.85 m/day) and flow velocity ( Fig. 9; V = 33-242 m/day).
Groundwater flow velocities in correspondence of fractures also vary over time, due to variations of hydraulic gradient during the period of the groundwater recharge (Fig. 9, indicated by the horizontal bars). However, this effect is only maximum 15% of the mean value, whereas velocities differs by several times (> 10 in BH2 and 3) between the upper and lower screened intervals.
The Reynolds (Re) number found using Eq. (4) varies, in proportion with flow velocities (Fig. 10), from 0.01 up to 1.1. Thus, laminar flow occurs under natural conditions. Notably, the Reynold number in the BH1-P1 screened interval, which is located only~2 m below the water table, partially exceeds the lower limit of the field of laminar/turbulence transition (Fig. 10). G. Medici, et al. Journal of Contaminant Hydrology 222 (2019) 1-16 K core-plug~1 0 4 ) with flow dominated by sub-horizontal (< 4°) bedding plane discontinuities with absence of large karstic features (> 1.5 mm) as schematically represented in Fig. 11. Inferred hydraulic fracture apertures and hence groundwater flow velocities in this dolostone aquifer of Permian age show strong reduction with increasing depth below the water table. Similar reduction of aperture of discontinuities with depth has been observed in quarries in the Jurassic dolostones of the Venetian Alps in northern Italy and in the Carboniferous limestones of the County Claire in western Ireland (Williams, 1983).
Our findings for the Cadeby Formation of NE Yorkshire are consistent with trends in transmissivity (T = 15-4300 m/day; n = 23) from historical pumping tests in that T values do not show any correlation with well depths up to~50 m (Allen et al., 1997), indicating that the permeability with the formation usually resides close to the water table. Such permeability distributions typically arise from dissolution of calcite and dolomite occurring in the shallowest part of aquifers immediately below the water table, which strongly enhances permeability and fracture apertures (Allen et al., 1997;Göppert and Goldscheider, 2008;Hartmann et al., 2014). The driving force of carbonate dissolution in carbonate aquifers is represented by elevated pCO 2 recharge water percolating through the soil zone (Atkinson, 1977;Palmer, 1991). This is an active process in the Critical Zone at the University of Leeds Farm site. In fact, research on soil geochemistry at this site has revealed that agriculture activities increase decomposition rates of soil organic matter (SOM). This leads to an increase in the release of carbon dioxide from the soil into the groundwater with consequence pCO 2 increase (Holden et al., 2018).
As a consequence, the most hazardous part of the Critical Zone for contaminant dispersal in these aquifer types is the uppermost highly permeable~10 m below the water table as depicted in the conceptual model in Fig. 11. Reynolds number is up to~1 in this highly permeable   (Allen et al., 1997;Cooper and Lawley, 2007).
At our field site, the hydraulic gradient and hence groundwater velocity does not significantly vary during the hydraulic year (Figs. 7C, 8A; 9). This likely arises from the boundary conditions of the aquifer system at the University of Leeds Farm site, i.e. streams coming from the adjacent hills to the west partially recharge the aquifer, maintaining a stable hydraulic gradient (Aldrick, 1978). These findings indicate how, at least under conditions of relatively minor seasonal hydraulic variation, changes in groundwater flow velocities with depth are much stronger than seasonal fluctuations. In such conditions, steady state groundwater flow models might be adequate to compute solute

Future research
Future research on fracture-flow carbonate aquifer-types is needed in several directions. We envisage efforts on both characterization of hydraulic conductivity versus depth, and flow and contaminant transport modelling. Hydraulic conductivity can be better constrained in open wells using flow logging (Day-Lewis et al., 2011;Sawdey and Reeve, 2012), FLUTe liner profiling (Quinn et al., 2015;Ren et al., 2018a) and nuclear magnetic resonance (NMR) testing (Walsh et al., 2013;Ren et al., 2018b). All these techniques that have been recently developed for fractured aquifers at shallow depths allow permeability characterization along the entire well (rather than specific intervals).
This research has also shown how the limit of the laminar/turbulent transition field can be reached even in non-faulted blocks of carbonate rocks, in absence of karstic cavities. Thus, flow regime is likely to change as we approach areas of more developed carbonate dissolution, such as extensional faults, which are characterized by high transmissivities, or systems of cavities connected to springs as seen in several carbonate aquifers across the world (Galdenzi and Menichetti, 1995;Raeisi and Karami, 1997;Gallegos et al., 2013;Bauer et al., 2016). Groundwater flow models need to be developed that can be implemented in correspondence of such structures, i.e. to include likely occurrence of turbulence and higher groundwater flow velocities (Hill et al., 2010;Saller et al., 2013).

Conclusions
Contaminant transport occurs at relatively high velocity in the saturated parts of shallow fractured aquifers, and this may be very important in areas of arable and livestock farming activity. Here, we report an investigation of the Critical Zone of a fractured dolostone aquifer (Permian dolostone of the Cadeby Formation) beneath the University of Leeds Farm, Yorkshire. A combination of fluid logging, slug tests and hydraulic head monitoring was used. Outputs from these techniques were applied to a parallel plate model to determine the spectrum of flow velocities and establish occurrence of laminar or turbulent flow.
Fluid temperature and conductivity logs identified the bedding plane fractures as principal flow pathways which is a common hydraulic feature in fractured carbonate aquifers. Hydraulic testing (e.g., slug tests) allows determination of the hydraulic fracture aperture by applying the cubic law. Determination of hydraulic gradient from hydraulic head monitoring then allows the determination of the spectrum of groundwater velocities in time and space. Indeed, this workflow was applied at shallow (<~40 mBGL) depths in a non-faulted part of the Permian dolostone of the Cadeby Formation. It produced fracture flow velocities ranging from 13 up to 242 m/day. Depth variations of fracture flow velocities are found much more significant than seasonal variations.
The upper~10 m immediately below the water table were most highly hydraulically conductive (K h = 0.30-2.85 m/day) and characterized by elevated groundwater flow velocities (V = 33-242 m/day) due to the processes of fracture aperture enhancement by carbonate dissolution which are typical of fractured carbonates. Hence, thesẽ 10 m represent the most hazardous part of the Critical Zone for contaminant transport. As a consequence, we envisage further efforts on characterization of permeability with depth in the Critical Zone of fractured carbonate aquifers.
Another key finding is that within unfaulted blocks in absence of large karstic cavities the Reynolds number may reach~1. Thus, turbulent flow is likely to arise as we approach tectonic structures (e.g., normal faults) which are characterized by enhanced carbonate dissolution and macro-karstic features. Groundwater flow models must be developed that represent flow as turbulent in correspondence of such tectonic structures.