Fluid flow from matrix to fractures in Early Jurassic shales

The potential of shale reservoirs for gas extraction is largely determined by the permeability of the rock. Typical pore diameters in shales range from the μ m down to the nm scale. The permeability of shale reservoirs is a function of the interconnectivity between the pore space and the natural fracture network present. We have measured the permeability of the Whitby Mudstone, the exposed counterpart of the Posidonia Shales buried in the Dutch subsurface and a possible target for unconventional gas, using di ﬀ erent methods and established a correlation with the microstructures and pore networks present down to the nanometer scale. Whitby Mudstone is a clay rich rock with a low porosity. The permeability of the Whitby Mudstone is in the range of 10 − 18 m 2 – 10 − 21 m 2 . 2D microstructures of the Whitby Mudstone show no connected pore networks, but isolated pore bodies mainly situated in the clay matrix, whereas 3D data shows that connected pore networks are present in less compacted parts of the rock. A closely spaced interconnected fracture network is often required to speed up transport of ﬂ uids from the matrix into a producing well. For ﬂ uids within the matrix the nearest natural fracture is on average at a distance of approximately 10 cm in the Whitby Mudstone. The combination of the permeability data and the porosity data with natural fracture spacing of the fractures present in outcrops along the Yorkshire coast (UK) resulted in new insights into possible ﬂ uid pathways from reservoir to well.


Introduction
Gas in shales is trapped in poorly connected micro pores by adsorption on and in particles of organic material and clay minerals in the matrix of the host rock. Pore diameters in shales typically range from the μm down to the nm scale, where the majority of the pores have diameters in the sub-μm scale (e.g.: Desbois et al., 2009;Loucks et al., 2009;Keller et al., 2011;Chalmers et al., 2012;Klaver et al., 2012;Houben et al., 2013;Hemes et al., 2013). The transport of fluids through low permeability rocks is controlled by the structure of available transport pathways (Keller et al., 2011). Shale reservoirs and pathways therein are highly heterogeneous and span from the nano-scale to the field-scale (Clarkson et al., 2016). At the nm-cm scales, transport pathways exist of often poorly connected pore networks present in the clay matrix with pores and pore throats often below the μm, whereas at the cm-m scales, an open network of natural fractures control the permeability. Together, the natural fractures and the pore network form flow pathways for fluids through unconventional reservoirs (Curtis, 2002;Gale et al., 2007Gale et al., , 2014Amann-Hildenbrand et al., 2012;Loucks et al., 2012). A closely spaced interconnected network of open fractures is often required to form high permeable pathways for transport of fluids from the tight shale matrix into a producing well (Curtis, 2002;Gale et al., 2014). Characteristic gas production curves from shales show a high initial peak in the production and a rapid decline, before the gas production develops into a sustained low level of production (Hyman et al., 2016). There is no consensus on the exact mechanisms controlling this phenomenon. It has been suggested that the initial peak is related to a fracture network permeability reaching the easy accessible fluids (either stored in the fractures or close to the fractures), whereas the sustained low production level is related to the matrix permeability of the shale (Hyman et al., 2016), where closure of the fractures and pores with time due to decreasing fluid pressures also plays a role.
The research presented here focuses on the permeability of the Whitby Mudstone. The Jurassic (Toarcian) Whitby Mudstone is the UK time equivalent of the European mainland and offshore Posidonia Shale ter Heege et al., 2015), which is one of the main hydrocarbon source rocks (Tissot and Welte, 1978). Interest in the Posidonia Shale and Whitby Mudstone has increased after it was recognized as a possible source for unconventional oil and gas (Herber and de Jager, 2010). This paper shows that the matrix permeability of the rock can be linked to the microstructure of the rock. The fluid pathways present in the Whitby Mudstone are revealed by visualizing the actual pore network in the matrix. The permeability and pore pathway data is combined with natural fracture spacing within the Whitby Mudstone observed in outcrops along the Yorkshire coast (UK). This combination of pore space and fracture network characterization provides new insights into the possible fluid pathways from reservoir to well.

Materials
The Whitby Mudstone consists of the Alum, Mulgrave and Grey shales, where the Mulgrave shale is the lateral equivalent of the Posidonia Shale in Northern Europe (Littke et al., 1991;Powell, 2010;Ghadeer and Macquaker, 2012). The Whitby Mudstone (WMF) samples used in this study originated from outcrops and were collected during fieldwork along the cliff coast North of Whitby (UK) near the villages Runswick Bay and Port Mulgrave (see also Zhubayev et al., 2016). All samples originate from the organic matter rich section of the Whitby mudstone, the bottom eight meters of the Mulgrave Shale member, the Jet rock ( Fig. 1, e.g.: Hesselbo et al., 2000;Powell, 2010;Ghadeer and Macquaker, 2011;Imber et al., 2014;Houben et al., 2016a). Samples were dried in air at room temperature during transportation and storage, and hence dried samples were used in the experiments. Sample blocks taken in the field were used for both permeability measurements (cores and powders) and Scanning Electron Microscope (SEM) imaging. Different subsets were subdivided based on the method and were identified in the following with the letters: S for SEM samples, C for cores used for permeability measurements, and P for powders used for permeability measurements. Different samples within one subset originating from the same hand specimen investigated with the same method are indicated with small letters (a-z) (see also Houben et al., 2016a).

A combination of Precision Ion beam Polishing (PIPS; Precision Ion
Polishing System model 691, Gatan) and Scanning Electron Microscopy (SEM; XL30 SFEG, FEI Company) and a Focused Ion Beam -Scanning Electron Microscope (FIB-SEM; Nova Nanolab 600 Dualbeam/Helios Nanolab UC-G3 Dualbeam, FEI Company) have been used to image microstructures of the Whitby Mudstone samples. Maximum sample diameter of the PIPS-SEM samples was 8 mm in diameter, since 8 mm was the maximum diameter that could be loaded into the Precision Ion Polishing System (Houben et al., 2016a). The samples have been PIPS polished perpendicular to the bedding (Houben et al., 2016a,b). Low resolution mosaics of the entire 8 mm diameter PIPS polished samples were obtained with an SEM, followed by detailed high resolution mosaics of selected areas. SEM imaging was performed using Back Scattered Electron (BSE) and Secondary Electron (SE) imaging modes to image mineralogy as well as the porosity (Houben et al., 2016a). High resolution mosaics were made of 100 single images with a pixel size of 25-50 nm (Houben et al., 2016a,b), and were stitched together with the Microsoft Image Composite Editor (ICE) into one high resolution mosaic covering an area of > 300 × 300 μm 2 (Fig. 2). Image segmentation methods as described in Houben et al. (2016a,b) were used to segment both porosity and mineralogy. Information regarding the 3D distribution and interconnectivity of the pores has been obtained using the FIB-SEM tomography method (De Winter et al., 2009;Liu et al., 2016), and this resulted in a stack of SEM images of consecutive FIBmilled cross sections. Avizo 9 software (FEI, Visualization Sciences Group) has been used to align the image stacks, and FIJI/ImageJ (Abramoff, 2004) has been used to enhance image quality. Both Avizo 9 and FIJI/IMAGEJ were used to perform pore segmentation of the aligned image stacks. For image segmentation a combination of thresholding and marker-based watershed transformation was used in combination with manual correcting of the polygons in Avizo 9. Avizo 9 and FIJI/ImageJ were used for 3D visualization of the pores and quantitative data analysis.

Permeability measurements
3.2.1. Crushed sample test -GRI method (10 − 4 -10 − 3 m scale) using gas expansion pycnometry Permeability measurements have been performed on crushed sam- M.E. Houben et al. International Journal of Coal Geology 175 (2017) 26-39 ples (see also Cui et al., 2009), where the method is generally called the GRI method after the Gas Research Institute. Different subsets from hand specimen WMF 4 were prepared and have been crushed and sieved resulting in three powders with different nominal fragment sizes: ca. 3 mm, 0.5-1 mm, 0.25-0.5 mm. The fragments were oven dried (50°C) prior to the experiment. After drying the sample was weighed and a known mass of the sample was placed in a container (volume V2 GRI ) which was connected to a gas reservoir (volume V1 GRI ) by means of a valve. After both volumes were evacuated for about 1 h the two volumes were disconnected from each other. Ar gas was introduced to V1 GRI and the pressure was increased until the desired pressure had been reached (P1, t1; Fig. 3). Opening the valve between the two ). An enlarged image of the matrix rich areas in two mosaics shows the difference in packing and porosity with the matrix in mosaic 1 having high porosity and loose packing and the matrix in mosaic 3 having low porosity and tighter packing. Below the graph the schematic illustration shows the permeability experimental apparatus used for crushed samples (GRI permeability measurements by gas expansion pycnometry). Ar gas was entered in V1 GRI at t 1 and the pressure was increased until the desired pressure had been reached (P1). At t 2 the valve between the two volumes was opened and that enabled the gas to flow into V2 GRI causing an immediate pressure drop due to filling of the space outside the grains with gas. At t3 the pressure dropped to a new equilibrium pressure, filling V2 GRI minus the volume of the grains, plus the volume of pores inside the grains.
M.E. Houben et al. International Journal of Coal Geology 175 (2017) 26-39 volumes then enabled the gas to flow into V2 GRI and the pressure immediately dropped due to filling of the space in between the fragments (void space) (P2, t2; Fig. 3). The pressure continued to drop until a new equilibrium pressure was reached, filling V2 GRI minus the volume of the fragments, plus the volume of pores inside the fragments (P3, t3; Fig. 3). The methods as described by Cui et al. (2009) and Civan et al. (2011) have been followed to calculate the apparent permeability values, based on the measured pressure decay over time and a number of parameters; porosity, sample radius, gas compressibility, the derivative of the adsorbate density with respect to the gas density, and gas viscosity.
3.2.2. Core samples -pressure step decay (10 − 2 -10 − 1 m scale) The Ar-gas permeametry apparatus at Utrecht University ( Fig. 4; Cui et al., 2009) has been used to measure the Ar-gas permeability of WMF samples, making use of the pressure transient step decay method (e.g.: applied to Rocksalt - Sutherland and Cave, 1980;Peach and Spiers, 1996;applied to Clay -Cui et al., 2009;Ghanizadeh et al., 2013a,b;Rutter et al., 2013). The samples used were 2.5 cm diameter cores, cored perpendicular and parallel to the bedding and the sample length varied between 1.7 cm and 3.3 cm. The sample to be measured was placed in between two plastic cylindrical sealing blocks with central gas-feed pipes and jacketed in a rubber sleeve sealed to the blocks by wire binding. The jacket was pressurized with a jacket pressure (Pj) to prevent undesirable gas flow around the sample either by using Ar gas for the low confining pressure experiments or using silicon oil for the high confining pressure experiments. The Ar-gas permeametry apparatus has an upstream gas reservoir (V1 PS ) and a downstream reservoir (V2 PS ) which are connected together through the sample (V s ) (Fig. 4). All volumes (V1 PS , V2 PS , V s ) were evacuated for at least 1 h before the experiment started to ensure the system only contained Ar gas during the measurements. The jacket pressure and the pressure in V1 PS and V2 PS were increased to the desired values. Low confining pressure experiments were conducted at an effective confining pressure of 0.8-0.9 MPa. For the high confining pressure experiments the effective confining pressure varied during the permeability tests. The system was left to equilibrate for at least 1 h enabling the Ar gas to fully permeate the sample. The upstream (V1 PS ) and downstream (V2 PS ) volumes were disconnected by means of closing valve number 4 (Fig. 4). Afterwards the Ar gas pressure in V2 PS was decreased so that the gas pressure was 0.2 MPa lower than the pressure in V1 PS . The exact pressure and temperature at both sides of the sample were measured every 5 s using two absolute pressure transducers (Keller PR33X, 2 MPa range, precision 0.01%, accuracy 0.05%, 16 bit precision data transfer, and temperature compensation; P1 and P2 in Fig. 4), monitoring the amount of time needed for the gas pressures in V1 PS and V2 PS to equalize by penetrating through the sample. Pressure equilibrium depends on the sample material, jacket pressure, kind of gas that has been used for the flow through and pressures in V1 PS and V2 PS . The raw data was used to calculate the corresponding apparent permeability value using the method originally proposed by Brace et al. (1968). Parameters used to calculate the permeability were: gas viscosity, pressure decay time constant, compressibility of the gas, length of the sample, cross-sectional area of the samples, and the up-and downstream reservoir volumes (Sutherland and Cave, 1980;Peach and Spiers, 1996). After finishing a run the confining pressures and pore pressures were changed to measure the permeability of a core again, measuring permeability at different mean gas pressures was done to correct for the gas slippage effect (Klinkenberg effect;Klinkenberg, 1941). Multiple permeability experiments were performed on different cores consistently showing a decrease in permeability with increasing mean gas pressure, i.e. the measured permeability must be corrected for the Klinkenberg effect (e.g.: Klinkenberg, 1941;Wu et al., 1998;Tanikawa and Shimamoto, 2006).
3.3. Natural fracture network (10 − 3 -10 2 m scale) Natural fractures present in outcropping pavements along the UK coast near Port Mulgrave were imaged on a photo mosaic from five selected interpretation domains. The interpretation domains were in the same location from where the shale samples were collected for porosity and permeability measurements. Imagery was acquired at mm resolution and with approximately 100 m of lateral coverage using an UAV (Unmanned Aerial Vehicle) and with a camera attached to a two meters elevated hand-held stick for local detailed photo collection. The top-down 'bird eye' photos were taken with at least 50% overlap since photogrammetry requires overlap between neighboring images. The processing of the photos was done with AgiSoft PhotoScan software (e.g. Verhoeven, 2011) and a 3D surface was generated from which large mosaics of images were obtained that were orthogonal-rectified onto a 2D horizontal reference surface. The orthogonal-rectification and georeferencing of the photo is based on ground-control-points for which UTM coordinates were measured with GPS and absolute cross-distances with measuring tape. The spatial extent is sufficiently large to depict the largest lineaments up to 10 2 m and the resolution sufficiently high to distinguish lineaments down to a length of 10 − 3 m, covering 5 orders of magnitude. The pavements form a valuable source of information on fracture network arrangements. Some areas of the mosaicked photo are excluded due to occlusions (e.g. boulders, gravel or seaweed) or because of limited quality as shown by the grey areas in the interpreted domain in Fig. 5. Our method distinguishes between a structural acquisition and an analyses phase (Hardebol and Bertotti, 2013). Detailed observations and measurements which required 'handson' access, in search for kinematic indicators, cement infill and compass measurements, were made in the field. Fracture network statistics, such as spacing, were not measured in the field, but were derived from computational geometric analyses that follow detailed digitization of the network geometry with Geographic Information System (GIS) Fig. 4. Schematic illustration of the Ar-gas permeametry apparatus at Utrecht University used for permeability measurements on core samples. The sample to be measured was placed in between two plastic stems and jacketed in a rubber sleeve sealed to the stems. The jacket was pressurized with a jacket pressure (Pj) to prevent leakage flow of gas around the sample either by using Ar gas for the low confining pressure experiments or using silicon oil for the high confining pressure experiments. All volumes (V1 PS , V2 PS , V s ) were evacuated for at least 1 h before the experiment started to ensure the system only contained Ar gas during the measurements. The jacket pressure and the pressure in Volumes 1 PS and 2 PS were increased to the desired values, valve number 4 was opened and the experiment was running until the pressures in V1 PS and V2 PS equalized.
M.E. Houben et al. International Journal of Coal Geology 175 (2017) 26-39 software. The fracture network geometries were captured and analysed using our in-house DigiFract software (Hardebol and Bertotti, 2013). We consider fracture geometries as 2D curvilinear features captured in a 2D quasi-planar slice of a 3D fractured rock volume.

Ion beam polishing and Scanning
Electron Microscopy (10 − 9 -10 − 2 m scale) Results of PIPS-SEM measurements on WMF samples are described in detail in Houben et al. (2016a,b), relevant results for this paper are summarized here. The WMF samples show microstructures with siltsized grains embedded within a fine-grained matrix which is interlayered with organic matter, where matrix is including clay minerals as well as all other minerals with diameters < 2 μm. Samples originating from above the Whale Stones in the Jet rock section show a microstructure dominated by porous carbonate fossils with porosities up to 11% (Houben et al., 2016b). Bedding orientation in the upper half of the Jet rock section is only visible due to the preferred alignment of the clay minerals in the matrix and the clay matrix pores, organic matter and micro cracks that are mostly aligned parallel to the bedding. Samples originating from the bottom three meters of the Jet rock formation show a mm-scale mineralogically layered microstructure with alternating carbonate poorer and carbonate richer layers (Houben et al., 2016a). The high resolution PIPS-SEM mosaics covered areas larger than 300 × 300 μm 2 and are therefore representative for the microstructure present in one layer of the sample, where the Representative Elementary Area for the microstructure is in the order of 200 × 200 μm 2 (Houben et al., 2016b).
SEM porosity in the investigated samples ranges between 0.5% and 2.5%, where pixel sizes in the single SEM images varied between 25 nm and 50 nm (Houben et al., 2016b), hence pores with diameters smaller than 25 nm were not visible. The distribution of the surface area of the pores in the clay matrix ranges between 10 4 nm 2 and 10 8 nm 2 and the distribution follows a power law distribution with a power law exponent of 2.1 ± 0.06 and a constant of proportionality of − 0.86 ± 0.354 (Houben et al., 2016a). Samples WMF4S and WMF15S were selected to be used for FIB-SEM to investigate whether or not there is a visible 3D pore network present in the clay matrix. The FIB-SEM experiments were performed on areas selected in the PIPS-SEM sample overviews (Houben et al., 2016a(Houben et al., , 2016b. Sample WMF4S showed a low porous clay matrix. Sample WMF15S showed layers with a low porous clay matrix similar to sample WMF4S, but also some layers with a less compacted clay matrix with a higher porosity. The PIPS-SEM mosaics of sample WMF4S (Houben et al., 2016a(Houben et al., , 2016b showed a total visible porosity between 1.4% and 1.9% and clay matrix porosity values ranging from 1.1% to 1.4%. A volume of 4 × 2 × 3 μm 3 of sample WMF4S has been analysed by FIB-SEM. The total porosity in the FIB-SEM volume is 1.1% and no inter-connected pore network has been found (Fig. 6). High resolution mosaics were made in three mineralogically different layers of sample WMF15S. The mosaics show total porosity values ranging from 1.4% to 2.5%, while the clay matrix porosity varies between 1.5% and 4.5%. Two volumes were analysed with the FIB-SEM from sample WMF15S. The first volume (7 × 3 × 5 μm 3 ) imaged the clay matrix with a higher porosity and the second volume (9 × 11 × 12 μm 3 ) imaged the pores in the clay matrix with lower porosity. The FIB-SEM data shows a total porosity of 4.8% (low porous clay matrix) and 7.3% (high porous clay matrix). Connected porosity of the low porous clay matrix has a value of 2.1% in M.E. Houben et al. International Journal of Coal Geology 175 (2017) 26-39 the direction perpendicular to the bedding and 0% and 2.2% parallel to the bedding (Fig. 7). For the high porous clay matrix the connected porosity is 6.3% and this is in all directions the same. Comparison of the 2D FIB-SEM pore area size distributions to the 2D PIPS-SEM pore area size distributions shows that the pore size distributions are similar for the same equivalent pore sizes, while the FIB-SEM data shows additional smaller pores since the FIB-SEM pixel size used for imaging was smaller (< 10 nm) than the PIPS-SEM pixel size (> 25 nm) Crushed and powdered samples (WMF4P) with different fragment sizes (Fragment diameters: circa 3 mm, 0.5-1 mm, and 0.25-0.5 mm) were used to measure the permeability using Ar gas. Published WMF porosities for Ar are in the range of 0.1-1.6% for all different nominal fragment sizes, where total porosities for the WMF measured with He are in the order of 1-7% (see Houben et al., 2016a). The density of the gas used in the measurements increases with increasing gas pressure. The skeleton density depends on the applied gas and whether the gas adsorbs onto the samples' pore surface. Ar gas permeability values for all fragment sizes range from 10 − 19 m 2− to 10 − 21 m 2 when intrafragmental porosity values are assumed to range from 1 to 7% for all measured samples. A fragment diameter of 3 mm results in higher permeability values (10 − 19 m 2 -10 − 20 m 2 ) than the samples with smaller fragment diameters. Samples with fragment diameters of 0.5-1 mm and 0.25-0.5 mm show similar permeability values, ranging between 10 − 20 m 2 and 10 − 21 m 2 (Fig. 9).
4.3. Core samples -pressure step decay (10 − 2 -10 − 1 m scale) Permeability measurements were performed on WMF core samples taken throughout the Jet rock section of the Whitby Mudstone Formation. Permeability measurements at low confining pressures (Effective confining pressure: 0.8-0.9 MPa; Fig. 10a) were performed on five drill cores with a 2.5 cm diameter cored perpendicular to the bedding. The Klinkenberg corrected permeability values measured are in the range of 2•10 − 19 m 2 < Ar κ WMF < 1• 10 − 17 m 2 (Fig. 10a). Two samples originating from the same sample block both cored perpendicular to the bedding show similar permeability results (WMF 6; Fig. 10a). Permeability was also measured parallel to the bedding on two sample cores, although both samples originated from the same sample block (WMF 44), they do show about 1.5 orders of magnitude difference in permeability (3 • 10 − 19 and 7• 10 − 18 m 2 ). Furthermore, experiments performed at a range of confining pressures performed on different WMF cores, cored both parallel and perpendicular to the bedding, show that the permeability during the first increase of confining pressures to a maximum value around 40 MPa (pore pressure = 2 MPa, effective confining pressure = 39 MPa) drops over 2-3 orders of magnitude. After, during effective confining pressure cycling, the permeability at different effective confining pressures stays more or less constant (Fig. 10b). In addition, there is one order of magnitude difference between the two measured samples. The sample measured perpendicular to the bedding shows lower permeability (WMF 44Cc; κ = 7• 10 − 21 m 2 at 10 MPa effective confining pressure) than the sample measured parallel to the bedding (WMF400C; κ = 8• 10 − 20 m 2 at 10 MPa effective confining pressure). Permeability values show a difference of four orders of magnitude depending on the effective confining pressure, the sample and the method used.
4.4. Natural fracture network (10 − 3 -10 2 m scale) The natural fracture networks mapped in the multiple subdomains of the pavement display two dominant orientation sets with (i) NeS striking fractures that are typically between 0.5 m and 5 m in length and (ii) E-W striking fractures that are typically < 1 m in length (Fig. 5). The fractures display an orthogonal network geometry with narrowly spaced long NeS systematic fractures and a 2nd set formed by the smaller E-W cross-fractures that abut against the first systematic fracture set (Figs. 5, 11a), see similar fracture networks in Rawnsley et al. (1992) and Rives et al. (1992). In such an orthogonal arrangement, orientation, fracture length and termination or intersection topology are interdependent network properties. For instance, the orthogonal E-W cross-fractures are frequently bound by the larger systematic NeS fractures and the length is thus a function of the spacing between the systematic NeS fractures. Most relevant for this study, which addresses possible flow paths through the tight shale matrix to a nearest fracture, is the fracture spacing and its directional variance (Fig. 11). Maximum travel distance from matrix to fractures is half of the average fracture spacing since the fluid can travel in all directions to the nearest fracture and the average travel distance is half  of the maximum travel distance (Fig. 11b). The spacing between the fractures depends on the fracture orientation and is on average in the order of 25 cm in the E-W direction and up to 55 cm on average in the NeS direction. This implies a shorter travel distance from matrix to fracture in the E-W directions. The mean fluid travel distance from matrix to fracture is one fourth of the fracture spacing and is 0.07 m in the E-W direction and 0.13 m in the NeS direction (Fig. 11c). The NeS fractures were on average longer than the fractures present in the E-W set. Fractures smaller than 0.4 m are dominantly E-W oriented, and fractures with lengths larger than 1.6 m are dominantly NeS oriented. Fractures in the range of 0.4 m to 1.6 m show a mix of E-W/N-S orientations where the NeS orientation is dominant (Fig. 11). About 90% of the matrix is within 30 cm from the nearest fracture in either E, W, N or S directions and the maximum distance of matrix to fracture is about 80 cm in a few areas with minor fracturing. The timing of fracture development could not be established, relative time indicators with respect to burial and exhumation are scarce. One indication that fractures did form before surface exposure comes from veins that are filled with carbonate cement (Imber et al., 2014). Secondly, numerous fractures show abutment relations and linkage with fault zones, which are thought to be formed under peak burial conditions (Imber et al., 2014).

Microstructure and porosity
WMF samples are matrix rich and show typical clay microstructures with silt-sized grains embedded within the fine-grained matrix which is interlayered with organic matter. Pores in the WMF are mainly present in the clay matrix (Houben et al., 2016b). WMF porosity in the 2D PIPS-SEM mosaics is in the order of 0.5%-2.5% (Houben et al., 2016b), which is similar to porosity values reported by Klaver et al. (2012Klaver et al. ( , 2016 and Mathia et al. (2016) for Posidonia Shale on 2D SEM mosaics (respectively 0.8%-2.8% and 1.1%-1.5%). Grathoff et al. (2016) found Posidonia shale FIB-SEM porosities of 1.5% and 2.4% which is on the low end when compared to the WMF FIB-SEM porosity of 1.1%-7.3% reported here. At the micrometer scale, connected pore networks are present in the clay matrix in sample WMF15S, whereas sample WMF4S does not show a connected pore network on the scale used for imaging ( Figs. 6 and 7). In addition, the clay matrix that showed a higher porosity in 2D (WMF15) also shows a higher connected porosity in 3D. The discrepancy between WMF and PSF porosity values and in porosity values of different WMF samples, could be due to a different compaction rate in the different layers/samples. The WMF15S high porous layer shows a higher carbonate content than the other WMF15S layers (Houben et al., 2016b) and the carbonate and silicate minerals are larger in this particular layer, influencing the compaction rate of the Fig. 10. Permeability results of the 2.5 cm diameter core samples. a. Samples measured under low effective confining pressure at different absolute mean gas pressures, cores were cored both parallel and perpendicular to the bedding. b. Samples measured at different effective confining pressures samples were measured both perpendicular and parallel to the bedding. M.E. Houben et al. International Journal of Coal Geology 175 (2017) 26-39 clay matrix causing higher porosity in the clay matrix in this particular layer. The slight difference in porosity values is comparable to a porosity difference found by Keller et al. (2013a) in the compacted clay matrix versus pressure shadow clay matrix in Opalinus Clay. Pore size distributions show that, for pores with pore areas between 10 4 nm 2 and 10 8 nm 2 , the pore sizes follow a power law size distribution that is similar for both the FIB-SEM and the PIPS-SEM samples (Houben et al., 2016a), where the power law exponent is 2.1 ± 0.06 and the constant of proportionality is − 1.9 ± 0.35. Although Houben et al. (2016a) expected that smaller pores follow similar pore size distributions as the ones visible in the PIPS-SEM mosaics, the FIB-SEM data shows that there is a decline in the power law exponent for pores with areas in the range of 10 2 nm 2 and 10 4 nm 2 , whereas the area pore size distribution steepens again for pores with areas smaller than 10 2 nm 2 (Fig. 8). The FIB-SEM size distribution trend is comparable to the combined trend of PIPS-SEM and Ar gas adsorption data (BET-Ar) presented by Houben et al. (2016a). Houben et al. (2016a) calculated associated permeability values that can be expected for the WMF based on porosity-permeability relationships as published by Yagiz (2009), Yang andAplin (2007), and Yang and Aplin (1998). The calculated permeability values are in the order of 10 − 24 m 2 -10 − 18 m 2 (Houben et al., 2016a). Although the range of modelled permeability values overlaps with the range of permeability values measured for the WMF samples, the modelled permeability also shows values that are 3 orders of magnitude lower than the lowest permeability values measured (lowest permeability measured = 10 − 21 m 2 ). One of the reasons could be the underestimation of the slip effect in the modelled permeability values versus the measured permeability values (Afsharpoor and Javadpour, 2016).

Pore network
The PIPS-SEM data shows that pores are not connected in 2D and that a connected pore network, when presented, should be within the matrix (Houben et al., 2016a). FIB-SEM data of pore networks in the Posidonia Shale (Grathoff et al., 2016) did not show a percolating network in the investigated samples when a voxel size of minimal 40 × 40 × 25 nm was used. Our findings for WMF samples show a percolating network of pores in some layers, where the pore connectivity is directionally dependent. The connected pore network could play an important role in gas storage and production. Naturally, lots of connected pores store more gas and make production easier, but the typical slit like pores in the clay matrix are more prone to collapse due to an increase in effective stress during production (Clarkson et al., 2016). Based on low pressure Ar and N 2 adsorption and PIPS-SEM measurements, Houben et al. (2016a) expected that most of the pore connections in the WMF are < 10 nm. Therefore, most of the porosity seems to be unconnected at the imaging scale used for PIPS-SEM. In addition, Keller et al. (2013a) show that due to dilatation or fracturing non-connected pore bodies could easily grow together and that the majority of bridges needed to connect pore tips are on the order of 100-300 nm in Opalinus Clay. There is no evidence found in the PIPS-SEM images that significant volume of isolated unconnected porosity is present if it is assumed that most porosity in the clay matrix is connected below image resolution, as is strongly suggested by the FIB-SEM data. What has to be noted is that the FIB-SEM cubes were aimed to image pore networks in the matrix and that the imaged volumes are not representative for the entire microstructure of the WMF (Houben et al., 2016b), at the FIB-SEM scale used for imaging the porosity is influenced by error that decreases asymptotically with increasing the imaged area (Keller et al., 2013b).
Gas in shales is stored within the microstructure of the rock; adsorbed to the organic matter or clay minerals, or stored in pores or fractures (Clarkson et al., 2016). The investigated connected pore network in the clay matrix together with the non-clay silt-sized grains play an important role in the transport of fluids through the rock because of the enhanced porosity along the boundaries between the matrix and the silt-sized grains (Keller et al., 2013a). For the WMF this is true since the permeability values where highest for the sample with most and the largest silt fragments.

Permeability
The GRI method shows permeability values in the order of 10 − 21 m 2 to 10 − 19 m 2 . The larger the fragment size of the powder measured the higher the permeability. Since the Representative Elementary Area (REA) for the microstructure of the WMF is in the order of 200 × 200 μm 2 (Houben et al., 2016b) all the fragment sizes used in the GRI experiments were large enough to represent the microstructure. Samples do contain larger fractures and macro pores, both of which are easily destroyed during crushing, hence samples with a larger fragment size (and core samples) more likely exhibit more larger fractures and macro pores and hence have a higher chance of a high permeability than samples existing of a smaller fragment size (Cui et al., 2009). The GRI method does not take into account the confining pressure experienced by the rock in the subsurface. Therefore, values are best compared to the low confining pressure core experiments showing a Klinkenberg corrected permeability parallel and perpendicular to the bedding in the range of 1 • 10 − 17 m 2 to 2• 10 − 19 m 2 (confining pressures of 0.8-0.9 MPa). The GRI method clearly shows lower permeability values than the core permeability values measured at low confining pressure. Higher confining pressures (up to 38 MPa) show permeability results down to 10 − 19 m 2 to 10 − 20 m 2 . The discrepancy between low and high confining permeability values could be due to sealing problems where the gas is more likely to have bypassed the sample at lower confining pressures (Sutherland and Cave, 1980;Ramakrishnan and Supp, 2015), but could also be due to the closure of large pore and cracks under higher confinement. GRI permeability values are partially overlapping with the permeability values measured at high confining pressures. Work performed on other shales also show that pulse-decay permeability values are equal to or higher than the permeability values measured with the GRI method (Handwerger et al., 2011;Cui et al., 2009;Heller et al., 2014;Ghanizadeh et al., 2015).
The permeability data shows that permeability does not solely depend on the microstructure of the sample. Controls on the permeability are linked to microstructural complexity of the samples, which is an interplay of porosity and mineralogy, but also the exact arrangement of the minerals (grain boundaries), pores and possible micro fractures at different length scales. The FIB-SEM data shows clearly that sample WMF15 contains a connected pore network down to the nm-scale. On the contrary, sample WMF4 does not contain a connected pore network, confirming that a better connected pore network exists in parts of the sample displaying higher permeability values.
Other organic-rich shales (Heller et al., 2014;Chalmers et al., 2012;Ghanizadeh et al., 2013aGhanizadeh et al., ,b, 2014Rutter et al., 2013) show permeability values in the order of 2• 10 − 17 m 2 -5•10 − 20 m 2 (see Fig. 12a). These values are similar to the permeability values measured on the WMF samples when extrapolated to effective stresses of 0.8-0.9 MPa, where the CH 4 permeability values are comparable to Ar permeability values according to Ghanizadeh et al. (2014). Although the permeability values for other shales fall in the same range as measured for the WMF one has to keep in mind that most of these measurements were performed on samples parallel to the bedding (Heller et al., 2014;Chalmers et al., 2012;Ghanizadeh et al., 2013aGhanizadeh et al., ,b, 2014. Published Posidonia Shale Ar gas permeability values were measured parallel to the bedding at confining pressures ranging from 19 to 38 MPa (Ghanizadeh et al., 2014), extrapolating this data down to a confining pressure of 0.8-0.9 MPa results in permeability values in the range of 5• 10 − 16 m 2 -9•10 − 16 m 2 (Fig. 12b), where these permeability values are higher than what we measured for the WMF at low confining pressures. Rutter et al. (2013) show Ar permeability of WMF samples both parallel and perpendicular to the bedding, where they confirm a permeability dependence on confining pressure also measured in the research presented here (Fig. 13a). They measured a difference in permeability for samples perpendicular and parallel to the bedding slightly larger than the difference that we have measured, where the permeability perpendicular to the bedding is lower than the permeability parallel to the bedding. It is clear from the stress dependence data that the low confining pressure permeability values (0.8-0.9 MPa) could easily drop 1 or 2 orders of magnitude once exposed to higher confining pressures thus resulting in a permeability value on average in the order of 1 •10 − 20 m 2 , which is similar to the GRI permeability values measured. The drastic permeability drop of 1-2 orders of magnitude is probably due to closing of pores and cracks that close progressively with increasing confining pressure, changing the perco-  (Heller et al., 2014;Chalmers et al., 2012;Ghanizadeh et al., 2013aGhanizadeh et al., , 2013bGhanizadeh et al., , 2014. b. Ar gas permeability of the Posidonia shale on samples cored parallel to the bedding as published by Ghanizadeh et al. (2014). Fig. 13. a. Modified after Rutter et al. (2013), the graph shows the dependence of WMF permeability data on effective confining pressure for both samples measured by Rutter et al. (2013) and sample measurements reported on in this paper. Measurements were made both perpendicular (blue colors/squares) and parallel (red colors/circles) to the bedding. b. Modified after Yang and Aplin (2007). Mudstone permeability-porosity dataset, where the different symbols represent the range of clay contents in the samples investigated by Yang and Aplin (2007) and the red (online version) dot represents an average of the WMF data (clay content 60-80%). The permeability has been measured perpendicular to the bedding. c. Mudstone permeability-porosity dataset as measured by Yang and Aplin (2007) measure parallel to the bedding, where the WMF dataset is represented by the pink (online version) ellipse. Clay content for the WMF data is higher than clay contents measured in the Yang and Aplin (2007) dataset. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) M.E. Houben et al. International Journal of Coal Geology 175 (2017) 26-39 lating network hindering fluid flow through the samples.
To check whether a porosity-permeability relation exists in the WMF samples more WMF data is needed, preferably the permeability should be measured on samples with higher porosity. The porositypermeability data in general can be compared to findings of Yang and Aplin (2010) who showed for a large data-set of clays that there is no single and simple relationship between permeability and porosity in mudstones and that the variability is mainly controlled by the clay content. Yang and Aplin (2007) measured a permeability-porosity relationship on cores from the North Sea and the Gulf of Mexico, were they used liquid (30,000 mg/L NaCl solution) instead off gas to measure the permeability. Comparing our data to Yang and Aplin's (2007) measured permeability-porosity data, shows that the Whitby Mudstone has really low porosity compared to the other mudstones investigated and a higher clay content (WMF; 60-80%). The permeability values of the WMF can only be compared to Yang and Aplin's (2007) data when we can compare the gas permeability to a liquid permeability, where we assume that the gas permeability is about 1-2 orders of magnitude higher than liquid permeability (see Busch and Amann-Hildenbrand, 2013). Permeability perpendicular to the bedding is following the same trend as the Yang and Aplin (2007) data (Fig. 13b). Parallel to the bedding the Yang and Aplin (2007) data does not show a clear trend, but the WMF permeability shows similar permeability values showing less dependency of permeability on clay content when permeability is measured parallel to the bedding. Comparing our data to Yang and Aplin's (2010) measured and modelled permeability-porosity data shows that the Whitby Mudstone has a similar permeability when compared to the other mudstones with a 60-80% clay content, again assuming that gas permeability can be converted to liquid permeability by decreasing the permeability with 1-2 orders of magnitude (Busch and Amann-Hildenbrand, 2013).

Fluid transport from matrix to fracture
Permeability measurements at low confining pressure using different samples show that permeability values vary at least one order of magnitude using samples of different fragment sizes (up to mm-scale) and vary more than two orders of magnitude when core samples (cmscale) were used. Schematically this is illustrated in Fig. 14, where matrix transport depends on the pore network present and can either be relatively fast when the network is well connected and more large pores are present or relatively slow when the pore network is poorly connected and pores and pore connections are relatively small. The larger the sample size the higher the chance that a well-connected pore network is present within the volume and that it is represented by an average higher permeability value measured on larger samples (Ghanizadeh et al., 2015). Better connected pore networks are present in parts of the clay matrix as was shown by the FIB-SEM cubes presented here (e.g.: WMF15). A larger number of silt-sized grains floating in the matrix contribute to a better connected pore network because the silt-sized grains limit compaction of the clay matrix (Bobko and Ulm, 2008;Houben et al., 2014). No evidence was found in the PIPS-SEM samples that more porous bands in the form of porous fossils, framboidal pyrite or porous veins exist. Microcracks, mainly parallel to the bedding, are present in the PIPS-SEM samples but it is not clear whether these are due to sampling or sample preparation. It is assumed that these microcracks are not naturally occurring in the field, especially not when rocks are buried and subject to high confining pressures (Houben et al., 2016a,b). Hence, we found no direct evidence that smaller fractures exist than the ones that form the orthogonal network of fractures as captured by our high resolution mapping of the pavements. The most likely and fastest fluid pathways through the microstructure of the WMF are areas with less compacted matrix in between silt-sized grains. In layers where the amount of silt-sized grains is relatively high porosity and connected porosity in the matrix is higher. where the actual permeability depends on the pore network present on the < mm scale, the permeability of the small orthogonal fractures and the permeability of the large fractures. c. Lowest permeability is present in the matrix and the permeability of the system depends on the distance the gas has to travel through the matrix into the small orthogonal fracture network.
M.E. Houben et al. International Journal of Coal Geology 175 (2017) 26-39 During production fluids will reach the first small orthogonal fractures on average after about 10 cm of travel distance through the matrix (Fig. 15a) and will move into the large fracture network on the 10 m scale where after the fluids will reach the borehole (Fig. 14). All scales demonstrate a characteristic permeability value and the only known value at the moment is the matrix permeability up to the centimeter scale (10 − 19 m 2 to 10 − 20 m 2 ). When assuming an average matrix permeability of 10 − 20 m 2 , and an open orthogonal fracture system, the significance of these matrix permeability values for residence time of fluids in the matrix can be calculated assuming Darcy flow. In reality the fractures could be blocked by mineral infill (Imber et al., 2014) or could have limited apertures under effective compressive stresses normal to the fracture walls, blocking direct pathways from the matrix into the well. Assuming open fractures, an average fluid pressure difference of 10 MPa between the matrix and the orthogonal fracture network, and Darcy type flow q = κA μ P x ∆ ∆ [m 3 /s] (where q is the total discharge, κ is the permeability, A is the cross-sectional area of fluid flow, μ is the dynamic viscosity, ΔP is the pressure difference between the matrix and the fractures, and Δx is the average travelling distance from matrix to fracture), the residence time of fluid in the matrix before it enters the nearest fracture can be calculated. Since we know the distance to the nearest fracture (Fig. 15a) from the field we can calculate the residence time of the fluid in the matrix (Fig. 15b, c), following the previous stated assumptions the residence time of fluid in the matrix is in the order of 200 days for 95% of the fluid. As we assume that the matrix permeability is the determining and limiting factor, 95% of the fluids could travel to the well within a year assuming none of the parameters change during fluid extraction from the well, because of the relatively small distances between neighboring fractures. In addition, Fig. 15d shows that residence time in the matrix can change drastically when the matrix permeability is not 10 − 20 m 2 , but higher or lower, and when distance from matrix to nearest fracture is larger. For instance at a Fig. 15. a. This graph shows a fracture network as present in the pavements along the coast near Port Mulgrove (UK), where the colour code tells you whether the matrix is close to the fracture (blue colors) or far from the fractures (red colors). b. Residence time of fluids in the matrix if an average permeability value of 10 − 20 m 2 can be used for the matrix permeability, and an average fluid pressure difference of 10 MPa between the matrix and the orthogonal fractures is assumed together with Darcy type flow (viscosity 1.10 − 5 Pa s). c. Residence time of fluids in the matrix where the normalized fraction of fluids is plot against the days, implying that when the same parameters were used a in b > 90% of the fluids have left the matrix after 200 days. d. Dependency of the speed of flow on the permeability of the matrix. Travel time through the matrix is directly related to the permeability value and the residence time of fluids in the matrix also depends on fracture spacing. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) M.E. Houben et al. International Journal of Coal Geology 175 (2017) 26-39 matrix permeability of 10 − 18 m 2 fluids travel a meter within 10 days, whereas at a matrix permeability of 10 − 21 m 2 one meter of fluid transport through the matrix will take a little bit over 27 years. Time frames of days to several tens of years for gas to travel one meter imply though that the gas could also have travelled out of the formation during geological history depending on when the formation was formed and how thick the overlying package of tight rocks is. Herber and de Jager (2010) comment on this by stating that gas saturation levels are higher in formations more recently uplifted.

Conclusions
The Whitby Mudstone is a matrix rich rock with silt-sized grains floating within the matrix. Porosity of the Whitby Mudstone is < 8% and varies in 2D, where lower 2D porosity values result in lower 3D porosity values and worse connected pore networks in that parts of the rock. The differences in matrix porosity present within the Whitby Mudstone could be due to a different compaction rate in the microstructural different layers. A larger number of silt-sized grains present in the matrix of the Whitby Mudstone caused a higher porosity and a higher connected porosity in that layer, where samples that exhibit these layers also displayed a higher permeability. Permeability does not solely depend on the microstructure of the sample but is linked to the microstructural complexity which is an interplay of porosity and mineralogy, but also the exact arrangement of the silt-sized grains (grain boundaries), pores and possible micro fractures at different scales. Both methods used to measure permeability point towards low matrix permeability values of 1 • 10 − 20 m 2 , where permeability is confining pressure and direction dependent. For our studied high density fracture networks, fluids reach the first small orthogonal fractures on average after about 10 cm of travel distance through the matrix where after they will move into the large systematic fracture network on the 10 m scale. No evidence was found that a network of smaller fractures exist when these rocks are buried in the subsurface. The fastest fluid pathway through the Whitby Mudstone and highest permeability of the Whitby Mudstone are probably found in the areas with less compacted matrix in between silt-sized grains. Using our best estimate for the matrix permeability together with assuming open fractures, an average fluid pressure difference of 10 MPa between the matrix and the orthogonal fractures, and Darcy type flow, the estimated residence time of fluid in the matrix is in the order of less than one year.