Verifying Pore Network Models of Imbibition in Rocks Using Time‐Resolved Synchrotron Imaging

At the pore scale, slow invasion of a wetting fluid in porous materials is often modeled with quasi‐static approximations which only consider capillary forces in the form of simple pore‐filling rules. The appropriateness of this approximation, often applied in pore network models, is contested in the literature, reflecting the difficulty of predicting imbibition relative permeability with these models. However, validation by sole comparison to continuum‐scale experiments is prone to induce model overfitting. It has therefore remained unclear whether difficulties generalizing the model performance are caused by errors in the predicted filling sequence or by subsequent calculations. Here, we address this by examining whether such a model can predict the pore‐scale fluid distributions underlying the behavior at the continuum scale. To this end, we compare the fluid arrangement evolution measured in fast synchrotron micro‐CT experiments on two rock types to quasi‐static simulations which implement capillary‐dominated pore filling and snap‐off, including a sophisticated model for cooperative pore filling. The results indicate that such pore network models can, in principle, predict fluid distributions accurately enough to estimate upscaled flow properties of strongly wetted rocks at low capillary numbers.


Introduction
Imbibition is the process where a wetting fluid replaces a nonwetting fluid in a porous medium. It occurs in nature in the topmost layer of soil during rainfall, during the flow of immiscible pollutants in ground water resources, the geological storage of carbon dioxide or hydrogen, and the management of petroleum reservoirs (Atteia et al., 2013;Bickle, 2009;Blunt, 2017). To understand the link between pore structure, wettability, and upscaled multiphase flow parameters, a large amount of effort has been put into simulating the pore-scale fluid distribution and calculating the corresponding macroscopic properties Meakin & Tartakovsky, 2009).
The key challenge in this field is to describe the intricate arrangement of the two fluids in the pore space during the fluid invasion Herring et al., 2013;Zhao et al., 2016). Models struggle to capture the effect of the many fluid-fluid interface movements on the fluid flow over sufficiently large length and time scales (Meakin & Tartakovsky, 2009). While Navier-Stokes solvers are becoming mature enough to perform useful numerical experiments on domains tens of pores across over relatively short time scales (Alpak et al., 2018;Ferrari & Lunati, 2013;McClure et al., 2018;Shams et al., 2018), it remains desirable to define a more conceptual model which describes the emergence of structures at the scale of fluid clusters for upscaling (Hilfer et al., 2015). This is important to upscale multiphase flow properties, especially in heterogeneous porous media where representative elementary volume considerations may be crucial.
Pore network models (PNMs) conceptualize the pore space as a network of larger void volumes, pores, bounded by restrictions, called throats ( Figure 1). In Earth science cases, the global flow rates within reservoirs are typically so low that the capillary force is often assumed to be dominant at the pore scale (Blunt, 2017). The most computationally efficient type of pore network models therefore adopts the quasi-static assumption: Only capillary forces are taken into account and the fluids are assumed in capillary quasi-equilibrium at each moment (Blunt et al., 1992;Wilkinson & Willemsen, 1983). The fluid invasion process is abstracted into filling rules which aim to produce realistic emerging fluid distribution patterns (Blunt et al., 2013). The pore space is invaded as a sequence of invasion events of three types: snap-off, piston-like displacement, and cooperative pore filling (Lenormand et al., 1983). Predicting the evolution of the pore-scale fluid distribution therefore reduces to calculating the entry capillary pressure for each of these events in all available pores and gradually filling the pore space with the invading fluid in order of this pressure. Dynamic pore network models, on the other hand, explicitly capture the viscocapillary force balance of immiscible transport (Joekar-Niasar & Hassanizadeh, 2012;Lenormand et al., 1983), but come at the price of increased complexity and computational demands. While a number of early studies have showed that quasi-static PNMs can match experimental relative permeability curves in specific cases (Blunt et al., 2013;Øren et al., 1998;Valvatne & Blunt, 2004), substantial efforts over several decades have not led to satisfying predictive capabilities in a general sense (Bondino et al., 2012). It has been unclear whether the shortcomings originate from uncertainties in the extracted network, the pore-filling rules, or whether it is adequate to use a quasi-static description at all. Newly available fast time-resolved X-ray micro-CT experiments (Andrew et al., 2015;Berg et al., 2013;Schlüter et al., 2017;Singh et al., 2017;Youssef et al., 2013) have challenged the quasi-static assumption, as viscous and inertial effects such as ganglion dynamics and intermittency have been observed Gao et al., 2017;Reynolds et al., 2017;Rücker et al., 2015). However, it is currently not well understood how strongly this affects the distribution of fluids in the pore space. The validation question is complicated by the fact that model simplifications lead to a large number of internal microscopic degrees of freedom, but the output consists of relatively few macroscopic parameters that can be used for validation. Therefore, there is a risk of overfitting the model when tuning it to experimentally measured continuum scale flow properties, for example, relative permeability and capillary pressure-saturation functions (Aghaei & Piri, 2015;Bondino et al., 2012;Jerauld et al., 2017;Masalmeh et al., 2015). To address this, the problem of validating multiphase flow simulations can be split into two questions . The first is whether the PNM provides a good prediction of the locations and connectivity of each fluid in the pore space. Once this has been established, the second question is how to translate this into relative permeabilities, that is, calculating saturations and pressure or flow fields of each fluid in the PNM. This general strategy is motivated by the fact that the main contribution of the flux in capillary-dominated flows comes from connected flow pathways of each fluid .
While there have been a number of dedicated studies on how to best extract a PNM in order to estimate flow rates and saturations (e.g., Joekar-Niasar et al., 2010;Raeini et al., 2017;Ryazanov et al., 2009), validation of the fluid distributions predicted by PNM has remained relatively underexplored and mainly limited to 2-D Figure 2. The workflow for validation of pore-scale multiphase flow models proposed in this paper. After extracting a pore network model from a micro-CT scan of the sample, we generated a series of fluid distributions in the pore network model by either numerical simulation (independent of flow experiments) or experimental micro-CT imaging. We then investigate whether the simulation succeeds at filling pores in approximately the right order compared to the experiment, particularly with regard to predicting upscaled flow properties. micromodels (e.g., Joekar Niasar et al., 2009;Zhao et al., 2019). In a previous study, we showed how PNM predictions can be mapped to fluid distributions measured with micro-CT to facilitate direct pore-by-pore comparisons in rock samples (Bultreys et al., 2018). Here, we extend this to time-resolved micro-CT data sets to validate the filling sequence rather than a static snapshot, resulting in a more complete picture. The experimental data in this study are two publicly available fast synchrotron micro-CT data sets of waterflooding in strongly water-wet Ketton limestone  and Gildehauser sandstone (Rücker et al., 2015). This study investigates whether a quasi-static assumption can result in realistic filling sequences for strong imbibition at low capillary numbers. Contrary to previous studies (e.g., Berg et al., 2016), the models used here implement snap-off and cooperative pore-filling processes on experimentally based geometries (Ruspini et al., 2017;Valvatne & Blunt, 2004). We used a realistic cooperative pore-filling algorithm (Ruspini et al., 2017) and minimized the influence of the specific PNM extraction and of the boundary conditions on the analysis. The general workflow is depicted in Figure 2. We find that the quasi-static model generally predicts fluid distributions that reproduce mm-scale averaged properties, at least in strongly water-wet samples at low capillary numbers.

Experiments
Details of the Ketton and Gildehauser unsteady-state (constant injection flow rate) imbibition experiments are described in Singh et al. (2017) and Rücker et al. (2015). Both data sets are publicly available Singh et al., 2018). Ketton and Gildehauser (the latter a local variety of Bentheimer), respectively, have a porosity of 23% and 20%; and a permeability of 2.8 × 10 −12 m 2 and 1.48 × 10 −12 m 2 (Andrew et al., 2014). In the experiments, KI-brine was used as the wetting phase and decane as the nonwetting phase. As the samples were clean quarry stones, they were water-wet. Cylindrical samples with a diameter of ≈4 mm and a height of 10 mm were inserted into flow cells. First, a high-quality dry micro-CT scan was obtained. Next, the samples were fully saturated with brine and drained by injecting decane. This ensures that wetting films and layers are present in crevices of the pore space prior to imbibition, which is the case in most realistic reservoir scenarios. Then, imbibition was initiated by injecting brine at very low capillary numbers (1.26 × 10 −9 in the Ketton case, 1.8 × 10 −8 in the Gildehauser case). During imbibition, each sample was partly imaged along its height by fast micro-CT (Ketton: 400 scans at 38 s/scan, voxel size 3.28 μm, image height 3.28 mm; Gildehauser: 40 scans at 45 s/scan, voxel size 4.4 μm, image height 2.49 mm). Imaging was performed at I13-2 (Diamond Light Source) and at TOMCAT (Swiss Light Source). As detailed in the aforementioned references, the images were treated with a nonlocal means filter and segmented with the watershed algorithm in Avizo 9 (ThermoFisher/FEI). The final Gildehauser image contained an oil cluster connected to the top and the bottom of the image, that is, the oil in the sample may still have been mobile.

Quasi-Static Modeling
The pore network models used here were extracted from the segmented dry scans using the maximal ball method (Raeini et al., 2017). This method employs the distance map of the pore space image to find pores and throats as constrictions and dilations. The radii of these elements were found by looking for the maximum inscribed spheres at these locations. Shape factors were used to describe the angularity of pores and throats as described in Bultreys et al. (2018). The pores and throats had circular, triangular, and square cross-sections, based on their shape factor (Valvatne & Blunt, 2004).
After network extraction, imbibition was simulated using a quasi-static invasion percolation approach, based on Øren et al. (1998), Valvatne and Blunt (2004), and Ruspini et al. (2017). Following these works, entry pressures were calculated for each pore and throat which could be accessed by the invading fluid and from which the defending fluid had an escape path (i.e., was not trapped). The pressure calculations are summarized in the supporting information. Quasi-static models do not have a time scale or injection rate associated to them and only predict the order in which pores and throats are invaded. The wetting phase was set to enter through the top and bottom face, and the nonwetting phase was allowed to escape through both, to account for backflow in the porous medium outside the field of view. Fluid clusters which were not connected to either of these faces, either through the centers of pores or throats or through fluid layers in corners of triangular or square network elements, were deemed to be trapped. We used the fluid occupancy at the start of imbibition in the experimental data to assign the initial condition of the model. This was done by mapping the experimental fluid occupancy to the PNM, as explained in section 2.3. The initial capillary pressure at the start of imbibition was set to the invasion pressure of the smallest throat that is oil-filled in the micro-CT image after drainage. We list here the PNM assumptions and user-defined parameters which can influence the simulated filling sequence: • The segmentation of the micro-CT scan affects the PNM extraction. In the simulations, this can cause poor predictions of the entry pressures in pores and throats. In this work, the large pores and the good image quality in the dry scans were deemed to lead to high-quality segmentations. We used the segmented images provided in the public domain data sets, without tuning the segmentation to improve the PNM results. • The detection of pore and throat centers as local minima and maxima of the distance map. The main user-defined parameter is the minimum inscribed radius for a point to be called a pore center, set to 1.5 voxel lengths. This does not affect the results for the well-resolved, rather simple pore structures in this work. In any event, the same network extraction was used to analyze both the experiment and the model. • Pore and throat shape factor determination, and the simplification of pore shapes to (predominantly) triangles with the same shape factor and inscribed radius (Øren et al., 1998). This induces an approximation which may cause deviations in the predicted entry pressures, but it is fully defined by the detection of pore centers and throat surfaces. • Advancing contact angle assignment. These were randomly assigned in each pore or throat in the range between 35 • and 45 • for Ketton (Singh et al., 2016;Scanziani et al., 2017) and 35 degrees and 55 degrees for Gildehauser (Khishvand et al., 2016), based on in situ contact angle measurements. Varying the advancing contact angle influences the relative importance of the different filling processes, thereby impacting the filling sequence and the amount of trapping in the simulations. The sensitivity within likely values (several ranges spanning 25 • to 65 • ) was found not to influence the conclusions in this paper, as shown in the supporting information. • The capillary pressure reached during drainage, prior to imbibition. This is inferred from the radius of the smallest pore or throat that is oil-filled in the micro-CT image after drainage. This impacts the thickness of water layers initially present in the PNM, which can influence the amount of snap-off in the simulations. • One parameter in Ruspini et al.'s cooperative pore-filling algorithm, which determines when the interface becomes unstable while entering the throats. This parameter (set to 0.5 as suggested in the original paper) has a very limited influence on the simulations results (Ruspini et al., 2017).

Data Analysis
To allow consistent comparison of the experiments to PNM simulations, we used the PNM as an image analysis tool on the imbibition data sets. As detailed in Bultreys et al. (2018), the inscribed spheres at each pore and throat center in the dry scan were determined during network extraction. These were used to define the fluid occupying the pore or throat center. For each time step in the imbibition experiment, a pore or throat was called wetting-filled if more than half of the voxels of the associated inscribed sphere are segmented as wetting fluid. The simulations were run on the same PNM used to analyze the experiments, enabling a like-for-like comparison (Figure 1). The measures shown in this paper were calculated on the PNM structure, with fluid distributions mapped from the experiments, or from simulations on the PNM.

Filling Size
Pores and throats in the experiment were considered to be invaded at time step t if they were filled with wetting fluid in image t and with nonwetting fluid in image t − 1. If they changed filling state more than 3 times, they were not listed here, as these were a small number of mainly narrow throats affected by image noise. The filling time was then translated to the filling sequence rank number (the rank of the pore or throat filling event when ordered according to the filling time). This was necessary because the simulations do not have a time scale, but only predict in which order pores and throats fill. To compare the models to the experiments, we calculated the running average of the pore and throat radii over 10 pore/throat fillings in the Ketton case, and over 20 fillings in the Gildehauser case. We also calculated the standard deviation on this running average. In Gildehauser, the rate of filling was relatively large compared to the time resolution of the imaging, resulting in many fillings being detected per time step. Therefore, filling sizes at certain time steps in this experiment were averaged per time step rather than per 20 filling events.

Connectivity
To characterize the connectivity, we calculate the Euler characteristic of the shape formed by the throats (conceptualized here as cylinders) and pores (conceptualized as spheres) that are filled with the nonwetting phase from either the experiment or the simulation. We first remove isolated nonwetting phase-filled throats to reduce the dependency on imaging noise. Then, we modify this shape by adding a sphere to the end of any cylinder which does not already end in a sphere. This does not change the Euler characteristic as it is invariant under continuous deformation of the shape. The Euler characteristic can be calculated using the Betti numbers i : Here, 0 equals the number of connected components and 2 is the number of holes in the shape, the latter by construction equaling zero. 1 is the number of loops in the shape , equalling N e − N v + 0 , with N e the number of edges and N v the number of vertices of the corresponding graph (Berge, 2001). Therefore, the calculation of the Euler characteristic reduces to Note that the calculation does not depend on the exact pore and throat shapes. Detecting the connected components in the nonwetting phase graph allows to calculate the amount of clusters ( 0 ) and their sizes, after which the number of loops can be calculated as 0 − .

Flow Paths
The relevance of discrepancies in the filling sequence to relative permeabilities is a question of key concern . However, the calculation of volumes and fluid flow conductances for each pore and throat often induces dependence to nonphysical tuning parameters (Sorbie & Skauge, 2011). Since we are only interested here in the fluid displacements and not in improving conductivity or saturation calculations, we circumvent this problem by calculating relative permeabilities for both the experimental and the simulated fluid distributions on the same PNM, using the same (possibly imprecise) conductivities and volumes. This yields a flow-weighted connectivity metric which is a proxy to (but not necessarily equal to) physical relative permeability measurements.
For the experiments, relative permeabilities were calculated by mapping the experimentally observed fluid distribution on the same PNM as used in the imbibition simulations. The Hagen-Poiseuille conductivity model was used to calculate the flow rates through the wetting and the nonwetting phase separately, using the OpenPNM pore network modeling package (Gostick et al., 2016). We neglect corner flow conductivity by assuming that throats filled with one phase have zero conductivity for the other phase. Flow rates for the simulations were calculated on the same PNM with the same conductivities, and discrepancies compared to the experiment therefore came solely from differences in the fluid distribution. A good match means the flow distribution in the model is approximately correct to predict the upscaled properties. However, to match experimentally measured relative permeabilities (based on pressure measurements), it may still be necessary to improve the volume and conductivity assignment of pores and throats. This is outside the scope of the study.

Results and Discussion
Previous studies have shown that the pore-scale fluid distribution in experiments is not fully reproducible (Bultreys et al., 2018;Ferrari et al., 2015;Ling et al., 2017). The averaged properties are, however, much more reproducible, as different realizations of the fluid distribution can lead to consistent averaged results if they are statistically similar. Therefore, simply comparing the filling sequence on a pore-by-pore basis does not offer a relevant answer (Bultreys et al., 2018). Recognizing that a good prediction of the upscaled multiphase flow properties needs to capture the geometry and the connectivity of the fluid phases throughout the filling sequence, we use the following criteria to validate the model: • In terms of size (radius), which pores and throats are filled throughout the filling sequence, and where are these pores in the sample? • How does the connectivity of the nonwetting phase change during the filling sequence? • How do the flow paths of the wetting and the nonwetting phase change throughout the filling sequence?
As our goal was to evaluate physical concepts rather than model implementation, the influence of PNM nonuniqueness on the validation was reduced as much as possible. Parameters that affect the filling sequence were set based on physical principles and were sensitivity-checked.

Filling Size and Location
The PNM used in this study explicitly defines pores as local dilations and throats as local constrictions (Raeini et al., 2017). In terms of the associated capillary energy, pores should locally be less favorable places for the invading wetting phase to pervade, while throats should be more favorable. Therefore, we tracked the invasion of the wetting phase in these locations in the experiment and the model, using the same network extraction on both. We define filling in the experiment as events where the image shows a change in occupancy, such that the voxels associated with the maximal sphere centered on a pore or throat first become occupied principally by the wetting phase in the segmented image. Figure 3 shows the radii of pores and throats that were invaded during the experiment as a function of their rank order in the filling sequence. In general, pores and throats of similar sizes were filled throughout the filling sequence in the experiment and the model. In both cases, filling did not strictly happen in order of increasing radius due to the occurrence of the three types of filling (snap-off, piston-like advance and cooperative pore filling) and the variation in pore and throat shapes. A small number of pores and throats with large radii were filled in the experiment, but not in the simulations. This could be related to geometry simplifications in the network model, to local variations in wettability , or to viscous and/or inertial effects which are not captured in the quasi-static PNM. A decrease of the local capillary pressure caused by the latter two effects may have led to the invasion of pores with larger characteristic size than predicted by the quasi-static model. The discrepancy between model and experiment may be larger in the Gildehauser case than in the Ketton case due to the observation of ganglion dynamics in the former experiment. Imaging at higher temporal resolution would be useful to clarify this. Despite the importance large pores can have on saturations and flow, the relatively limited number of wrongly invaded pores did not appear to result in significant deviations in the flow paths, as shown later in this paper.
Aside from the pore and throat filling size, the location of fillings over time can also be compared (Figure 4). In the Ketton experiment, the invasion was more front-like than in the PNM. This is an effect related to viscous forces and flow boundary conditions, which thus cannot be reproduced in the quasi-static model. The PNM invaded pores of similar sizes, yet in different locations. This front-effect was not observed in the Gildehauser case, where filling happened throughout the sample. Figure 5 shows the size distribution of the water-filled pores and throats in the experiments and the model after imbibition, and indicates pores and throats filled with a fluid in the model different from the experiment, similar to Bultreys et al. (2018) and Øren et al. (2019). In the Ketton experiment both model and experiment reached residual saturation. There is a very close agreement, with a filling discrepancy of 9% in the pores and 6% in the throats (the percentage of pores or throats that were filled with a different fluid in the model as in the experiment).
The Gildehauser experiment did not reach the residual saturation in the experimental field of view, and therefore, we compared the simulated fluid distribution with the same saturation as the final state in the experiment. The model overpredicted the number of water-filled pores and throats, resulting in a filling discrepancy of 30% in the throats and 38% in the pores. This is explained by the fact that the simulation was stopped at the same saturation as the experiment, while a number of large pores were filled earlier in the experiment than in the model (Figure 3). A previous study showed filling discrepancies of 18% and 21% for PNM simulations compared to steady-state Bentheimer experiments at the end of imbibition (Bultreys et al., 2018), and a standard deviation of up to 25% on the pore occupancy for the intermediately sized pores in repeated experiments. A different study on Paaratte sandstone found filling discrepancies of 15% and 13% (Øren et al., 2019).
The validation of pore and throat filling sizes presented here improves on the approach in Bultreys et al. (2018) and Øren et al. (2019) because it took the dynamics of the filling process into account. Overall, a more favorable agreement was found here than in Bultreys et al. (2018) due to the use of a different cooperative pore-filling algorithm (Ruspini et al., 2017) and initial and boundary conditions that more closely resembled those of the experiment. The cooperative pore-filling algorithm developed by Ruspini et al. (2017) is a significant improvement over the previously available algorithms because it accounts for the local geometry of the surrounding throats to assign the capillary invasion pressures. In repeated drainage-imbibition experiments on Bentheimer, the standard deviation of the filling state was found to be up to 25% for the most variable pores, while the large-scale behavior did not vary significantly (Bultreys et al., 2018). The filling discrepancy between the PNM simulation and the experiment is thus of similar magnitude to the difference between two repeat experiments that are macroscopically nearly indistinguishable.

Connectivity
The upscaled flow characteristics are strongly influenced by the phase connectivity (Hilfer, 2006), which can be characterized by the Euler characteristic Schlüter et al., 2016). We compared the evolution of the Euler characteristic during imbibition in the model and the experiment ( Figure 6). A consistent comparison to the simulations was achieved by mapping the experimental distributions on the PNM and calculating the Euler characteristic according to the method described in section 2.3. The experimental trend is consistent with direct image-based calculations (Alpak et al., 2018). Figure 6 shows that magnitude and trends of the Euler characteristic are correctly represented by the model. There are, however, some discrepancies. In the Ketton simulation the model shows a local minimum at the start of imbibition. This is caused by an image-processing artifact: small pores that were wrongly identified as nonwetting-filled after drainage. The effect of this on the Euler characteristic was partly compensated by redundant loops in the modeled nonwetting phase distribution being removed too quickly. The latter may have been caused by a mismatch in the amount of cooperative pore filling compared to snap-off between the simulation and the experiment, due to geometric simplifications in the model. Imbibition in Ketton is more sensitive to cooperative pore filling than most other rock samples, due to Ketton's pore structure being similar to that of a bead pack (Ruspini et al., 2017). The relatively small variability in pore size and low constriction ratios in such samples increases the importance of the occupancy in a pore's neighbors, whereas the magnitude of pore and throat radii may have a more dominant role in other samples. In the Ketton case, the model may therefore overestimate the amount of capillary trapping.
In the Gildehauser case, the Euler characteristic increased more quickly in the experiment than in the model. In the experiment, there was a change in the slope of the Euler characteristic and the number of loops after ≈1,600 pores and throats were filled. This was likely related to the occurrence of two different filling regimes observed by Rücker et al. (2015): first, pathway flow and film swelling, and then a regime which includes cluster dynamics. This was not reproduced in the model. Furthermore, the faster increase in the experimental Euler characteristic is partly due to misidentification of the fluid occupancy in smaller pores and throats caused by imaging noise in the Gildehauser experiment. This resulted in an excessive amount of small nonwetting phase clusters being detected. However, the average sizes of multithroat clusters (excluding single-throat clusters to reduce imaging noise) were in fair agreement between the model and the experiment, as were the sizes of the largest nonwetting phase clusters which dominate the flow (see supporting information). As will be shown in the next section, the discrepancies discussed here had a limited influence on the flow paths through each fluid.

Flow Paths
The filled pore/throat sizes and the connectivity together control the flow paths which are accessible to each fluid and consequently the relative permeability. We compared the flow paths in the model and the experiment by calculating relative permeabilities on the imaged and the modeled fluid distributions. This relative permeability can be seen as a flow-weighted connectivity of the fluid distributions. Since both the experimental and the modeled results used the same conductivity and volume assignments in pores and throats, mismatches in the presented relative permeabilities can only come from discrepancies in the modeled fluid distributions.
The resulting relative permeabilities indicated a close match between the model and the experiment ( Figure  7). This implies that the quasi-static model provided a good estimation of the fluid distribution with respect to its impact on the upscaled flow properties. The discrepancies in Ketton may be related to the front-like filling in the experiment, which was not reproduced by the PNM. This, together with a small mismatch of cooperative pore filling versus snap-off and the unknown boundary conditions in the experiment, may have also induced discrepancies near the end of imbibition. While ganglion dynamics have been seen to play an important role near the residual nonwetting phase saturation in the Gildehauser experiment (Rücker et al., 2015), there is no direct evidence of this in the Ketton data set. We do not draw conclusions about the behavior near the end point of imbibition in the Gildehauser case here, as the experiment did not continue until residual nonwetting saturation.

Conclusions
A key question concerning imbibition in permeable materials is how to formulate conceptual models that decrease the degrees of freedom compared to an explicit description, while still capturing the salient behavior of the process. Quasi-static PNMs abstract the process by simplifying the local pore shapes and by assuming that capillary forces are dominant, allowing imbibition to be described as a sequence of simple pore and throat filling processes. However, one of the main problems has been validation: the complex relation between the fluid occupancy of the system and the experimentally measured properties (e.g., relative permeability) typically used for validation makes it difficult to understand the origin and generality of any observed (mis-)match (Bondino et

10.1029/2019WR026587
The analysis we have developed here allowed the comparison of fluid distributions measured at the pore scale to a quasi-static PNM with minimal influence from anything but the fluid distributions. This indicated that the quasi-static assumption provided a good approximation of the order in which pores with different sizes were invaded. The model did misinterpret the filling behavior of some larger pores, and overestimated the amount of loop-breaking in the nonwetting phase. However, the pores and throats causing these discrepancies showed to be of minor importance for the flow paths, which closely matched the experiments.
The current study only focuses on the basic capability of quasi-static PNM to predict fluid distributions in relatively simple cases. The conclusion is that an abstracted, conceptually simple model based on quasi-static physics can give a good approximation of continuum scale properties such as relative permeability, at least in strongly wetted rocks with relatively homogeneous pore structures and at low capillary numbers. Whether or not it does, depends on model implementation (e.g., volume and conductivity assignment) and on accurate input concerning pore geometry and advancing contact angles. These issues remain the most important challenges to improve PNM as generally applicable predictive models. This is especially important in rocks with more complex pore spaces (e.g., higher degree of diagenetic alteration), where resolving the pore space is more challenging. There remains a general need to understand and improve the limitations with regard to image quality and resolution. The validation workflow should thus be extended to more complex rock types and possibly to multiscale pore network models before more general conclusions on the predictiveness of PNM can be made.
In summary, PNMs as a concept appear a valuable starting point for building conceptual models. However, it must be stressed that more validation needs to point out whether this holds in circumstances where more influence of ganglion dynamics has been observed or can be expected, for example, near the end of imbibition (Rücker et al., 2015), in intermediate/mixed-wet samples (Zou et al., 2018) and in samples with very wide pore size distributions.

Data Availability Statement
The experimental data in this manuscript are freely available via the references in section 2. The simulation data are made available through the Open Science Foundation (https://osf.io/h6sqt/?view&urluscore; only=10a17e75ba78461da4e84fbb20f65ee1).