In this work, porous carbonate rocks from Maiella Mountain (Abruzzo, Italy) and Favignana Island (Sicily, Italy) were investigated. The Favignana grainstones (Early Pleistocene in age) consist of well-preserved bioclasts composed of Vermetus, Serpula, bivalves, echinoids, red algae, and corals ranging in size from submillimeter to centimeter (Tondi et al., 2012). The host rock is poorly cemented with the cement limited to the grain contacts, around echinoids, or within intragranular pores (Tondi et al., 2012). Its pore-network is mainly composed of well-connected intergranular pores characterized by a dominant diameter ranging between 100 and 400 μm (Tondi, 2007). In addition, isolated porosity is represented by intragranular space within the bioclast (Zambrano et al., 2017). The Favignana grainstones experienced a shallow burial depth of approximated 30 m (Tondi et al., 2012; Antonellini et al., 2014). The carbonate grainstones from the Maiella Mountain (Campanian to Maastrichtian in age), here after called Orfento grainstones, mainly consist of rudist fragments interpreted as a package of proximal bioclastic turbidites (Mutti, 1995). Detailed X-ray micro-CT analyses indicated that the pore space is characterized by both micro- and macro- pores (Tondi et al., 2006; Ji et al., 2015; Zambrano et al., 2017), however permeability is mainly controlled by pores with diameter larger than ca. 10 μm (Zambrano et al., 2018). These rocks experienced a maximum burial depth between 0.5 and 3 km (Ori et al., 1986; Graham et al., 2003; Rustichelli et al., 2016).

In general, both rocks have similar total porosity values (near 30%), however their permeability differs by two orders of magnitude (for details see Table 2). Zambrano et al. (2017, 2018) have explained these differences in permeability as related to diverse values of connected porosity, connectivity, pore size distribution, and tortuosity. In fact, the connected porosity observed by means of X-ray micro-CT (pixel size 9 μm) for the Favignana grainstones is near 27%, whereas the Orfento grainstones is less than 15% (Zambrano et al., 2017). Also, connectivity density (a value representing the number of redundant connections between pores) is higher for the Favignana grainstones. On the other hand, the specific surface area is greater (indicating a smaller grain size) for the Orfento grainstone (near 79.6 mm^–2) in comparison to the Favignana one (18.1–38.1 mm^–2). Moreover, the tortuosity is almost double for the Orfento grainstones in front of the Favignana one (Zambrano et al., 2018).

For the experiments, cylindrical samples were prepared with a diameter of 10 and 20 mm in height (Figure 1). These dimensions were chosen so that the cross-sectional water thickness would not be greater than 5 mm in order to get suitable neutron transmission. Considering the previous reported porosity (Table 2), it was estimated a maximum water thickness of 4.5 mm. Samples were placed in cylindrical holders made from Polytetrafluoroethylene (PTFE), commercially known as Teflon, with 12 mm of inner diameter closing the remaining space between sample and holder with Teflon tape. The choice of Teflon, for both sample holder and tape, was due to its suitability for both neutron and X-ray imaging (Domanus, 1992¹). The water container was attached to a rotation stage by means of an adapter that allows the connection of a pipe for water without interfering with the experiment.

The experiments were performed at the beamline for Imaging with COld Neutrons (ICON) at Swiss spallation neutron source (SINQ) of the Paul Scherrer Institut (PSI²) in Switzerland. Details of the ICON facilities and available instrumentations have been described by Kaestner et al. (2011). Bimodal imaging was performed using neutrons and X-rays, sequentially. The schematic layout is given in Figure 2a (Kaestner et al., 2017) and a photo of the setup is given in Figure 2b (Kaestner et al., 2016).

Time-lapsed or dynamic NR was performed for two different samples pertaining to Favignana and Orfento grainstones. Similar to Cnudde et al. (2008), we have performed a dynamic NR experiments to investigate the position of the wetting front as a function of time. Our experiments differ from Cnudde et al. (2008) in that we did not use the gravity as force for triggering the flow, instead we used only constant flux rate (forced imbibition). In order to introduce deionized water into the system, a dual syringe pump (Pharmacia P-500) was used with a flow rate set to 1 ml/h. The pumping device can provide a maximum pressure of 5.0 MPa. Since the studied samples are highly porous and pervasive, it was expected (and tested) that the threshold pressure was not be reached during the experiments.

The output data consist of a time series of 2D NR images. With the aim to map the dynamic flow front in each sample, at the selected flow rate, about 26 neutron radiographs (one radiograph every 90 s) were acquired until the sample was filled with water. For corrections, 5 images of dark current, 5 flat field beam images, and a radiograph of the dry sample were taken before the experiment started. To indicate the presence of water, the original dry condition image is subtracted to radiograph to the n-time radiograph. During the processing, the noise was removed with a median filter by using the freeware ImageJ (Schneider et al., 2012).

For estimating a pseudo water saturation trough time, a series of profiles were made in each photograph. The first intention of this analysis is to roughly estimate the advancement of water flow and to distinguish important heterogeneities within the samples. The pseudo water saturation was estimated by processing the raw data X_awr voxel by voxel following the equation (partially modify) proposed by Riaz et al. (2007):

Where the subscripts ar and wr refer to air-filled and water-filled rock, respectively. The first case corresponds to images taken before the flow started (dry condition), whereas the water-filled stage was arbitrary selected as the time when the water started to spill out of the sample. The rate of water volume injected into the sample is given by the flux rate of the pump. In comparison, it was estimated the fraction of water within the invaded volume of rock a certain time. This data is useful to distinguish when the imbibition is dominated by capillarity forces or the constant flux rate.

For the static tomographic experiments, a single sample corresponding to Favignana grainstones was selected. This approach consists of three steps: (1) XCT of a sample in dry condition, (2) NCT of the same sample in dry, and (3) wet conditions (water-filled).

An XCT scan was recorded with a pixel size of 13.5 μm. A set of 375 projections were taken over 360 degrees with air on both sides of the sample. The exposure time per projection was 30 s. For the NCT at dry and wet conditions, the pixel size was also set to 13.5 μm. A total of 625 projections were planned over 360°; however, for the dry case the acquisition was stopped after about 280 degrees due to time limitations. The exposure time was 80 s per projection. Before every CT scan, 5 dark current and 5 open beam images were acquired.

The reconstruction was done by the Octopus 8.8 software (Vlassenbroeck et al., 2006). The beam was chosen as having a cone-beam geometry although it was nearly parallel, but this option was better as the beam is slightly divergent. Then, the noise was removed using a 3D Median filter in ImageJ. Ring artifacts on the reconstructed slices were removed using the modified Sijbers & Postonov algorithm (Brun et al., 2011) which is part of the Pore3D software library (Brun et al., 2010; Zandomeneghi et al., 2010). For the NCT with incomplete projections, the MuhRec tomographic reconstruction software was used (Kaestner, 2011). The filtering of the data was done using the ISS filter in KipTool (the in-house developed image processing tool at ICON).

For the segmentation of the X-ray images in dual components, voids, and solid phase (framework grains and cement), the automatic multiphase k-means clustering algorithm (Hartigan, 1975; Hartigan and Wong, 1979) was used. The algorithm was set to 4 classes and implemented twice in order to get most of the visible pores. Results are binary images composed of voids and grains. After this, the tool ‘Find Connected Structures’ of the open source Fiji software (Schindelin et al., 2012) was used for dividing the pore space into two components: (i) connected pores and (ii) isolated pores. For the simulations, only the connected pore networks were used for facilitating the computation.

Considering the attenuation factors for the different studied phases (water, air, and carbonate rock; Table 1), a rough indicator of the water content was founded by using the difference of NCT wet-dry images. In the gray color intensity scale, wet zones tend to have a more whitish color in comparison to air or carbonate rock. Therefore, warmer colors indicate a higher saturation of water.

Since the laboratory experiments were performed on dry samples, Lattice-Boltzmann simulations were useful to investigate the flow at saturated condition. The simulations were executed by means of the open-source computational fluid dynamics software PALABOS (Latt, 2009) using the methodology introduced by Degruyter et al. (2010) and later modified by Zambrano et al. (2018). The methodology consists of imposing a single-phase fluid flow through the segmented 3D images by maintaining a fixed pressure gradient between the inlet and outlet faces of the volume, while the rest of the faces were padded. The interface pore-voids was converted to bounce-back boundary conditions. The main difference with the methodology proposed by Degruyter et al. (2010) is the replacement of the collisional operator BGK by an MRT (d’Humières et al., 2002) with a D3Q19 lattice, which guarantees results viscosity independent (Zambrano et al., 2018). The steady state condition of the simulation is considered after the standard deviation of the average energy falls blow 10^–4 after 1000 time-steps, Degruyter et al., 2010). After that, the permeability component parallel to the imposed flow can be computed applying the Darcy’s law,

where, δP/δx is the pressure gradient, μ the fluid kinematic viscosity, and U the average fluid velocity per unit of area. All the variables are handled in lattice units before the permeability computation, results are transformed to SI units multiplying by the effective length of the voxel side in meters.

The results of the experiment consist of a series of time lapse neutron radiographs obtained at a time interval of 90 s indicating the flow front trough time (Figures 3, 4). In addition, longitudinal profiles, flow front progression and water saturation graphs were plotted to show the results in a more quantitative manner (Figures 5, 6).

In the case of the Favignana grainstones, the studied sample seems homogeneous at initial dry conditions (NR 90s). However, the time-lapsed neutron radiographs (Figure 3) show that the flow front is slightly inclined (NR 90s -NR 990s). Initially the water invaded the rock sample rapidly, however the invaded volume is not highly saturated (NR 90s -NR 990s). The saturation increment is more evident in the Figures 5A, 6A. After the sample is totally invaded by water, the water saturation increases (NR 1290s – NR 1800s) until the water spilled out (NR 2250s). The profiles confirm both the initial rapid migration of the flow front and continuous saturation of the sample trough time until the water spilled out (Figures 5A, 6B). The profiles also indicate the heterogeneities in the sample associated with local porosity variation. In general, these heterogeneities are kept in time but their general intensity changes as an indicator of an increment of the water saturation.

In the case of the Orfento grainstone, the initial image in dry conditions (NR 90 s) shows a sample with smaller grain/pore size with respect to Favignana grainstone. Whereas, the time lapsed neutron radiographs, once the flow started (Figure 4), rapidly show the presence of lineations (fractures) characterized by an earlier saturation. These fractures could be created during the preparation of the sample or likely of natural origin. These features control an anisotropic distribution of the flow. The fractures are also recognized as important heterogeneities controlling the fluid flow in the longitudinal profiles (Figure 5B). In this graph, an ‘anomalous’ rapid increment in water saturation is observed at the time 1260s in correspondence to the presence of fractures. A similar migration of the flow front and increment of the relative saturation is showed by the profiles through time (Figure 6B).

We have observed that water saturation was driven by to mechanism: capillarity and constant flux rate. The limit between the two flows mechanism correspond to the time when the water reaches the top of the rock samples. After that, the saturation starts to increase uniformly according to the constant flux rate provided by the pump. This limit seems to be preceded by an inflection of the fraction of water curve. During the capillarity driven flow, the waterfront velocity for the Favignana samples is slightly lower (0.0202 mm/s) than for the Orfento sample (0.0208 mm/s). In the Figure 6, the flow front vs. square root of time presented some scattering and variability with the expected best fit line. This variability is more important for the Orfento grainstones, which is characterized by the presence of fractures.

The maximum fraction of water absorbed by the samples seems a good indicator of the connected porosity. In the case of Favignana sample, the fraction of water reached during the capillarity driven flow is near 12.53%, whereas the maximum, after the saturation of the sample, is near 32.4%. For the Orfento sample the fraction of water reached during the capillarity driven flow is near 11.5%, and the maximum at the end of the experiment is near 20.45%.

Imaged from both methods were aligned and presented with the same nominal resolution (pixel size equal to 13.5 μm) for comparing and integrating the results (Figure 7). The XCT data provides a good evaluation of rock components (framework grains, matrix, cement, and pores). Particularly the resolution of the images allows to distinguish and segment the macropores (diameters between 50 and 800 μm) contained in the sample (Figure 7A). It is expected that some pores below the nominal resolution (pixel size equal to 13.5 μm) are present in the sample.

The strength of neutron imaging is its sensitivity to hydrogen and hence to the fluid in the sample. As shown in the Table 1 neutrons give a very high contrast between the three main phases in the studied samples (i.e., water, air, and carbonate rock). The limitation of the neutron image is that the pores/grains are not very well-defined (images are blurry in comparison to XCT). However, the method proved useful for comparing the wetting-phase saturation after the inhibition experiment. In fact, NCT images at dry and wet conditions look very different (Figures 7B,C, respectively). In the use gray color bar, wet zones tend to have a more whitish color in comparison to air or carbonate rock. This is likely due to the sensitivity of NR to hydrogen and the different attenuation factors for the different studied phases (water, air, and carbonate rock). A rough estimation of the relative concentration of water considering the difference of NCT wet-dry images (Figure 7E). In this figure, warmer colors indicate a higher saturation of water.

Some differences between the pores extracted from X-ray images and the pores filled by water (differences of NCT wet-dry images) were observed (Figures 7, 8). The first case some detected pores are not filled by water, this could be related to the presence of isolated intragrain (e.g., bioclast) pores previously reported by Zambrano et al., 2017). Also, this could be attributed to the presence of air bubbles trapped due to a ‘rapid’ and low-pressure water flow. This fact may explain the presence of macropores partially filled with water. In addition, the water could have taken preferable pathways like the borders of the sample even though it was well-sealed by Teflon tape. Another important difference, correspond to apparent zones filled by water that do not correspond to the observed pores. This is more likely due to the presence of micropores (sub-resolvable pore by XCT) that are contributing to the water flow. The structure of these pores may be not clearly defined in the images, however the water presence within the pores can be mapped.

The fluid flow simulations provide information about the possible behavior of the flow through time after saturated conditions. In steady state conditions we compared the velocity field (Figure 7F) with the connected pore network (Figure 7D) extracted from the X-ray micro-CT images. Results indicate that some pores do not contribute to the flow in simulation; however, in reality these pores are connected (Figure 8, label 2). This is likely an issue of image resolution of the XCT. Similarly, some pores below the resolution (XCT) contributing to the flow (NCT) are not detected and therefore not included in the simulation (Figure 8, label 3). Also, pores located at the lateral edges could be excluded if they are not connected to the top and base of the sample. Therefore, the pores contributing to flow in the simulation are only macropores clearly connected (Figure 9). In these regards, the computed permeability (1.29×10⁻¹³ m²) is slightly lower than the values obtained by Zambrano et al. (2018) using X-ray micro-CT at higher resolution (see Table 2). In fact, the calculated permeability tends to be underestimated at low spatial resolution of pore network model (see Zambrano et al., 2018).

The dynamic 2D NR permits the evaluation of the front of the fluid flow as a function of time. However, we detected an initial front likely driven by capillarity forces and later increment of the saturation on the wetting phase. The latter is likely driven by a forced imbibition caused by the constant flux rate. The limit of these domains is given by the time when the waterfront reaches the top of the sample. Also, it was possible to observe the influence of sample heterogeneities on the fluid flow. These heterogeneities can be the intergranular pores or fractures. Particularly, some fractures observed in the Orfento sample could be natural or related to the sample preparation. In contrast to some sandstone samples investigated by Cnudde et al. (2008) we did not observe any ‘abnormal’ dual-waterfront due to these heterogeneities. Nevertheless, the NR methodology provides information in two-dimensions resulting in a smooth boundary water/air. As a consequence, the detail control on the surface geometry and roughness (e.g., Alava et al., 2004) was not possible. This issue could be solved by implementing time-step NCT in further investigations. We found that final values of water absorbed by the samples are near and slightly higher that the reported values of connected porosity using XCT (Zambrano et al., 2017). Therefore, the reached water saturation in the experiments should be close to the available porosity. However, some air bubbles should remain trapped in the rock.

The dynamic 2D NR experiments showed differences regarding pore structure between the two evaluated rocks that have been reported by previous authors (Tondi et al., 2016; Zambrano et al., 2017, 2018). In the case of Favignana grainstones a high heterogeneity of the gray values in the image profiles was observed, whereas the Orfento grainstones were more regular and homogeneous. This is related to the huge differences in texture, grains and pores sizes being near 10 times larger in the Favignana grainstones in comparison to the Orfento ones. This difference in texture (specific surface area, connectivity, tortuosity) is responsible that Favignana grainstones have a permeability two orders of magnitude higher than the Orfento grainstones. Since both rocks have similar porosity, but different pore-throat sizes. The capillarity forces driving the flow at this first stage are expected to be very different. In fact, lithofacies with larger pore-throat sizes are expected to have small initial displacement pressures, and therefore the air phase can be easily displaced by the injected water. Major initial displacement pressures are expected in lithofacies with smaller pore-throat sizes, which may buffer the displacement of air by water phase. Therefore, differences in the flow velocity in the capillarity driven flow domain, may be related to both the difference in permeability and pore geometry (diameter, connectivity, tortuosity) between the two rocks as have been reported by Zambrano et al. (2017, 2018).

The combined use of XCT and NCT methodologies can provide information regarding the pore space architecture (X-rays) and the fluid contents (neutrons). According to our results, X-rays can define the geometry of the solid components of the rocks (framework grains, cements, and matrix). Neutrons, however, show images more diffused and only a few macropores are distinguished. On the other hand, comparing NCT in dry and wet conditions allows to discriminate the presence of water within pores. In fact, the method seems to provide information on micropores contributing to fluid flow even below the actual spatial resolution of the images. On the other hand, the NCT at wet conditions reveled air trapped in the pores. This issue may be related to the initial unsaturated and low-pressure conditions of the experiment.

This approach using NR and integrated XCT an NCT is relevant for the study of petroleum systems, where the capillarity governs the initial distribution of the different phases of fluids within the reservoir (Brown, 1951), control the relative permeabilities, pressure and saturation of different phases during production (Jerauld and Rathmell, 1997; Christiansen, 2001). Particularly, the implemented approaches were useful for mapping fluid flow as a function of time, wetting phase (water) saturation dynamically and at the end of the experiments. In these regards, the methodology finds uses to the study of spontaneous and forced imbibition gas recovery (e.g., Di Stefano et al., 2017), irreducible water estimation useful for petrophysical evaluations and oil recovery planning (e.g., Treiber and Owens, 1972; Salathiel, 1973; Morrow and Melrose, 1991). Another important application of this method could be the evaluation of moisture inside building stones (e.g., Cnudde et al., 2008). Here the absorption of water by capillarity can cause important degradation of in building materials (e.g., carbonate rocks).