Predicting Fluid Flow Regime, Permeability, and Diffusivity in Mudrocks from Multiscale Pore Characterisation

In geoenergy applications, mudrocks prevent fluids to leak from temporary (H2, CH4) or permanent (CO2, radioactive waste) storage/disposal sites and serve as a source and reservoir for unconventional oil and gas. Understanding transport properties integrated with dominant fluid flow mechanisms in mudrocks is essential to better predict the performance of mudrocks within these applications. In this study, small-angle neutron scattering (SANS) experiments were conducted on 71 samples from 13 different sets of mudrocks across the globe to capture the pore structure of nearly the full pore size spectrum (2 nm–5 μm). We develop fractal models to predict transport properties (permeability and diffusivity) based on the SANS-derived pore size distributions. The results indicate that transport phenomena in mudrocks are intrinsically pore size-dependent. Depending on hydrostatic pore pressures, transition flow develops in micropores, slip flow in meso- and macropores, and continuum flow in larger macropores. Fluid flow regimes progress towards larger pore sizes during reservoir depletion or smaller pore sizes during fluid storage, so when pressure is decreased or increased, respectively. Capturing the heterogeneity of mudrocks by considering fractal dimension and tortuosity fractal dimension for defined pore size ranges, fractal models integrate apparent permeability with slip flow, Darcy permeability with continuum flow, and gas diffusivity with diffusion flow in the matrix. This new model of pore size-dependent transport and integrated transport properties using fractal models yields a systematic approach that can also inform multiscale multi-physics models to better understand fluid flow and transport phenomena in mudrocks on the reservoir and basin scale.

Slip factor D f (-) Fractal dimension D b (m 2 /s) Diffusion coefficient of diffusing species in bulk fluid D c (m 2 /s) Gas diffusion coefficient D pd (m) Present-day burial depth D eff (m 2 /s) Effective diffusion coefficients D im (m 2 /s) Intermediate diffusion coefficients D Kn (m 2 /s) Knudsen diffusion coefficients D 0 (m 2 /s) Molecular diffusion coefficients D G (-) General topological dimension D τ (-) Tortuosity fractal dimension dV/dlogD (cm 3 /g) Logarithmic differential pore volume distribution g (m/s 2 ) Gravitational acceleration I; I(Q) (cm −1 ) Scattering intensity K app (m 2 ) Apparent permeability K D (m 2 ) Darcy permeability K exp (m 2 ) Experimental permeability K L (m 2 ) Intrinsic permeability K t (m 2 ) Total matrix permeability K n (-) Knudsen number Kn (-) Average Knudsen number L (nm) Tortuous length of capillary tubes L 0 (nm) Straight length of capillary tubes M (g/mol) Atomic mass of the mixture; gas molecular weight m (-) Slope; power-law exponent N (-) Cumulative number of capillary tubes N A (mol −1 ) Avogadro's number P (Pa) Pressure P p (Pa) Intrinsic pore fluid pressure p (Pa) Mean gas pore pressure Q (Å −1 ) Scattering vector Q J (Kg m 3 /s) Total gas diffusion flux Q t (m 3 /s) Total gas flow rate q g (m 3 /s) Gas slip flow rate q L (m 3 /s) Gas Darcy flow rate q J (Kg m 3 /s) Diffusive gas flux R (J/K/mol) Gas molecular constant T (K) Temperature VR r (%) Vitrinite reflectance Greek Letters (-) Tortuosity exponent ΔC (Kg) Concentration difference ΔP (Pa) Pressure gradient (nm) Mean free path (-) Ratio the straight length of capillary tubes to the average pore size (-) The exponent for the VSS model (-) Viscosity index (-) Intermolecular collision coefficient for the VSS model (Å) Wavelength g (Pa s) Gas viscosity (m/s) Mean molecular velocity (nm 2 ) Squared maximum pore size p (g/cm 3 ) Pore fluid density (-) Tortuosity (-) Average tortuosity (-) Porosity (nm) Pore size or pore diameter max (nm) Maximum pore size ; mean (nm) Mean (average) pore size min (nm) Minimum pore size

Introduction
Technologies utilising the subsurface are impacted by the presence and properties of mudrocks. This includes the energy industry evaluating top seals for hydrocarbons or the properties of shale gas reservoirs, but also applications relating to the energy transition like permanent storage of CO 2 , or intermittent storage of H 2 or CH 4 (Amann-Hildenbrand et al. 2013;Beckingham and Winningham 2020;Ilgen et al. 2017).
In addition, mudrocks have been identified as a potential host rock for the disposal of radioactive waste, where H 2 can be generated from anoxic corrosion of stainless-steel waste containers and from water radiolysis reactions caused by alpha decay (Charlet et al. 2017;Sellin and Leupin 2013). To assess the feasibility of mudrocks for these (geo)technical applications, it is necessary to characterise their pore structures Rutter et al. 2017). The study of porosity in mudrocks has improved through the (combined) application of standard to advanced techniques, such as fluid immersion, gas adsorption, mercury intrusion porosimetry, electron microscopy, nuclear magnetic resonance, or X-ray and neutron scattering (Anovitz and Cole 2015;Busch et al. 2016Busch et al. , 2017Leu et al. 2016).
However, our understanding of how fluid flow regimes and transport properties (e.g. permeability and diffusivity) are controlled by the pore structure in mudrocks across different scales is limited. The pore structure of mudrocks consists of inter-and intra-particle pore space related to organic and inorganic matrix components (Chalmers et al. 2012;Curtis et al. 2010;Loucks et al. 2012;Nelson 2009). Pore sizes generally range over several orders of magnitude, including macropores > 50 nm, mesopores 2-50 nm, and micropores < 2 nm according to the International Union of Pure and Applied Chemistry (IUPAC) pore size classification (Sing et al. 1985). Intrinsic permeability is a function of topology and morphology of pores (Day-Stirrat et al. 2011;Kuila et al. 2014;Loucks et al. 2009), even though the permeability in mudrocks is also stress dependent (Cui et al. 2009). In addition to traditional Hagen-Poiseuille or Darcy-type viscous flow descriptions, slip flow governs transport phenomena in mudrocks that encompass pores from macrometer to micrometer scales (Amann- Hildenbrand et al. 2012;Gensterblum et al. 2015;Ilgen et al. 2017;Sakhaee-Pour and Bryant 2012). It has been shown that transport in mudrocks varies at different characteristic time and length scales (Amann-Hildenbrand et al. 2012;Gensterblum et al. 2015;Ghanizadeh et al. 2014b;Javadpour 2009;Javadpour et al. 2007). In this context, the Knudsen number (K n ), defining the ratio between the molecule mean free path length and the pore size, allows characterising the pore size boundaries for fluid flow regimes (Knudsen 1909). In fact, it relates dominant flow regimes to the corresponding range of pore sizes in the matrix: free molecular/Knudsen diffusion flow (K n > 10), transitional flow (0.1 < K n < 10), slip flow (0.001 < K n < 0.1), and continuum/Darcy flow (K n < 0.001) (Colin 2014;Tartakovsky and Dentz 2019). The pore structure of mudrocks accommodates the rock-fluid interactions controlling transport of elements associated with hierarchical pore morphology (Bahadur et al. 2014;Busch et al. 2017). This leads to a scale dependence of effective permeability, which brings about different fluid flow mechanisms at the corresponding pore size (Amann-Hildenbrand et al. 2012;Mehmani et al. 2013).
In this study, we use a combination of small (SANS) and very small-angle neutron scattering (VSANS) to characterise the pore structure of mudrocks across the full pore size range, covering sizes between 2 nm (2E−9 m) and 5 µm (5E−6 m). We quantitatively apply this method on a wide range of mudrocks, representing reservoirs and seals for various geoenergy applications, such as radioactive waste storage, shale gas production, or subsurface gas storage. The combined SANS/VSANS data provides a multiscale and multi-parameter characterisation, including fractal dimensions, specific surface area (SSA), porosity, and pore size distribution (PSD). In combination with Knudsen numbers for single gas flow, this allows the determination of pore size boundaries for possible fluid flow regimes in mudrocks. Given that transport phenomena are pore size-dependent, we show, for the first time, how different fluid flow regimes are controlled by different pore size ranges at different reservoir depths (corresponding to different temperatures and pressures) and how they are related to total porosity. We further introduce matrix permeabilities and gas diffusivities that are coupled to dominant fluid flow regimes. This provides an advanced understanding of permeability and diffusion in mudrocks to assess caprock leakage or unconventional reservoir production. We expect that this novel systematic experimental-analytical approach will inform multiscale, multi-physics models dealing with hydraulic processes in mudrocks on the reservoir and basin scale.

Samples
Experimental work to characterise the pore structure was carried out on two groups of mudrocks, with the first group consisting of 40 organic lean and the second group of 31 organic-rich samples. The mudrocks studied differ in lithology, mineralogy, age, depositional environment, and burial depth (Table 1).

Mineralogy and Geochemistry
Bulk mineralogical compositions were derived from X-ray diffraction patterns of randomly oriented powder preparates of Opalinus, Carmel, Entrada, Posidonia, Carboniferous, Bossier, Haynesville, Eagle Ford, Newark, and Jordan samples on a Bruker D8 diffractometer using CuK -radiation produced at 40 kV and 40 mA. The mineralogy of Boom samples was obtained from Jacops et al. (2017), mineralogy of Våle shale samples was kindly provided by Norske Shell, Norway. Total organic carbon (TOC) content data were measured on powdered samples with a LECO RC-412 Multiphase Carbon/Hydrogen/Moisture Determinator. Vitrinite reflectance ( VR r ) data were obtained using oil immersion (ne = 1.518) on a Zeiss Axio Imager microscope. Details of the mineralogical and geochemical compositions of the mudrocks and details of the measurements are provided in Supporting Information (S2.1 and S2.2).

Very Small-and Small-Angle Neutron Scattering Experiments
VSANS and SANS experiments at ambient pressure and temperature conditions were conducted at the FRM-II facility at the Heinz Maier-Leibnitz Zentrum (MLZ) in Garching, Germany. Air-dried mudrock samples were cut parallel to bedding, fixed on quartz glass slides, and polished to a thickness of ~ 0.2 mm. SANS measures the scattering intensity I(Q) as a function of momentum transfer Q, the resulting scattering curves contain statistical information that allows pore structure interpretations based on a shape model, e.g. the polydisperse spherical (PDSP) model (Melnichenko 2015). We used the KWS-3 instrument, operated by the Jülich Centre for Neutron Science (JCNS) at MLZ, to obtain very small-angle neutron scattering (VSANS) data of all samples, which detects pore sizes of ca. 5 µm-250 nm. Data at KWS-3 were collected at wavelength of λ = 12.8 Å (with a wavelength distribution of the velocity selector Δλ/λ = 0.2), and a sample-to-detector distance of 9.5 m, covering a Q-range from 0.0024 to 0.00016 Å −1 (Pipich and Fu 2015). SANS data was obtained using the KWS-1 instrument, operated by the JCNS at MLZ, covering pore sizes between 250 and 1 nm. Measurements at KWS-1 were performed using wavelengths λ = 5 and λ = 7 Å with a 10% spread at sample-to-detector distances of 1.2, 7.7, and 19.7 m, covering a Q-range of 0.002-0.35 Å −1 Frielinghaus et al. 2015). Data correction, normalisation, radial averaging, and background subtraction were carried out using the QtiKWS software (Pipich 2006), following standard procedures of the instruments. The data processing and analysis were carried out using our MATSAS software ). The analysis of scattering profile yields fractal dimensions, specific surface area (SSA), porosity, and pore size distribution. Full experimental and analytical information are provided in Supporting Information (S2.3).  Olsen et al. (1996), Fink et al. (2018) Predicting Fluid Flow Regime, Permeability, and Diffusivity… 1 3

Fractal Models
Mudrocks are often characterised by a fractal geometry (Liu and Ostadhassan 2017;Radlinski et al. 1996;Zhang et al. 2017); a detailed discussion is available in Radlinski (2006). The geometry of pores in nature can be described by the single number D f , the fractal dimension, representing pores are self-similar but over a limited pore size range (Teixeira 1988;Wong et al. 1986). The deviation from self-similarity across scales can be due to variations in essential mineral constituents, e.g. clay minerals, see Krohn (1988) for further discussions. Mudrocks are foliated, which is why the sample thickness should be as thin as its constituting clayey lamina to satisfy self-similarity. However, one lamina-thick sample does not provide adequate statistical information for adequate pore structure interpretations. Therefore, we argue a trade-off for the sample thickness (~ 200 µm) must be considered to allow for a sufficient Q-range for measuring the relevant slope (m), thereby D f . Given that fractal dimensions provide morphological information on the surface roughness of pore networks, fractal models can be used to predict the matrix permeability (Miao et al. 2015;Yu and Cheng 2002) as well as diffusivity Liu and Nie 2001;Zheng et al. 2018). Based on the pore structure information obtained from SANS, we developed three fractal models to predict: (i) Darcy permeability for continuum flow, (ii) apparent permeability for slip flow regimes and (iii) effective diffusion coefficient (diffusivity) for diffusional flow regimes in mudrocks. It should be noted that this fractal approach is limited to single gas flow and transport obeying the ideal gas law.

Permeability Fractal Models
The fluid pathways of mudrocks are associated with (micro)-fractures as well as matrixhosted pore bodies and throats associated with inorganic and organic compounds (Ghanizadeh et al. 2014b). Depending on the dominant fluid flow regime, the matrix permeability is subdivided into two main categories: Darcy permeability (no slip flow boundary condition; K n < 0.001) and apparent permeability (slip flow boundary condition; 0.001 < K n < 0.1) (Javadpour 2009). We developed analytical fractal solutions that relate permeability to three pore characteristics including fractal dimension ( D f ), tortuosity fractal dimension (D τ ), and porosity ( ) associated with a dominant pore size which ranges between min and max . The fractal model to predict Darcy permeability is based on the Hagen-Poiseuille equation representing flow in a unit cell consisting of a bundle of tortuous capillary tubes with circular cross-sectional area and a fractal distribution of pore sizes (Yu and Cheng 2002); the relationship between pore size and pore number can be described using the general power scaling law leading to fractal dimension (Katz and Thompson 1985;Mildner and Hall 1986): where N is the cumulative number of capillary tubes ≥ , m = 6 − D f is the slope of the pore density distribution directly obtained from the neutron scattering profile (Radlinski 2006). The fractal dimension ( D f ) is determined from the slope (m), which is defined by the technique. Based on a box counting method in image analysis, Yu and Cheng (2002) defined D f = −m , which varies between 1 and 2. In SANS, D f = 6 − m (Radlinski 2006) and represents a surface fractal. These are bulk objects with a rough surface within a certain range of sizes. The surface area is proportional to r D f , with 2 ≤ D f < 3 in 3D and r is the linear length scale. Since the intensity I(Q) = constant × Q −m and 3 < m ≤ 4 (Radlinski 2006), D f must follow D f = 6 − m to vary between 2 and 3. Yu and Cheng (2002) suggest that the number of capillary tubes between and + d can be derived by differentiating Eq. (1): The negative sign in Eq.
(2) implies that the density of capillary tubes decreases with an increase in pore size, and −dN( ) > 0 (Yu and Cheng 2002). In addition to the pore size of capillary tubes, the tortuous pathways have fractal characteristics. Yu and Cheng (2002) argued that the connection between the tortuous capillary size and its length satisfy the same fractal scaling law; this has been verified for sandstone (Chen and Yao 2017), carbonates (Wang et al. 2019) and shales (Sheng et al. 2016;Zhang et al. 2018). Using a modification by Wheatcraft and Tyler (1988), the quantitative relationship between pore size and pore length within a bundle of capillaries is described as and where L 0 (tortuosity = 1 ) and L( ) (tortuosity > 1 ) are the straight and tortuous lengths of capillary tubes between the start and end points of the fractal path. The range of D is 1 < D < 3; D = 1 corresponds to a straight capillary and D = 3 represents a highly tortuous capillary in 3D (Wheatcraft and Tyler 1988). The average tortuosity ( ) can be described as a transformation relationship between the general topological dimension ( D G ) and D (Wheatcraft and Tyler 1988): where is the ratio of L 0 to the average (mean) pore size ( ). For D G = 3 , the tortuous fractal dimension is thus expressed as where is a function of porosity ( ) and calculated from Xu and Yu (2008): The straight capillary tube is related to the total cross area; Wu and Yu (2007) propose that the total pore area ( A p ) can be obtained from: By substituting Eqs.
(2) in (8), the total cross-sectional area A of a unit cell perpendicular to the flow direction is (Yu and Li 2001). L 0 can be expressed as Under continuum flow conditions (K n < 0.001), the pore size is significantly larger than the mean free path length of gas molecules. This results in the dominance of the molecule-molecule collisions leading to viscous Poiseuille flow. The gas flow rate through a single tortuous capillary, q L ( ) , is given by modifying the well-known Hagen-Poiseuille equation (Cussler 1997): where g is gas viscosity (Pa s) and ΔP is the pressure gradient (Pa). The total gas flow rate ( Q t ) (m 3 s −1 ) can be obtained by integrating the individual flow rate q L ( ) over the entire pore size range for continuum flow in a unit cell: Substituting Eqs. (2), (4), and (11) into Eq. (12), the integration gives We assume that the maximum pore size ( max ) does not exceed the length of a capillary tube ( L 0 ), and is therefore described by a similar fractal scaling law (Yu and Cheng 2002). As a result, the intrinsic permeability for continuum flow ( K D ), [m 2 ] can be expressed according to Darcy's law: is transformed to max ∕L 0 8−D to allow both max and L 0 to follow the same fractal behaviour. According to Yu and Li (2001) Note that min and max are pore size limits of the continuum flow regime.
Moreover, the fractal model to predict apparent gas permeability is based on the gas slip flow rate through a single tortuous capillary q g that can be obtained by correlating the viscous Poiseuille flow ( q L ) and the Knudsen number ranging between 0.001 < K n < 0.1: where f K n is the correlation coefficient (Wang et al. 2019), which can be expressed as (Freeman et al. 2011): Here, is the mean free path [m] from kinetic theory: where p is the mean gas pressure (Pa); R represents the universal gas constant (J/(mol K)); T is the temperature (K), and M is the gas molecular weight (g/mol). The total gas flow rate for tortuous capillaries can be expressed as Substituting Eqs. (2), (4), (11), (15), and (16) into Eq. (18), results in By combining Darcy's law and Eq. (19), the apparent permeability ( K app ), is: where the intrinsic permeability ( K L ) is equivalent to: and the slip factor ( b ) (Pa): Note that min and max are the pore size limits of the slip flow regime, mainly depending on pressure and temperature. The total intrinsic permeability of the entire pore size range can be calculated by combining the intrinsic permeabilities associated with continuum flow and slip flow regimes; K t = K D + K L . K D and K L have similar expressions, but they are not necessarily equal since these are intrinsic for different pore size ranges. The fractal permeability models are valid for SANS-derived fractal dimensions ( D f and D ), only.

Diffusion Fractal Model
Information on diffusional flux or effective diffusion coefficients (D eff ) are crucial to analyse the dissipation of gases in the interconnected pore structure of mudrocks (Amann-Hildenbrand et al. 2012;. Fractal dimensions obtained from SANS data can be utilised to develop a fractal model for the estimation of effective diffusion coefficient, D eff ). The fractal model is based on Fick's law (Fick 1855), and relates the diffusive flux to the gradient of the concentration along the diffusing path. The gas flow rate through a single tortuous capillary, q J ( ) , is given by Fick's law (Zheng et al. 2018): where A( ) = 4 2 is the pore area, ΔC is the concentration difference, and L( ) is the tortuous length of a capillary tube that is obtained by Eq. (4). D c is the gas diffusion coefficient in the porous material, which is expressed as (Ghanbarian et al. 2013): where D b is the diffusion coefficient of diffusing species in bulk fluid (typically water or brine) and is average tortuosity. varies between 1 ≤ ≤ 2 , while 1 indicates smooth and 2 rough pores (Moldrup et al. 2001;Zheng et al. 2018). For a smooth pore system, −1 is termed the pore continuity (Moldrup et al. 2001). Using the Wheatcraft and Tyler (1988) modification, gas diffusion coefficient in a tortuous capillary is thus obtained by where D f is the fractal dimension.
The total gas flux for diffusion through a tortuous bundle of capillaries with the total cross section area A can be expressed as Substituting Eqs. (2), (4), (23), and (25) into Eq. (26), results in The diffusive gas flux across A can be obtained by Fick's law (Crank 1975): where A = L 2 0 . By combining Eqs. (27) and (28), the effective diffusion coefficient ( D eff ) based on the fractal model is Equation (29) represents the effective diffusion coefficient as a function of the diffusion coefficient of diffusing species in bulk fluid ( D b ), fractal dimensions ( D f and D ), and structural parameters including tortuosity ( ), minimum ( min ) and maximum ( max ) pore sizes, and a straight capillary tube (L 0 ) with τ = 1. D is calculated using Eq. (6), from Eq. (7), and L 0 from Eq. (10). If K n ≪ 1 , the molecular diffusion transport mode is advection-diffusion in which D eff = D 0 ≡ ∕3 where D 0 is the coefficient of molecular diffusion defined by the kinetic theory of gases and is the mean molecular velocity (m/s). If Knudsen diffusion is characterised by , where D im is the intermediate diffusion coefficient (Tartakovsky and Dentz 2019). Depending on PSD and K n , D eff can be one or a combination of these diffusion coefficients. The fractal diffusion model is valid for SANS-derived fractal dimensions ( D f and D ), only.

Application of Knudsen Number (K n ) to Transport Phenomena
We used K n to characterise fluid flow regimes in mudrocks. For an ideal gas and an inverse power-law collision model, the Knudsen number is defined as K n = ∕ , where is the mean free path (MFP) of a gas molecule and is the pore size. MFP is obtained by (Colin 2014): where g is gas viscosity (Pa s), R = 8.315 J/mol K being the universal gas constant, T is temperature (K), M is molecular weight (Kg/mol), and P is pressure (Pa). represents intermolecular collisions between gas molecules confined in the pore system. Matsumoto (1991, 1992) introduced the variable soft sphere (VSS) model, which corrects the MFP and the collision rate by expressing the deflection angle taken by the molecule after a collision. Accordingly, the intermolecular collision coefficient is obtained by: in which is the exponent for the VSS model and is the temperature exponent of the coefficient of viscosity (viscosity index) for a given gas. These exponents are available in Bird (1994) for a range of gases (e.g. H 2 , CO 2 , or CH 4 ).
√ 2 Table 2 summarises the pore size boundaries of different fluid flow regimes, and Fig. 1 presents the Knudsen number as a function of depth (pore pressure and temperature) for different pore sizes ranging from macropores to meso-and micropores. According to Rouquerol et al. (2014), the average pore size is 4× Pore Volume SSA . The Knudsen number decreases with depth for a given pore size, and it becomes progressively smaller towards larger pores at constant P-T conditions. Accordingly, at shallow present-day burial depth, organic lean mudrocks are dominated by transitional flow in micro-and smaller mesopores, followed by slip and continuum flow in larger meso-and macropores. Organicrich mudrocks at deeper present-day burial depth accommodate slip flow and continuum flow within small and large pores, respectively. Furthermore, Knudsen numbers calculated for the mean pore size of all mudrocks ( Kn ) show that slip flow is the dominant transport mechanism. This slip flow is taking place in the mesopore range, which has the highest population of pores (Fig. 1).

Darcy/Apparent Permeability and Effective Diffusion Coefficient (K D , K app , and D eff )
Pore characteristics of the individual samples are necessary to obtain K D , K app , and D eff . These include fractal dimensions (D f ), tortuosity fractal dimension (D τ ), pore volumes, porosities, average (mean) pore size, as well as minimum and maximum pore sizes (see Supporting Information S3.2). Table 3 summarises calculated mudrock permeabilities and diffusivities using the fractal models and compares them with experimental values on twin sample plugs obtained from published data. The results consistently show higher plug compared to fractal permeabilities (by about 0.25-1.0 order of magnitude). The Darcy fractal permeabilities appear significantly lower than apparent fractal permeability in Opalinus, Boom, and Våle Shale, whereas the difference between K D and K app is less significant in Carmel, Big Hole, Entrada, and most of the organic-rich mudrocks. D τ values are divided into two separate ranges (organic lean or organic-rich mudrocks). In both separate ranges, K D and K app permeabilities increase with increasing D τ . This finding is invalid when considering all mudrocks. Furthermore, K D is positively correlated to D f in organic lean mudrocks, only. K app shows a weak positive correlation with D f for all mudrocks. The fractal diffusion coefficients are in accordance with experimentally determined D eff values, although the relative deviations of fractal D eff values from experimental findings vary between 0.01 and 0.7 order of magnitude for Opalinus Clay samples using tritiated water (HTO) and between 0.01 and 0.93 order of magnitude for Boom Clay samples using HTO and CH 4 . Experimental findings for Boom Clay vary only slightly , suggesting that the pore structure is rather uniform between samples. Similarly, the diffusion fractal model suggests that Boom Clay samples are composed of a uniform pore structure as D f values differ slightly with values ~ 2.9. The wide range of fractal dimensions (2.0-3.0) obtained by SANS gives good confidence in the approach of using fractal dimensions obtained by SANS. This allows capturing the heterogeneity of the pore structures, which increases with an increase in fractal dimension. A detailed discussion of SANS based fractal model to understand diffusive transport parameter has been provided previously by . The authors showed that model findings match experimental results well and SANS data provide a reliable method to retrieve effective diffusion coefficients. The technique enables measurements at different scales and orientations, thus allowing to understand the relationship of transport properties (porosity, SSA, and PSD apart from D eff ) to other rock properties, such as mineralogy.

Table 2
Pore size boundaries of the fluid flow regime. Note that only one sample is presented for each mudrock but may not represent the entire sample set for that mudrock Intrinsic pore fluid pressure is expressed as P p = p gD pd ; P p is defined as the pressure averaged over the pore area of a representative elementary area (REA) where p is the density of the pore fluid (water, 1000 kg/m 3 ) and g is the acceleration due to gravity (Zhao et al. 1998

Dominant Flow Regimes in Mudrocks
Organic lean mudrocks that are currently at shallow burial depths of < 300 m (Opalinus, Boom, Carmel, and Entrada), correspond to low hydrostatic pore pressures of 2-3 MPa, low temperatures of 287-292 K and large mean free path length of 2.4-5.2 nm. This results in transitional flow within pores less than ~ 30 nm, continuum flow within pores larger than ~ 3 µm, and slip flow within pores between 30 nm and 3 µm ( Fig. 2A). Figure 2A suggests an increasing control of slip flow in meso-and macropores, with a contribution of > 50% of the relative pore volume (Fig. 2B). Within these rocks, transitional flow determines the flow regime in micropores and in a large fraction of the mesopores and contributes to ~ 20-40% of the relative pore volume. Organic-rich mudrocks (gas shales) are subject to greater present-day burial depth of > 1000 m, corresponding to higher pore pressures of 12-40 MPa (assuming hydrostatic conditions) and temperatures of 320-420 K, and lower mean free path lengths of 0.15-0.3 nm. This results in transition flow in micropores (< 2 nm), slip flow in meso-and smaller macropores (~ 2-250 nm), and continuum flow in larger macropores (> 250 nm) (Fig. 2C). Figure 2C suggests that organic-rich mudrocks are rather dominated by slip/continuum flow. Lower pore pressures for Posidonia and Carboniferous have caused slip flow to become more dominant and higher pore pressure resulted in continuum flow dominating > 50% of the relative pore volume in Bossier, Haynesville, Eagle Ford, and Newark Shales (Fig. 2D). Therefore, continuum flow will become increasingly important in smaller pores as the pore pressure increases and slip flow will become increasingly dominant for larger pores during pore pressure decrease which becomes relevant when depleting shale gas reservoirs. Although porosity in mudrocks can be high with values well above 10%, a large fraction of this porosity can be associated with pore sizes where molecule/surface  We assume different solutes diffusing in the aqueous phase. D b of different solutes were calculated from the model developed by Boudreau (1997). Full results are listed in Supporting Information (S3.3). interactions dominate and only diffusion or gas slippage is possible. In addition, pore orientation for high porosity mudrocks might be anisotropic due to the increased clay content, improving horizontal yet limiting vertical flux rates (Dabat et al. 2020). The high specific surface area associated with clay minerals and kerogen allows gases (CO 2 , CH 4 or H 2 ) to form a sorptive layer on the pore surfaces (Rother et al. 2007(Rother et al. , 2014, changing pore throat or pore body sizes. This results in an effective porosity reduction, thus a possible change in pore connectivity during production and/or storage. Therefore, average pore sizes and related distributions are the result of random aspect ratios (pore body/pore throat) over the entire pore size range (Busch et al. 2017). These are important controls on fluid flow, diffusion, and sorption mechanisms in mudrocks (Rezaeyan et al. 2019a, b, c;Seemann et al. 2019).
For the samples studied here, the MFP is close to the average pore size in organic lean but smaller than that in organic-rich mudrocks. As a result, transition flow occupying the micropore domain becomes more important, leading to a higher probability of intermolecular collisions that requires a molecular approach to solve the fluid flow in direct simulation Monte-Carlo and/or Lattice-Boltzmann models (Agarwal et al. 2001). In transition flow, the continuum approach and thermodynamic equilibrium assumptions of the Navier-Stokes equations are no longer valid . The slip boundary condition does not apply due to negligible collisions between molecules and the pore  wall (Li et al. 2011); however, the slip flow model may still partly be used in the transition regime particularly for the organic lean mudrocks with average pore sizes of ~ 30 nm. In slip flow, the thickness of the Knudsen layer that forms in the vicinity of the mudrock pore wall approaches the MFP in the meso-and macropores. This results in gas not being in thermodynamic equilibrium, leading to gas slippage in the interconnected pore structure (Dongari et al. 2011;Zhang et al. 2006). The Navier-Stokes equations remain applicable, provided the boundary conditions are modified in the expression of a velocity slip as well as a temperature jump at the wall of slip domain pore sizes (Colin 2011). In continuum flow, fluid flow is the continuity of temperature and velocity between the fluid and the pore wall in the macropores of mudrocks. Flow is solved by the compressible Navier-Stokes equations, the ideal gas equation of state (thermodynamic equilibrium), and classic boundary conditions in Lattice-Boltzmann models (Bird 1994;Colin 2014). According to above, we argue that transport phenomena in mudrocks are pore size-dependent. This suggests that the multi-aspect interaction between bulk volume flow, sorption and transport mechanisms must be adequately addressed in experimental and numerical investigations. By analysing the pore size distribution, total porosity, and specific surface area in relation to pore orientation, we can improve our understanding of transport phenomena and sorption relationships. Figure 3 illustrates Darcy permeability (K D ) for continuum flow and apparent permeability (K app ) and effective diffusion coefficients (D eff ) for diffusion flow in all mudrocks obtained by the fractal models. We find that the difference between K app and K D decreases with decreasing porosity, which can be related to mean pore sizes which decrease with an increase in present-day burial depth (Fig. 3A). Unlike, analytical or numerical models or laboratory tests, fractal models can define permeabilities for dominant flow regimes, which depend on pore size range, pressure, temperature, and molecular size. The permeability fractal models permit distinguishing between the two dominant fluid flow regimes, continuum and slip flow, if pore characteristics (e.g. , D f , D , etc.) are individually specified for each regime. As such, the incorporation of fractal features with pore size-dependent transport phenomena seems useful to allow for an improved prediction of permeability in mudrocks. We further argue that a unifying fractal model cannot be used for every technique, because the tortuosity fractal dimension ( D ) strongly depends on the general topological dimension ( D G ). The definition of D G is different depending on the method used, e.g. 3 for SANS and 2 for SEM. Moreover, as discussed previously, D f is differently determined, which depends on the specific technique applied. Therefore, it is misleading to inform a non-SANS-derived fractal model with fractal dimensions ( D f and D ) obtained from SANS. As an example, Fig. 3B shows that gas diffusivities obtained from a fractal model that is based on SANS data are in good agreement with experimental data in comparison to a non-SANS-derived fractal model. These fractal models developed, discussed, and applied in this study are the first SANS fractal models derived to predict matrix permeabilities and gas diffusivities in mudrocks.

Transport Properties for Dominant Flow Regimes
Micro-fractures can be one of reasons for the difference between K t and K exp . Sample plugs contain both matrix and fractures, while SANS fractal data represent matrix properties only. Matrix permeabilities on gas shales have been determined experimentally by Fisher et al. (2017) and Fink et al. (2017a) based on the pressure pulse-decay method on crushed samples. The matrix permeability of crushed samples was overestimated from K exp by up to 6 orders of magnitude. There are many reasons for the difference, e.g. errors in calculating permeability from pressure transients, suitability crushed rocks for permeability measurements Fisher et al. (2017). In comparison, the total fractal permeability (K total ) determined in this study varies by 0.1-1 order of magnitude from K exp , only (Table 3). Our results suggest that the tortuosity fractals combined with fractal dimensions capture tortuous pore structure with slit like cross-sectional shapes in mudrocks, allowing for an improved estimate of permeability. Nevertheless, predictions based on fractal models match reasonably well with experimental data for samples having porosities > 0.10, while the match is insignificant for lower porosity samples ( < 0.10). Similar findings have been reported by previous studies (Chen and Yao 2017;Xiao et al. 2014;Zhang et al. 2018;Zheng et al. 2018). We suggest pore size limits be constrained to the length in which fractal criteria are satisfied for the individual flow regime so both D f and D represent the heterogeneity of mudrocks.
Clay type/content and compaction (maximum stress) controls the porosity-permeability relationship in mudrocks. With increasing maximum burial depth, mechanical and chemical compaction result in a porosity reduction (Bjørlykke 2006), leading to a decrease in permeability. In contrast, the abundance of framework grains (mainly quartz and carbonates) can help preserving macropores in the absence of chemical compaction, resulting in increased permeabilities. To exemplify this, we can focus on Entrada and Carmel samples, originating from the same location with depth differences of few tens of metres only. Entrada consists of ~ 60 wt% quartz, dolomite, and feldspar and ~ 30 wt% illite. Carmel consists of ~ 15% wt quartz, dolomite, and feldspar with ~ 80 wt% illite. As a result, the average pore size of Carmel is ~ 5 nm which is significantly lower than for Entrada (~ 7.8 nm), resulting in a significantly lower Darcy permeability due to micro-to-mesopores associated with the illite-rich matrix (Supporting Information, S2.1 and S3.1). Mudrocks however accommodate higher apparent than Darcy permeabilities resulting in a greater total permeability (K t ) since slip flow is commonly associated with macropores as well as part of the mesopores (~ 25 nm-50 nm) (Supporting Information, S3.1). If a large fraction of macroporosity is interconnected throughout meso-to macropores, conductivity for flow increases in the pore network.
Furthermore, the permeability of mudrocks has been experimentally tested (Fink et al. 2017b;Gaus et al. 2019;Ghanizadeh et al. 2014a, b), clearly demonstrating that plug permeability and pore volume decreases with an increase in effective stress. In a uniform system, compaction is considered spatially constant; however, the compressibility of different minerals (e.g. clay versus quartz) can be quite different (Dautriat et al. 2011). Assuming that different pore sizes are associated with different mineralogy regardless of diagenetic history (e.g. clays with smaller, quartz/carbonate with larger pores), we can also speculate that the permeability dependence on stress varies over different pore sizes as well. While the permeability fractal model allows calculation of fluid pressure-dependent gas permeability (Zhang et al. 2017), it cannot reproduce effective stress changes since all SANS measurements were done on unstressed samples. Of the samples tested in this study, we expect significant differences in mechanical properties (especially bulk moduli) and therefore differences in stress relaxation of the pore space when bringing the samples to ambient conditions. This aspect cannot be addressed here and requires future work by potentially determining pore size distributions under applied stress. In contrast, stressed permeabilities conducted on sample plugs can only provide a bulk permeability assuming a certain flow regime that dominates. Constraining the fractal model to certain pore sizes provides pore size-dependent (pressure-dependent) permeabilities that can be integrated with the dominant fluid flow mechanism. Figure 3B shows that the effective diffusion coefficient is positively correlated with porosity. Organic-rich mudrocks are characterised by low effective diffusion coefficients Fig. 3 Pore size-dependent transport properties using fractal models. A Matrix permeabilities of the mudrocks calculated for continuum and slip flows (K D and K app , respectively) as well as plug permeabilities (K exp ), the insert plot shows permeabilities for individual samples with porosities ranging between 0.02 and 0.12; B effective diffusion coefficients (D eff ), the diffusion fractal model is compared with the Liu and Nie (2001) and modified Liu and Nie (2001) fractal model (Rezaeyan 2021); C methane diffusion coefficient versus permeabilities (K D and K app ). Non-linear curve fits are obtained using power functions with values on the order of ~ 1E−11 m 2 s −1 ; high clay/kerogen content along with relatively higher present-day overburden stresses result in lower diffusion coefficients. For permeability, the pore throat diameter is the determining factor. In diffusion, tortuosity is a key control, which is again controlled by pore throat and pore body sizes that can change upon changes in effective stress and diagenesis (Fathi and Akkutlu 2014). Yet, this does not invalidate the relationship of low permeability relating to low diffusivity and vice versa, as can be seen in Fig. 3C. Especially for high permeability and high diffusivity samples, a linear relation can be observed, indicating that pore throats are the key controls for both transport modes. We can assume that pore throat sizes are similar to the MFP and as such, concentration driven gas (e.g. CH 4 , CO 2 , H 2 ) diffusion is likely to control migration through these pores with pore throats (Amann-Hildenbrand et al. 2012;Gensterblum et al. 2015;Jarvie et al. 2007;Loucks et al. 2009;Ross and Bustin 2008). Therefore, we can argue that the pressure-driven volume flow and molecular diffusion tend to become distinguishable in extremely low permeability rocks by segregating dominant flow regimes based on effective pore sizes.

Conclusions
Small-angle neutron scattering resolves a wide range of mudrock pore sizes (2.5 nm-5 µm). Fluid flow regimes in mudrocks vary depending on the pore sizes as well as pressure and temperature conditions. For some of the organic lean mudrocks studied, originating from low hydrostatic pore pressures (2-3 MPa) due to shallow depth, gas molecules develop transitional flow within micropores and mesopores with sizes up to 30 nm, slip flow in pore sizes between 30 nm-3 µm, and continuum flow within pores > 3 µm. Most organic-rich mudrocks studied originate from depth associated with high pore fluid pressures (12-40 MPa). Because of the smaller mean free path length at larger depths, continuum flow is dominant in macropores > ~ 250 nm, slip flow in smaller macropores and mesopores, and transitional flow in micropores. With a reduction in pressure during reservoir depletion in gas shale reservoirs, slip flow becomes more dominant for larger pore. In contrast, when injecting gases into the subsurface and pressure is continuously increasing, continuum flow becomes increasingly dominant when gas is flowing through tight mudrocks. This shows that bulk volume flow related to pore pressure changes and pore size distributions needs to be addressed in experimental and numerical investigations. Further complexity relates to diffusion and sorption to understand bulk fluid migration in pore systems of mudrocks. By analysing the pore size distribution, total porosity and SSA in relation to pore orientation, we can inform fluid dynamic models to improve our understanding of these flow-diffusion-sorption relationships.
The study of gas transport in low permeability rocks revolves not only around the validity of dominant fluid flow regimes associated with different pore size ranges, but also their pore size-dependent transport properties. Fractal models calculate Darcy permeability for continuum flow and apparent permeability and effective diffusion coefficients for slip/diffusional flow for the relevant pore sizes in mudrocks. Most mudrocks are characterised by higher apparent permeabilities than Darcy permeability, since slip flow dominates a wide pore size range of ~ 25-250 nm with large pore volumes of up to 70%. If a large fraction of macroporosity is interconnected by meso-and macropores, this results in a higher conductivity to flow for the entire pore network. On the other hand, the increased nanoporosity with small pore throats results in high diffusivity. The pressure-driven volume flow and molecular diffusion tend to become distinguishable in low permeability rocks by segregating dominant flow regimes based on the effective pore sizes.