Interactions Between Imbibition and Pressure-Driven Flow in a Microporous Deformed Limestone

Neutron imaging is used for direct observation of evolving water–air and deuterated water–normal water exchanges in flow experiments performed on a laboratory-deformed, microporous laminated limestone, an extremely fine-textured rock altered by arrays of superposed fractures generated in a rock mechanics apparatus. The neutron images document significant, evolving, water speed and flow direction variability at the deci-micron scale and spatially complex patterns of both increasing and decreasing water saturation. We infer that capillarity-driven and pressure-driven water movement occurs concurrently, in close proximity and in competition, and that as local and global water saturations evolve these two drivers can change their dominance in both matrix and deformed elements. Thin sections are used to obtain sub-micron resolution SEM images that provide multi-scale information on the textural features’ spatial arrangements. The textural characteristics are consistent with the inferences made from the coarser flow imaging. Alternating lamina types provide the primary lithological heterogeneity, while the experimentally created deformations lead to quasi-planar zones of highly comminuted matrix and fracture-like voids, each with lengths ranging from sub-mm to cm. Together deformation features delineate a partially connected array. The interplay between fluid movement through deformation features, and flow into (and out of) the laminae, implies near-equivalence of local driving pressure- and capillary-related energies, with subtle shifts in this balance as water saturation increases. The insights gained invite a re-examination of common rules-of-thumb for multi-phase fluid flow often adopted in fractured, low-permeability microporous rocks. Neutron tomography enables texture and water transit visualisation: SEM provides associated microtexture characterisation Neutron tomography documents water movement implying concurrent imbibition and viscous flow in deformed carbonate rocks Spatio-temporal fluid motion analysis necessitates active bi-directional fluid exchanges between matrix and fractures Neutron tomography enables texture and water transit visualisation: SEM provides associated microtexture characterisation Neutron tomography documents water movement implying concurrent imbibition and viscous flow in deformed carbonate rocks Spatio-temporal fluid motion analysis necessitates active bi-directional fluid exchanges between matrix and fractures


Introduction
Porous-medium characteristics strongly affect movement of liquids and gases through a material's connected void networks. Capillary and viscous effects are widely understood to be the primary factors, along with gravitational effects related to fluid densities (Miller and Miller 1956;Blunt et al. 1992;Ide et al. 2007;Schlüter et al. 2016). In water-air fluid systems, the notion of a capillary fringe is often invoked (Ronen et al. 1997;Nielsen and Perrochet 2000), in which a locally 1D water saturation variation occurs. Such concepts, plus experimental observations from simple rock types (Prodanović et al. 2007), lead to a common view in which interaction of the wetting phase and the non-wetting phase yields a piston-like fluid interface displacement.
Here we describe experimental observations of fluid movements, using neutron imaging to detect spatial and temporal changes, in air-water and in water-water fluid systems, in a rock sample that exhibits significant spatial textural variation. The cylindrical sample is 38 mm diameter, with thin (± 1 mm), alternating layers (laminae) that have differing pore systems. This primary lithological arrangement is overprinted by deformation features produced during a rock-mechanics experiment. The locally tabular new features are nominally open, somewhat like fractures and somewhat like deformation bands (Aydin 1978), as revealed in post-experiment thin sections. Together, these textural types create a sample with multiple, space-delimited pore-system types. The fluid flow experiments provide evidence, consisting of both multi-phase and miscible processes, that allows us to infer the small-scale capillary and viscous properties of the textural components of the rock, and the assembly of these elements into larger systems.
The paper is organised as follows. The rock sample and its geological context and the experimental deformation are described and quantified. Next Scanning Electron Microscopy (SEM) images of selected sample slices are provided. Operationally, acquiring SEM images is the last step because the sample is cut to make thin sections, but it is logical to quantify the rock textures before assessing the fluid flow beam experiments. The SEM section is followed by the design and implementation of the fluid flow experiment performed in the neutron beam. The resulting fluid movement radiographs and tomographs are presented, together with previously acquired X-ray tomograph (XRT) images. The results are then analysed, primarily in terms of continuum-scale flow partitioning between matrix and fracture, with emphasis on the differences between how D 2 O liquid replaces air and how H 2 O liquid replaces the very similar (essentially miscible) D 2 O liquid. Finally, we discuss how these observations relate to the competing driving mechanisms of imbibitionand pressure-driven flow, particularly in rocks with the textural complexity present in this sample.

The Rock Sample
The limestone used here consists predominantly of 5-10 μm-sized particles, and the pore system sizes vary from a few μm to sub-μm diameters. When deformed in the rock mechanics lab to represent natural subsurface strains, this sample is distinctive in that it develops complex fracture arrays that are superposed on, are significantly governed by, and alter, the pre-existing texture and its pore networks. Such complex arrays in few centimetre diameter samples are uncommon in rock mechanics experiments, with most porous rocks of this size developing a single primary discontinuity. This unusual pattern of deformation features makes this rock particularly informative for studying the impact of multiple kinds of textural complexity on fluid flow in samples small enough for tomographic experiments in the neutron beam.
The limestone is a laminite from the Crato Formation that crops out in the Araripe Basin, northeast Brazil (Assine 2007;Neumann, 1999). This sample was collected from Nova Olinda Quarry. The Crato Formation is attributed to a shallow lake/lagoonal depositional setting (Catto et al. 2016). Within this setting, micro-crystalline calcium carbonate crystals were precipitated and locally transported (and sorted) via traction currents and gravitational slumping which, together with diagenetic alterations, created the observed primary textural complexity. The present-day micrite (micro-crystalline calcite) has a dominant particle size of a few μm and creates local, lamina-scale, porosities that range from a few percent to tens of percent.
The mechanical response of the Crato laminites has been assessed by laboratory experiments . Undeformed laminite blocks were cored and these cylindrical samples were subjected to what is normally called a triaxial compression test, i.e., axial shortening in a pressure vessel under constant radial stress (a uniform confining pressure). The samples develop small amounts of permanent strain through the growth of, and small displacement on, fractures and deformation bands . The full suite of experimentally deformed samples from these experiments exhibits arrays of small-scale fractures that show very good geometric correspondence to more widely distributed features (of larger size) that are observed in the Araripe Basin outcrops Miranda et al. 2018). This study focusses on just one sample which has the most complete set of post-deformation information.
The sample's laminae consist predominantly of aggregates of (variably) calcitecemented calcite crystals in 200 μm to 2-mm-thick regular and less regular layers that exhibit contrasting textures, and which are comparable to those described by Catto et al. (2016). When deformed under 10-40 MPa confining pressures with laminae normal to the axial load, the suite of samples develop arrays of short semi-planar fractures . These fractures tend to occur in clusters with a local common orientation so that, when seen in aggregate, they approximate a series of semi-connected, semi-planar, fractures. Typically, they are steeply inclined to the laminae. On closer visual inspection, the arrays are not obviously fully connected, suggesting that sample along-axis fluid movement may well involve interactions between the laminae pairs and the fractures.
Whole-sample porosities of undeformed laminite samples studied by Zihms et al. (2017) vary from 12 to 20%, typically with micro-Darcy bulk permeabilities, in contrast to a few milli-Darcies to Darcies for the bulk permeability of the deformed samples. Experimentally induced fractures substantially increase bulk permeability in some samples, and only raise it slightly in others. The variability will relate to (at least) the individual-sample connectivity of the induced fracture arrays. The laminae are divided into to two classes using SEM images. Figure 1a shows ~ 5-10 μm crystal size, heavily calcite cemented so the residual pore spaces are dominated by slit-like voids. Figure 1b shows more variable crystal sizes, including a common ~ 1-2 μm fraction and less cementation.
The L20 sample deformed at 20MPa confining pressure (Fig. 2) exhibits a complex array of newly generated, mostly short, semi-planar features visible on its surface (Fig. 3). Some new features are visibly open at benchtop conditions while others appear to be closed. Most are steeply inclined to the laminae, and in the lower part they are arranged in a quasi-radial pattern. A thin band, positioned approximately lower-middle in the sample height, shows a few clear open fractures, and indications of semi-continuous openings parallel to the laminae are visible using a hand lens. Above this band the open fractures appear to be fewer, steep and with a small included angle between conjugates.
The visual interpretation of the fracture array is reinforced by XRT (Fig. 4) acquired after the rock mechanics experiments but before the fluid flow experiments in the neutron beam. These data were acquired with a Zeiss XRM520 at the 4D Imaging Lab, Lund University, with a cubic voxel dimension 40 μm. The images in Fig. 4 have been grey-scale thresholded to highlight lower-density voxels in the reconstructed volume which reveals semi-planar lower-density zones. They are interpreted as open fractures containing fragment aggregates which lead to the locally reduced but variable densities. Figure 4 confirms the complexity of the open fracture system inferred from visual observation, including its semi-radial pattern in the lower part of the sample and a moderately well-developed set of steep fractures. The XRT images also support the visual observations that the ~ 5-mm-thick band in the lower-mid part of the sample consists of high-porosity laminae with a reduced occurrence of open fractures. In the upper part of the sample, the XRT images confirms the lack of a radial fracture array, and the set of steeply dipping fractures with a small included angle. Note that some white-filled fractures, or portions of fractures seen by eye on the sample exterior are not seen in the thresholded XRT images. This supports the initial supposition that they are closed shear fractures, or at least that they have minimal dilation accompanying shear offset. The details of the textural effects of these induced features are assessed in Sect. 3.

Detailed Examination of Rock Texture Using SEM
The selected images are all taken from SEM scans of thin sections cut from the L20 sample and from the immediately adjacent undeformed part of the source block. The SEM acquisition uses an FEI Quanta 650 FEG SEM. The settings are low-vacuum, 20 kV, spot size 3.5-4.5 and a working distance approximately 10 mm. The MAPS montage method is used, which acquires slightly overlapping, highest-possible resolution 'tiles' over the extent of the entire thin section (Buckman et al. 2017). This scheme assists in the task of relating any detailed feature to its local and global situation.

Undeformed Rock
A thin section was prepared from the block from which L20 was cored. Figure 5 shows a composite SEM image of the whole undeformed rock thin section, approx. 35 mm by 60 mm. At this composite resolution many tens of laminae are just discernible (Fig. 5) though their texture and the nature of the transition between adjacent laminae is not clear. The more obvious pattern is laminae groups varying from about 3 mm to 20 mm thick. Starting from the base these are divided into the four zones as shown in Fig. 5. See also Fig. 6 and Supplementary Materials (SM)1. The Base Zone and the Upper Zone are very similar, both having the alternating laminae types (Fig. 1). The Upper Middle (UM) Zone has wider, more diffuse and inclined laminae that are otherwise texturally similar to other parts of the sample. The Lower Middle Zone looks like a diagenetically modified version the Base Zone. In this study the focus is on the textural characteristics of the laminae at the sub-μm scale upwards. Figure 6 illustrates the dominant calcite crystal texture and the differing cementation levels. See also SM1. There is little evidence of natural damage to the rock. Figure 6a, b shows the rock texture to be almost exclusively crystalline, with angular crystal boundaries, poor expression of calcite crystallographic shapes and highly variable degrees of cementation. Locally porosity can vary from 3 to 26% with rare areas greater than 50% (Buckman et al. 2017). Typical sizes are -10 μm and pore shape varies from equant to elongate along parallel crystal boundaries (Fig. 6b, d). Significantly, sub-μm pore throats are common. Figure 6a shows the character of each of the two laminae types each  Figure 6c shows a rare highly porous zone with the morphological character of a natural open fracture. However Fig. 6d shows little evidence of the rock damage expected in a fracture, either natural or induced.

Deformed Sample L20
Four polished thin sections were made from L20, chosen for the best oriented cuts through the fracture network seen in XRT. For the purpose of this paper, the individual polished thin sections are treated together. One full sample length composite image SEM is shown in Fig. 7. Also see SM5. The rock and its deformation features are described from bottom to top. This order is followed in the remainder of the paper.

The Base Zone
The Base Zone is the lowest ~ 1 cm of L20 (Fig. 7). Laminae are visible as very low dip alternating lighter and darker bands with indistinct boundaries, plus some lamina-parallel elongate darker regions. There are also fractures, deformation bands or possibly a combination of both features. Figure 8 shows representative sub-regions, moving from left to right with some distributed darker lamina-parallel zones not obvious in Fig. 5 and semivertical fractures. See also SM1. These fractures potentially connect more porous zones within, or bounding, individual laminae via what appears to be several openings along wavy lamina boundaries. It is unclear if the fracturing and the separation along lamina boundaries was induced by triaxial testing or by subsequent sample handling. While the steeper fractures are consistent with known experimentally caused fracturing, laminite samples have a strong tendency to split along lamina boundaries at all stages of handing and such unintended splitting may well explain the laminae-parallel openings. Figure 8b shows more detail of the darker zones seen in Fig. 8a and reveals them to be of similar  Fig. 8c shows a fracture array that was also observed by Zihms et al. (2017). Figure 8d shows some fractures at a much higher magnification. Typical fracture widths are between 10 and 30 μm, equivalent to a three to five crystal diameters. At this magnification a few very small-scale examples of fracture-associated damage within individual crystals can be seen, the fractures themselves being developed primarily along crystal boundaries.  Figure 9 shows example SEM images. Also see SM1. Figure 9a illustrates the abundant steep dominantly shear-offset deformation bands common in parts of the LM Zone, plus one dark lamina-parallel unit. Figure 9b-d shows these features at progressively higher magnifications. At least some LM Zone laminae are highly porous, a result of significant dolomitisation, and they have been offset by shear movement on the deformation bands. Some bands are filled with much finer comminuted material. The rock mass outside the bands has a very similar crystalline character to the Base Zone. Figure 9b also confirms that the darker zone in (a) is considerably more porous: Buckman et al. (2017) calculate 19-26% porosity in certain laminae in this zone, in stark contrast to the rest of sample L20 where porosities above 20% are unusual and the most common porosity is 12-14%. Figure 9b, c confirms that the deformation bands are typically a few microns (μm) wide, have shear offset and contain highly comminuted material. Figure 9c also confirms that the crystal-crystal boundaries and some crystals between the bands are slightly damaged, though any openings are of significantly sub-μm width. Most deformation bands are a few mm long and both right-dipping and left-dipping, in roughly equal numbers. The result is a very complex network of shallow lamina-parallel high porosity zones and steeply inclined deformation bands with variable laboratory-induced void networks.

The UM and Upper Zones
The UM and Upper Zones make up more than half the sample. The laminae in the UM Zone approach 2 mm thick, thicker than in the Base Zone with less distinction in character between adjacent laminae. They are also inclined (dipping) at about 20° from horizontal most probably from sedimentary slumping. The Upper Zone has a very similar undeformed textural characteristics to the Base Zone. Both open fractures and deformation bands are developed (Fig. 10).
The SEM images of fractures and deformation bands that accommodate part of this flow system are shown in Fig. 10. Also see SM1. Most are dip steeply ( Fig. 10a-g) and . Porous zone is more wavy than in undeformed laminite: its offset by deformation bands filled with much finer grained material is clear; c detail of porous zone with diffuse boundary and floating particles (as seen in 2D); d detail of (c) confirms damage at crystal-crystal boundaries is common are a few mm to a few cm long. They tend to terminate at lamina boundaries and enechelon arrangements with partial connectivity is common (Fig. 10a). Some features have boundary-normal (mode 1) opening ( Fig. 10c, d, e) while most display shear offset ( Fig. 10c, e, f), with or without opening. Some features have virtually no contained material while others are full of highly comminuted material. Yet others have mixed fill types. The features that most resemble deformation bands typically show somewhat diffuse damage in the surrounding rock. A few features (Fig. 10e, f, g) show deformed textures very like well-developed en-echelon fault zones but occur on only one side of the open fracture. These features have also been noted by Lewis et al. (2019) for the same sample.

Neutron Beam Flow Experiment
The complex spatial arrangements of L20 voids, comprising the natural pore system and induced elements, with their expected differences in textural characteristics, provide an excellent opportunity to investigate how such variability of material textures impacts fluid-flow. This sample especially presents the possibility to observe fracture-matrix fluid interchanges. Based on published studies observing fluid movement through porous rocks via neutron imaging (e.g., Perfect et al. 2018;Tengattini et al. 2021), plus anecdotal information, we choose this method. The option to use X-ray imaging was not selected, due to challenges associated with distinguishing fluids of comparable densities, where dopants are often employed (Singh et al. 2017a, b).

Selection of Neutron Imaging
While neutron radiography has been used to image fluid flow in rocks since the 1980s (Jasti et al. 1987), the comparatively low availability of neutron sources and the reduced flux, as compared to a synchrotron collinear X-ray beam, have generally limited the use of neutron imaging to radiography or to tomography in quasi-static conditions. Viggiani and Tengattini (2019) and Tengattini et al. (2021) summarise more recent uses of neutron tomography, including its deployment in combination with XRT. The selection of neutron imaging for this sample is also influenced by the outcomes of flow experiments in similar materials (Hall 2013;Tudisco et al. 2019).
Radiographic studies, using X-rays or neutrons, have been judged suitable for assessing flow processes where their projection onto a single plane normal to the beam's line of flight doesn't obscure the observations. This is particularly true when fluid moves with time. The longer times typically required to acquire tomographic data have also restricted the study of inherently very rapid processes because the spatial distribution of the features of interest can change significantly during the acquisition leading to significant blurring in the derived tomographic images. However, recent developments have allowed significantly faster acquisition of data to create neutron tomographs (NTs), as reviewed in Tengattini et al. (2021). For example, rapid neutron tomography has been performed to study the evolving fluid front in rock samples, both for quite homogeneous rocks and for those containing clear deformation features, such as a shear band (e.g., Tudisco et al. 2019). More complex fluid flow patterns have also been studied using both radiographic and post-deformation tomograph acquisitions for a number of porous materials, including shales (Stavropoulou et al. 2019;Roshankhah et al. 2018) and concrete (Kanematsu et al. 2009).
This current work extends the capabilities by acquiring images of complex fluid flow interactions with this sample. We acquired the images used in this paper at the NeXT (D50) beamline at the Institut Max von Laue-Paul Langevin (ILL) (Tengattini et al. 2021). A uniquely high neutron flux is available at NeXT in the form of a neutron beam passing through a 30 mm aperture placed 10 m ahead of the sample position. Beyond the sample, a scintillator-based detector provides 2D images (projections) of the neutrons that pass through the sample and, by subtraction from a no-sample image, the neutron attenuation due to the sample is calculated.
Unlike X-ray methods, in which the beam of photons is primarily affected by the atomic number of the material it interacts with, neutrons exhibit strong interactions with some light elements and are isotope sensitive. They have little interaction with many other materials including most metals and the mineral components of rocks (Tengattini et al. 2021). The strong interaction of neutrons with hydrogen, which scatters a large fraction of neutrons, making water highly visible due to the attenuation, is used. The experiment is designed to obtain images of the entry into and water movement through the originally air-dry rock. Since air has very low attenuation to neutrons, we can readily distinguish water movement through the sample as it replaces air. Deuterium (D), an isotope of hydrogen with an extra neutron in its nucleus, has a much smaller interaction with neutrons and so much less scattering, than does hydrogen (H). Therefore, H 2 O and D 2 O (heavy water) are readily distinguishable in neutron imaging while being chemically and physically very similar. This provides an opportunity to track movement of the interface (or the gradient) between two fluids that can otherwise be considered to form one single fluid phase. D 2 O has lower neutron attenuation than H 2 O, but it is still readily distinguishable from air in neutron tomography images. Because of the highly variable void size distribution in this multiply fractured and laminated sample, calculating a simple fluid front and a simple fluid velocity was recognised as unlikely. But the emergent patterns were expected to inform on the roles played by undeformed and deformed rock volumes, especially when combined with the detailed, measurable SEMderived textures.
Comparable with the more-familiar XRT, high-speed neutron tomography (e.g., Tudisco et al. 2019) acquires multiple 2D projections while a sample is rotated about its central axis in the neutron beam. Thereafter, a 3D reconstruction of the sample is derived from the information in the multiple projections. 3D representations are ideal for capturing the spatial distributions of water, or of the interface between two forms of water, and are used extensively in this work. However, under certain circumstances, 2D radiographs of a nonrotating sample can be almost as useful, in particular because they are acquired more rapidly and can be interpreted in real time by the operator. Using selected 2D radiographic projections, for example when a fracture is aligned with the flight direction of the neutrons, the composite attenuation effect can enhance the recognition of the water contained with the planar feature in contrast to the composite effect of water in the surrounding porous matrix. However, effective radiographic imaging is limited to special such cases and, for example, two planar features with different orientations would require 3D imaging (tomography) to identify their spatial relationships.

Experimental Setup
The fluid flow experiment on L20 was performed in the neutron beam while acquiring in-operando neutron radiography almost continuously, plus NTs at key points. The experimental configuration for this flow experiment is similar to the general description in Tudisco et al. (2019). The cylindrical sample is placed on a basal end-cup to provide good contact between the sample base and a small distribution chamber containing water supplied by the inlet tubing. This set-up also allows the sample assembly to be attached to the rotation stage. The sample was wrapped with PTFE tape, then a heat-shrink poly-fluoride (FEP) tube. These materials were chosen for their reduced hydrogen content. A cable-tie pressed the heat-shrink tube against an O-ring sitting in a groove in the basal end-cup. This method has proved practical for operating unconfined experiments with fluid injection pressures up to 80-100 kPa: approximately 10 kPa was used here. Fluids were injected using a constant head into the distribution chamber in the bottom end-cup and fluids exiting the sample were captured at the top.
Observations of D 2 O replacing air, and then H 2 O replacing D 2 O, were obtained. The neutron beam configuration, the sample size and position and the detector characteristics resulted in pixel images that each represent ~ 20 μm of the sample cross-sectional area: the sample is approximately 38 mm diameter. The fluid flow experiment involves injecting D 2 O into the sample base at a constant head and observing its transit through the connected void space, displacing the air in the voids of the originally room-dry sample, to the sample top. The achieved D 2 O saturation of the sample was estimated by volume replacement to be 70-85+%. More precise saturation calculations are not viable. The SEM images confirm a highly variable void size, and although these images could be used to calculate local porosity, they are of one slice through the volume. Buckman et al. (2017) using the same slice calculated highly variable porosity values across the sample. D 2 O injection was followed by injection of distilled H 2 O from a parallel header tank. This approach allows observations of if, and how, this miscible and essentially single-phase fluid replacement process follows preferred flow pathways under conditions where water had already moved through much of the void system, and so there should have been minimal new capillary effects.
The fluid flow was studied through a series of radiographs acquired in the reference sample orientation without sample rotation, supplemented by several acquisitions for tomographic reconstruction (with sample rotation). The suite of radiographs represent a time sequence of 2D images ('projections') of a non-rotating sample where the varying intensities reflect the absorption and scattering caused by the rock sample and the invading fluid types. Starting and stopping a tomograph acquisition requires a few minutes, so there are gaps of a few minutes in the otherwise continuous radiographic acquisitions. Ideally the tomographs capture the water position through the sample volume at an instant in time. However a reconstruction needs (typically) hundreds of radiographic images taken as the sample is rotated about the vertical axis. This takes time during which the water positions may well have moved. To best approximate a static fluid system, fluid injection was stopped during tomograph acquisition to reduce water movement in the sample as much as possible.
There were three tomograph acquisitions during the flow experiment. The first acquisition, Tomo1, was performed when D 2 O had almost completely transited the sample and before the injected fluid was changed to H 2 O. Its 5 min acquisition was chosen as a compromise between the image resolution and the risk of water movement and associated blurring at this stage of the experiment. A 45 min acquisition, Tomo3, which used longer radiographic acquisitions, was performed at the end of the experiment, after the H 2 O injection was completed and the sample had been left without fluid injection for about 30 min. The longer radiographic acquisition time and resulting higher resolution took advantage of the lower likelihood of water movement. An intermediate acquisition was also made in the early stages of H 2 O injection, but it is not used in this work.
The sequence of radial projections of the rotating sample obtained for each tomographic scan is subsequently reconstructed through a Feldkamp (FDK) back projection algorithm to create 3D datasets of neutron attenuation values. The nearly complete time-sequence of radiographs reveals water content changes (and permits determination of the advancing position of the air-water interface) that are interpreted as being the result of complex interactions between imbibition-dominated processes and viscous-dominated flow patterns. These are described in the following section.

Results
The radiograph sequences and the calculated tomographs are, of course, results in their own right, and are used extensively to analyse fluid movement pathways that are subsequently correlated with the pore and void systems seen in the SEM images. Because of the highly complex and spatially variable void and pore network in the deformed L20 sample, there is no simple 'front' visible at the macroscopic scale. Boundaries between fluid types are identified, but they are not the simple boundaries commonly seen in neutron beam flow experiments (Tudisco et al. 2019). Instead there is a complex intermixing of for example D 2 O entering a lamina while adjacent laminae are still air-filled, of D 2 O filled fractures with adjacent air-filled matrix or H 2 O filled fractures with adjacent matrix and fractures still being D 2 O filled.
Equally, the radiographs reveal time-dependent increases in neutron attenuation at many locations over the projection, which is compatible with changing fluid saturations along the neutron flight path to that point on the detector. But, with the clear textural evidence of many pore-system types in the sample, it is not feasible to make assumptions about the character of the void space along each flight path, and there exist an infinite number of possible distributions of fluid saturation to explain the attenuation at each pixel. The variability of void sizes along the path leads to an expectation that there would exist many different pore-scale saturations at all stages between air-filled and fully water saturated, and this is why the inversion is not possible.
The strength of this experiment is its discovery and verification of the co-occurrence of movement through the matrix pore system and the fracture void system at the same time. This Section is focused on characterising this fluid movement for both D 2 O replacing air, a phase change, and H 2 O replacing D 2 O which is an internal marker in the water system. The consequences of these discoveries for the driving mechanism follow in Sects. 5 and 6.

D 2 O Advancing into Air-Filled Pores and Voids
The experiment starts with a room-dry sample at atmospheric pressure and D 2 O is introduced at the base under a near-constant head of approximately 10 kPa. The sample top is at atmospheric pressure. Approximately one pore volume of D 2 O, as calculated from the sample's bulk porosity, is injected before changing to H 2 O injection which continues until the majority of the D 2 O has been swept. Continuous radiography allows live sweep monitoring: the radiography reference orientation aligns some significantly open fractures along strike (see SM2). Water filling these fractures should be quite discernible but water filling of the less well-oriented features will be much less easily seen. This is because apertures are at maximum a hundred microns wide, and even if fully saturated, they represent only a small net amount of water between the inbound beam and the detector. The general water movement is from base to top, but with common fracture and matrix exchanges, many lateral movements and some local movements towards the base. These exchanges include fracture-to-matrix flow, and flow from matrix-to-fracture, especially in the initial saturation stages.
Images from Tomo1, acquired when D 2 O reached the sample top are shown in Fig. 11. The D 2 O arrival time map (Fig. 12) was determined by fitting a sigmoidal function to the evolution of the grey value for each pixel of the radiography sequence. The arrival time  (Fig. 7). Numbers are equal intervals in the radiograph sequence. Each arrow shows the main progress of fluid movement over that interval also moving up along an increasing number of steep fractures ( Fig. 13 steps 2-5) into the upper laminae of the Base Zone ( Fig. 13 steps 6-10). New fractures are occupied at Fig. 13 step 8. Fluid initially moves faster along fractures (Fig. 12 right side Base Zone), and then stops while fluid moves into adjacent matrix. Fluid movement is occasionally down fractures backfilling previously by-passed lower laminae (Fig. 12). Filling is earlier (faster) on the right, reaching broad filling at 10 minutes in contrast to 20 minutes on the left. D 2 O reaches the LM Zone base on the right side after 5 minutes at the same time as the left side base starts to fill (Fig. 12). At all stages, steeply dipping, semi-tabular fractures are filling with D 2 O at the same time as a wide and rather diffuse lamina-parallel zone of filling is observed. Speeds of the nominal frontal advances, along these two pathway types, are similar, although this does not necessarily mean that the fluxes are similar. See also SM3, 4. It is not possible to determine from a radiograph if any image darkening that is not the characteristic tabular fracture shape is caused by D 2 O in a face-on fracture or by D 2 O in the matrix laminae. (SM3, 5). But NTs can resolve this. Distinct water-containing open fractures can be identified in any orientation, and regions of water-filled matrix can equally be distinguished by their laterally extensive semi-tabular shape.
Three NTs and one XRT are available (Figs. 4,11;SM6). (See also SM4). Tomo1 confirms that the Base Zone filling to the point of D 2 O reaching the sample top is correctly attributed to both steep fracture filling and to D 2 O entering and filling some but not all the individual laminae. The lack of any macroscopic (at decimicron resolution) open fractures in the XRT supports this observation as do the macroscopically fracture free regions of the SEM image.

D 2 O Advancing into LM Zone
The LM Zone has a large number of laminae of a different colour and general appearance and are patchily highly porous (Buckman et al. 2017). There are few visible through-going fractures, and it is unclear how much of their extent through this Zone is open and how much is closed. Figure 13 shows that D 2 O enters the laminae of the LM Zone in a distinctly sequential fashion (8-10+), filling one lamina before moving upwards to a higher lamina from a single connecting point with those upward jumps likely to be assisted by, perhaps very short, open fractures. Figure 12 also identified lamina-parallel filling with approximately alternating faster and slower D 2 O entry. Figs. 12, 15 both show a general right to left D 2 O movement, consistent with the dominant point of delivery of D 2 O to the base of the LM Zone. Also see SM3. Figure 12 shows a steep sharp boundary dividing faster right and slower left D 2 O entry. Tomo1 (Fig. 11: also see SM4) shows high water occupancy in the LM Zone laminae. But at the resolution achieved, water content in each voxel is far too high to differentiate what small-scale features the water occupies. So the best guide on what fractures are available is the XRT in the LM Zone (Fig. 14: also Fig. 4, SM6). A number of open fractures are indicated but with limited vertical extent. These open fractures appear to be part of the same set below and above the LM Zone, and are even seen on the photographs of the sample's perimeter (Fig. 3). They are not convincingly continuous open fractures and there are strong indications of en-echelon arrays where individual fractures terminate and a similarly oriented new fracture begins a millimetre or less away. Most of these fractures also lack convincing continuity laterally. However the sharp boundary seen in the LM and Base Zones in Fig. 12 suggests there is one set of connected open fractures from the Base Zone to the top of the LM Zone. Figure 14c, e shows candidate fractures.

D 2 O Advancing into UM and Upper Zone
The radiographs of Fig. 15 show D 2 O entering the air-filled laminae of the UM Zone from the LM Zone ( Fig. 15; see also SM3). They do not distinguish if the fluid exits the LM Zone via fractures or as a continuation of the lamina-jump process noted above. Figure 12 shows a slow upward movement, and possibly a hiatus, before much D2O enters the UM Zone. They also show D 2 O entry into the lowest UM Zone as slow. In either case the UM Zone fractures are certainly progressively D 2 O filled. Figures 12, 15 show movement up several open fractures and into the matrix with filling times increasing outwards from the fractures. The resulting pattern is of a widening region of partly filled laminae in a pattern approximately symmetric about the sample axis (SM3). Note that D 2 O entry into the lowest section of this zone is poor.
While initial fluid ingress into the Upper Zone is similar to the Base Zone, fluid advance into the Upper Zone is not well imaged in this experiment. This is probably caused by the fluid exiting the sample top along fractures, ponding and moving downwards. Nominal filling of the sample's pore space with D 2 O, even while recognising that saturations are not everywhere close to 100%, was the planned point to switch to H 2 O injection.  The dry XRT (Fig. 4) shows the steep open fractures in the UM and Upper Zones. Figure 11 confirms that D2O has entered these fracture and begun to move into the adjacent matrix, particularly towards nearby open fractures.

H 2 O Advancing into D 2 O-Filled Pores and Voids
When the injected liquid changes to H 2 O so that H 2 O replaces D 2 O and unswept air, the interface between these two very similar forms of liquid water remains distinct for the timescale of the experiment, with H 2 O exhibiting a more pronounced neutron flux attenuation. As a consequence fewer neutrons reach the detector, and the detector data shifts to a higher opacity, closer to the dark noise of the camera. This makes it more difficult to extract subtle effects though certain characteristics remain clear. H 2 O moves rapidly through the sample by moving up the fracture systems with very little flow into or out of the undeformed matrix void system. This is best seen in the Base and LM Zones.
Images from Tomo3, acquired at the end of the flow test after H 2 O reached the sample top, are shown in Fig. 16. See also SM7.
Base Zone In the early stages of H 2 O injection H 2 O has a strong tendency to occupy and move through the fracture system not the matrix (Fig. 17), leaving some matrix still occupied by D 2 O or air. This is in stark contrast to replacement of air by D 2 O where both fractures and matrix are important flow path components. Over timescales of an hour plus after H 2 O inflow stops, H 2 O moves from fractures into the matrix, as shown in Tomo3. Thin fractures have halos representing H 2 O that has moved into the matrix. This movement is possibly by pressure-driven flow, but given that no H 2 O is being injected, it is probably mostly by imbibition. Tomo3 shows attenuation values in the Base Zone matrix that are consistent with a mixture of D 2 O and H 2 O. The fractures stand out in the radiographs and the tomograph, and their attenuation is consistent with the fluid type being dominantly H 2 O, indicating a high level of replacement of D 2 O in the fractures by H 2 O in the Base Zone.
LM Zone In the LM Zone the details of the water movement pattern as H 2 O replaces D 2 O in the matrix laminae is rather difficult to determine. Figure 18 indicates very rapid transit of the LM Zone by the H 2 O/D 2 O interface while Fig. 17 shows H 2 O filling over a few minute interval. It also shows variations in filling time start for different laminae with typically very rapid filling along each lamina or laminae group. These observations are compatible with rapid flow up the fractures, and fast flow into suitable laminae that are Fig. 16 Images from Tomo3 cut by the fractures. This suggested process is similar to the behaviour observed for H 2 O movement through the Base Zone, but in the LM Zone the flow occurs much more rapidly. The high H 2 O saturations in the LM Zone are confirmed by Tomo3 (Fig. 16).
UM and Upper Zones The movement of H 2 O in the UM Zone is comparable to the Base Zone, with flow dominated by focused fracture flow and limited matrix. Because the experiment doesn't create any back pressure at the top, once the H 2 O reaches the top and experiences no further resistance, the short-term flow pattern settles into a fracture-dominated mode and no further changes are observed.
A final brief comment is warranted. For operational reasons, the flow experiment was left undisturbed for ten hours. Radiography obtained at this time reveals that the distinct H 2 O/D 2 O separations are no longer visible. The inference is that diffusion caused mixing of the different isotopic forms of water, leading to a characteristic attenuation value that is intermediate between those observed for the initial D 2 O ingress and those seen with the H 2 O replacement process. A short flow experiment, using H 2 O as the injecting fluid, revealed that the short-term flow pattern is almost exclusively using the fractures.

Summary of Important Flow Phenomena
This Section summarises the complex suite of fluid-flow observations detailed in the Results Section.

D 2 O into Air
Base The final distribution of Base Zone D 2 O is substantially the consequence of alonglamina flow through the laminae's pore systems, supplied by water that arrived via a sparse and partly connected network of induced fractures. D 2 O is also contained in these fractures. In some laminae, the D 2 O arrived from superjacent or subjacent laminae. Many 5 μm and smaller fracture-like features, making potential connections across the previously unfilled laminae, are seen in the SEM images. This pattern of water movement along both fractures and laminae in a non-regular fashion is very unlike the picture conjured for a piston-like advance.
LM Zone D 2 O enters the laminae, sequentially filling one lamina before moving to a higher one, apparently from a single point along the inter-lamina interface. Very short, open fractures are the likely explanation of the observed points of connection. D 2 O laminae filling tends to be alternately faster and slower and generally right to left, consistent with the dominant D 2 O delivery from the right of Base Zone. A steep, sharp boundary divides faster (right) and slower (left) D 2 O movement. A number of open fractures are visible in the XRT and Tomo1, but most do not extend fully across the LM Zone. Instead they form en-echelon arrays of closely spaced fractures, but separated by mm scale gaps. One set of open fractures connects the Base Zone to the LM Zone top, in the same position as the steep boundary in Fig. 12, suggesting that the pathway up an open fracture and into the matrix was more favourable than was a jump across an open fracture in that same lamina.
UM and Upper Zones D 2 O enters the air-filled laminae from the LM Zone, but it is unclear if the fluid arrives via longer fractures or using the lamina-jump process. Upward movement slows or pauses as the lowest UM Zone filling. D 2 O clearly moves up several open fractures and outward into the matrix with filling times increasing outwards from the fractures. The resulting pattern is of a widening region of partly filled laminae in a pattern symmetric about the sample axis.

H 2 O into D 2 O
Base Zone H 2 O mainly moves through the fractures not the matrix, leaving much of the matrix occupied by D 2 O or air.
LM Zone There is a very rapid H 2 O transit of the LM Zone. There are also variations in the start of filling time for different laminae, with typically very rapid filling along each lamina or laminae group.
UM and Upper Zones The movement of H 2 O in the UM Zone is comparable to the Base Zone, with flow dominated by focused fracture flow and limited matrix invasion.

De-Saturation
A reduction of laminae water saturation occurs in some places as D 2 O invasion continues. During Base Zone filling, the radiographs progressively darken as water enters laminae. But, then some regions start to lighten, indicating a reduction of net water content. We infer this to be de-saturation. In Sect. 4, we documented water movement along the larger fractures, visible ahead of the current lamina filling, occurs in rapid steps. Those rapid water advances into the fractures coincide with episodes when the laminae experience a reduction in grey-scale values. The correlation between the fracture-water flow and the decrease of water content in some laminae leads to the hypothesis that the capillary 'pull' within the fracture, or elsewhere, is capable of removing water contained in the laminae connected to that fracture. The implications for operative physics is examined in the Discussion.

Discussion
The small pore and void sizes and the context of water replacing air lead to the expectation that capillary effects are significant. The first sub-section explores the role of capillary physics to derive process explanations for the responses observed including the desaturation noted above. The second sub-section considers the approximately saturated fluid flow (H 2 O replacing D 2 O) through the fracture network, including why there is so little evidence of new water movement into the matrix, even into laminae which contain residual air. The final sub-section provides a perspective on possible next steps involved in testing the nascent hypotheses that emerge from the previous two sub-sections.

Interpretation of Partial-Saturation Flow Processes
Consider first the lamina 'jumps', where water initially invades only one lamina and then moves into a nearby lamina, transiting the interface(s) at a very localised entry point. This observed behaviour implies common texturally related impediments preventing wholesale water crossing of laminae interfaces. The SEM images show no relative lack of pores at lamina tops or bases but neither do they show fully open pathways. This is as is expected in a micrite. But the SEM images are of a slice and it is not possible to exclude connectivities, or lack of connectivities, that would be visible in a differently oriented thin section. Also the pore network could be significantly anisotropic, with fewer and poorer connections in the vertical direction.
We cannot completely exclude the possibility that the pore systems adjacent to the lamina boundaries possess both reduced vertical connections and extremely small pore throats. But such small throats are expected to be accompanied by extreme capillary forces. If that were the case, then wouldn't the water respond by moving across these lamina interfaces as a preference? Perhaps the difficulty here is the expectation that high capillary forces always cause an increase in fluid flux. This naïve supposition ignores the issue of rates. Even with a large driving force, fluid moving through a pore-system restriction can only do so when the fluid entering the restriction is replaced upstream by other fluid. In this experiment the water has not achieved full pore saturation when these jumps occur. Thus, the drive that would cause flow across the interface is not sufficient to cause a discernible flux. Recall that such a flux is in competition with water being drawn into the partially saturated lamina. This issue is one where the analysis cannot be based on state parameters only. Instead, it needs to be formulated in terms of the system energy changes in the wetting-phase fluid (the water), and other changes such as compression of trapped air, in a domain large enough to account for multiple-pore fluid movements.
The fluid jump observations indicate that single points are involved. Moreover, once a jump event begins, the adjacent lamina exhibits lateral filling with water, progressing away from that point. As noted above, fractures of a few-μm size are observed in the SEM images in locations suitable to a notional flow path across the lamina boundary. The differences in image resolution, and in registering data from multiple sources (neutron radiographs and SEM) are acknowledged. This limits our ability to associate flow observations with specific fractures that we expect to extend away from the thin-section plane and into the sample, where their role in jump events could be captured by the radiographs. Although our present inability to directly associated individual small fractures with the jumps is disappointing, the cause-and-effect relationship seems likely.
The fluid motion after the jump introduces a further perplexing issue. The gist of the argument in the preceding paragraphs introduced a speculation that the already-filled lamina could not deliver an increase in fluid flux, once a jump opportunity occurs in the pore system of the overlying lamina. But, the flow response when a fracture is breached indicates that flux is not overly limited in the subjacent lamina. Thus, any limits on fluid flux along the initial lamina (and thus the energy terms in the system analysis) are associated with conditions that are not strictly local.
The obvious candidate for the putative non-local control is the significant fracture that initially supplies the fluid into the lamina. How could this larger fracture fail to be able to supply increased flux into the partially filled lamina? As noted in Results, and further considered in the next sub-section, the fully saturated part of the experiment leads to significant flux that is concentrated in the main fractures. So, why would the fractures impose any limits on flow in the unsaturated case?
One possibility, which is speculative, is that the main fractures are capillary limited during the early part of the experiment, when the lamina jumps occur. We noted in the Results that we observe comparable, but larger, jumps during the initial filling of the main fractures, and this is compatible with the speculation. But, the fractures are connected to the sample bottom, and the adjacent water chamber. Surely, the fractures must be experiencing Darcy (pressure-driven) flow? Perhaps they are not. Or it may be only the lower portions of the fracture system that become sufficiently saturated to operate according that set of flow physics. The hydraulic head imposed at the base of the sample is quite low (10 kPa) and the time-average flux rate is very small. It is conceivable, but cannot be directly tested with our data, that fluid entry into the sample is so restricted that little or none of the sample reaches high enough saturation during the initial experimental step to convey the inlet fluid pressure into significant parts of the sample interior. If so, then the later stage of the experiment, where saturated conditions lead to localised flow in the fracture network, is not a good indicator of what is occurring during the preceding stage.
The scope of this paper does not encompass flow simulations that could provide insights into the matters discussed in this sub-section. Regardless, the observations, and then the inferences about the operative physics, point towards some very complex flow responses while the sample is only partially saturated. In developing potential explanations, the range of behaviours allows us to test a provisional hypothesis against flow observations from a different location in the sample. The outcome of the analysis is that the interplay between distinct textural types can lead to capillary-dominant flow systems that exhibit significant levels of complexity in both spatial and temporal terms. Understanding of this complexity seems to demand a full system analysis and simulation that avoids imposing state conditions at convenient points in the sample system.

Saturated Flow Patterns in Some Fractures
When H 2 O replaces D 2 O in near-saturated conditions almost the entire flow occurs in the fracture system, at least during the flow experiment. Many more fractures played a role in while achieving water saturation using D 2 O compared to the number of fractures involved in the saturated situation as revealed by H 2 O replacing D 2 O. Tomo3 shows that even within prominent fractures, the flow was concentrated along only portions of the available fractures.
It seems reasonable to suppose that the H 2 O movement took place in the most energetically favourable places within the system. The resulting flow paths do not appear to be strongly linked to fracture intersections or to other geometric characteristics that could be used to make pre-experiment flow path predictions. Some pathways were used in the saturation (by D 2 O) stage, while others are less obvious reusing those fracture segments. These facts point to a strong influence (in the two conditions) of self-organisation ideas, in which the fluids make use of those available paths that allow the most-favourable energy dissipation.

Perspectives on Next Steps
This experiment, which involves a micro-porous rock modified by imposed small inelastic strains reveals considerable complexity in both multi-phase and single-phase fluid systems. Additional experiments, involving similar as well as different rock types, have been completed, and those results (using methods like those described here) are to be the subject of other publications. However, the analysis of this sample demand a new suite of experiments. Some of the hypothesised pore-scale physics interactions might be suitable for study via manufactured 2.5D micro-models (Sohrabi et al. 2004), possibly allowing the demonstration of long-range interactions of the sort identified here.
A very obvious next step is to use digital-rock methods to simulate of pore systems' fluid-flow (Blunt et al. 2012;Andrä et al. 2013;Berg et al. 2013;Bultreys et al. 2016), and from which it is notionally possible to test ideas such as those posed here. The porenetwork method (McDougall and Sorbie 1997;Ryazanov et al. 2009) has been extended to allow the construction of models containing spatially limited pore systems of contrasting characteristic scales (Jiang et al. 2018), now with the ability to include fractures of arbitrary textural complexity. Although these methods are relatively well-advanced, they have not given strong consideration to slit-like micro-pores. This paper provides a motivation to undertake additional development of the digital-rock toolset.

Conclusions
Water flow experiments in the neutron beam introduce D 2 O into an air-filled, experimentally deformed, extremely fine-grained limestone, followed by H 2 O injected into the water-saturated sample. The three pore-fluids (air, D 2 O and H 2 O) each have distinct attenuation effects and allow tracking of the flow processes. Radiographs and tomographs reveal a complex spatial and temporal interplay between imbibition and pressure-driven flow at a spatial resolution of 40 μm. Analysis of the continuum-scale images leads to process interpretations that imply long-range interactions of capillarydominated movements, including local de-saturation associated with advances of water elsewhere. In the fully saturated portion of the experiment, fluid flow is focused sections of the connected fracture network, with little exchange into or out of the matrix, and with stagnant water regions in many fractures. SEM images add textural evidence that explain the local property variations implied by the process model. They document both micro-fractures that explain the fluid communication across laminae interfaces and macro-scale fractures that are transitional between open fractures and deformation bands. This work has implications concerning the complexity of multi-phase fluid flow in fine-grained rocks that have been subjected to only very slight magnitudes of inelastic strain, imposing a range of discontinuity types onto the previous rock. The results reported here have relevance in a wide range of extractive, storage and disposal activities involving single-and multi-phase fluid flows through micro-porous rock types.
Author Contributions HL conceived the research topic and acquired free-at-point-of use beamline access at the ILL facility, Grenoble France. HL and GC designed and created the experimental equipment with ET who contributed parts of her equipment and her protocols in the final experimental design. AT designed the ILL facility neutron and X-ray experimental protocol and performed the data collection and the initial data processing and analysis. GV advised on the experimental design and on the analysis of the neutron and X-ray results. AT, ET and ME ran the experiments with assistance from GC and HL. JB performed all the SEM acquisitions and processed all the data. HL and ET completed the final neutron and X-ray data analysis with advice from AT. SAH performed the XRT study that preceded the neutron beam flow experiments. The first draft of the manuscript was written by HL and all authors commented on previous versions of the manuscript. All authors read and approved the final manuscript.

Funding
The main the experiments were run using the instrument NeXT at the Institut Max von Laue-Paul Langevin (ILL). This work has been carried out in the framework of NI-Matters, a joint research unit of Universite Grenoble Alpes (UGA), the ILL, and Helmholtz Zentrum Berlin (HZB). No direct funding was provided but beamtime was free at the point of use.

Data Availability
The datasets generated during the current study are held by the ILL (Institut Max von Laue-Paul Langevin) and will be made available on request in due course according to their protocols.

Conflict of interest
The authors have no relevant financial or non-financial interests to disclose to the present.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.