The Origin of Non-thermal Fluctuations in Multiphase Flow in Porous Media

Core flooding experiments to determine multiphase flow in properties of rock such as relative permeability can show significant fluctuations in terms of pressure, saturation, and electrical conductivity. That is typically not considered in the Darcy scale interpretation but treated as noise. However, in recent years, flow regimes that exhibit spatio-temporal variations in pore scale occupancy related to fluid phase pressure changes have been identified. They are associated with topological changes in the fluid configurations caused by pore-scale instabilities such as snap-off. The common understanding of Darcy-scale flow regimes is that pore-scale phenomena and their signature should have averaged out at the scale of representative elementary volumes (REV) and above. In this work, it is demonstrated that pressure fluctuations observed in centimeter-scale experiments commonly considered Darcy-scale at fractional flow conditions, where wetting and non-wetting phases are co-injected into porous rock at small (<10−6) capillary numbers are ultimately caused by pore-scale processes, but there is also a Darcy-scale fractional flow theory aspect. We compare fluctuations in fractional flow experiments conducted on samples of few centimeters size with respective experiments and in-situ micro-CT imaging at pore-scale resolution using synchrotron-based X-ray computed micro-tomography. On that basis we can establish a systematic causality from pore to Darcy scale. At the pore scale, dynamic imaging allows to directly observe the associated breakup and coalescence processes of non-wetting phase clusters, which follow “trajectories” in a “phase diagram” defined by fractional flow and capillary number and can be used to categorize flow regimes. Connected pathway flow would be represented by a fixed point, whereas processes such as ganglion dynamics follow trajectories but are still overall capillary-dominated. That suggests that the origin of the pressure fluctuations observed in centimeter-sized fractional flow experiments are capillary effects. The energy scale of the pressure fluctuations corresponds to 105-106 times the thermal energy scale. This means the fluctuations are non-thermal. At the centimeter scale, there are non-monotonic and even oscillatory solutions permissible by the fractional flow theory, which allow the fluctuations to be visible and—depending on exact conditions—significant at centimeter scale, within the viscous limit of classical (Darcy scale) fractional flow theory. That also means that the phenomenon involves both capillary aspects from the pore or cluster scale and viscous aspects of fractional flow and occurs right at the transition, where the physical description concept changes from pore to Darcy scale.

Core flooding experiments to determine multiphase flow in properties of rock such as relative permeability can show significant fluctuations in terms of pressure, saturation, and electrical conductivity. That is typically not considered in the Darcy scale interpretation but treated as noise. However, in recent years, flow regimes that exhibit spatio-temporal variations in pore scale occupancy related to fluid phase pressure changes have been identified. They are associated with topological changes in the fluid configurations caused by pore-scale instabilities such as snap-off. The common understanding of Darcy-scale flow regimes is that pore-scale phenomena and their signature should have averaged out at the scale of representative elementary volumes (REV) and above. In this work, it is demonstrated that pressure fluctuations observed in centimeter-scale experiments commonly considered Darcy-scale at fractional flow conditions, where wetting and non-wetting phases are co-injected into porous rock at small (<10 −6 ) capillary numbers are ultimately caused by pore-scale processes, but there is also a Darcy-scale fractional flow theory aspect. We compare fluctuations in fractional flow experiments conducted on samples of few centimeters size with respective experiments and in-situ micro-CT imaging at pore-scale resolution using synchrotron-based X-ray computed micro-tomography. On that basis we can establish a systematic causality from pore to Darcy scale. At the pore scale, dynamic imaging allows to directly observe the associated breakup and coalescence processes of non-wetting phase clusters, which follow "trajectories" in a "phase diagram" defined by fractional flow and capillary number and can be used to categorize flow regimes. Connected pathway flow would be represented by a fixed point, whereas processes such as ganglion dynamics follow trajectories but are still overall capillary-dominated. That suggests that the origin of the pressure fluctuations observed in centimeter-sized fractional flow experiments are capillary effects. The energy scale of the pressure fluctuations corresponds to 10 5 -10 6 times the thermal energy scale. This means the fluctuations are non-thermal. At the centimeter scale, there are non-monotonic and even oscillatory solutions permissible by the fractional flow theory, which allow the fluctuations to be visible and-depending

INTRODUCTION
Multiphase flow in porous media is an integral part in many aspects of every-day life and plays a critical role in some of the most important processes and technologies from agriculture to energy. Examples range from the hydrology in the Vadose zone where the water table rises and falls and water and air occupy the pore space in soil, contaminant hydrology, hydrocarbon recovery, hydrogen storage, and carbon capture and sequestration (CCS) (Bui et al., 2018). Multiphase transport is also a rate-limiting step in gas diffusion layers in electrocatalytic devices such as fuel cells (Simon et al., 2017), electrolysis, and more novel concepts, where CO 2 is converted into base chemicals (Kondratenko et al., 2013).
In most applications, multiphase transport is described with the 2-phase Darcy equations, which are a continuum mechanic concept for relating transport, i.e., average phase fluxes to average pressure gradients applicable at the "Darcy scale." They are phenomenological extensions of Darcy's law from single to multiphase flow. One of the consequences of operating with a phenomenological transport equation is that it contains parameters, such as relative permeability, which cannot be predicted within the framework of the 2-phase Darcy equations. Since relative permeability is specific to the porous medium and its chemical interaction with the fluids present manifesting in characteristic wetting behavior of the system (Abdallah et al., 2007), they have to be determined for each situation individually. Common methods to experimentally determine relative permeability are core flooding experiments, which are typically conducted on porous media samples of few centimeters in size. For instance, in the petroleum industry, cylindrical rock samples from drilled cores are used, which are limited in size by the diameter of the core and are typically between 2.54 and 5.00 cm in length and diameter, which is commonly assumed to represent Darcy scale. The steady-state method is one of the most trusted methods as it provides a wide accessible saturation range and in combination with numerical simulation (Kokkedee et al., 1996;Masalmeh et al., 2014;Sorop et al., 2015) allows the correct treatment of experimental artifacts, such the capillary end-effect (Huang and Honarpour, 1998). In steady-state experiments, the two immiscible fluids phases, e.g., water or brine and oil or gas are co-injected at stepwise varied fractions of wetting phase flux over total flux, termed fractional flow f w . At each fractional flow step, fluids are co-injected until a "steady-state" has been reached, which is defined by saturation, pressure-drop, and electrical resistivity becoming time-independent, i.e., dS w /dt = 0, where S w is the space and time-averaged saturation. Relative permeability is then determined from the time-averaged pressure-drop and time-averaged saturation.
However, many of these fractional flow experiments show notable fluctuations (Datta et al., 2014a,b;Masalmeh et al., 2014;Reynolds and Krevor, 2015;Rücker et al., 2015a;Sorop et al., 2015;Gao et al., 2017Gao et al., , 2019Gao et al., , 2020Lin Q. et al., 2018;Alcorn et al., 2019;Clennell et al., 2019;Lin et al., 2019a;Spurin et al., 2019;Wang and Masalmeh, 2019;Menke et al., 2021), for instance in pressure drop but also saturation. In some cases, the magnitude of these fluctuations is comparable or even larger than the average values of the respective property, e.g., pressure drop, between the previous or following fractional flow. Historically, such fluctuation have been considered as indicative for failed experiments and often led to dismissal of the experiment. That is potentially the reason why despite extensive literature on steadystate relative permeability in only relatively few cases the raw data is shown. Even in the cases where pressure drop as a function of time is shown, it is often not clear whether this data is filtered or smoothed, operating under the assumption that the fluctuations are caused by noise.
Also, pressure fluctuations are more prominent for water-wet rock and much less visible for intermediate or mixed-wet cases (Jung et al., 2016;Lin Q. et al., 2018;Lin et al., 2019a), where the magnitude of capillary pressure is much less or even near zero (Lin et al., 2019a). Experiments conducted on twin samples of the same rock to exclude all other factors show noticeable fluctuations in the water-wet case (Lin Q. et al., 2018) but much less in the intermediate-wet case (Lin et al., 2019a), which has been achieved by aging with crude oil. In the intermediate-wet case, capillary pressure is near zero, because the mean curvature is near zero, which has been established from imaging the liquidliquid interfaces by micro-CT. But that does not mean the interfaces are flat, but rather form bi-continuous interfaces that allow both aqueous and oil phases to be continuous at the same time, i.e., fractional flow does not cause displacement events, which is the primary reason why pressure fluctuations are less. Wettability heterogeneity may also impact the magnitude of fluctuations due to the associated energy dissipation (Murison et al., 2014). Also, fluctuations are stronger for a gas-liquid system and more often reported explicitly (Alkan and Müller, 2008;Reynolds and Krevor, 2015;Xu et al., 2015;Spurin et al., 2020), which suggests that viscosity ratio may be an important factor, which has been also reported in (Spurin et al., 2019).
Pore scale experiments (DiCarlo et al., 2003;Moebius and Or, 2012;Berg et al., 2013Berg et al., , 2014Armstrong et al., 2014a;Reynolds et al., 2017;Lin Q. et al., 2018;Lin et al., 2019a) and numerical simulations (Lenormand et al., 1983;Raeini et al., 2014;Armstrong et al., 2015;Guédon et al., 2017;Alpak et al., 2019;Berg C. F. et al., 2020;Winkler et al., 2020) also exhibit fluctuations in pressure and saturation (Ramstad and Hansen, 2006;Pak et al., 2015), which are caused by pore scale displacement events, such as Haines jumps and coalescence (Rücker et al., 2015b), where the non-wetting phase replaces the wetting phase and snap-off and pistonlike displacement (Lenormand et al., 1983;Dixit et al., 1998), where the wetting phase replaces the non-wetting phase. These events lead to interruption and rearrangement of the connected pathways the respective phases flow through (Tuller and Or, 2001) and are described by a rigorous theoretical framework of pore-scale thermodynamics (Morrow, 1970). Depending on flow rate and other conditions (Avraam and Payatakes, 1995a;Lenormand and Touboul, 1998), there are different regimes that can cause fluctuations in pressure and saturation that have been systematically investigated in Avraam and Payatakes (1995a), Tsakiroglou (2019). The main regimes identified were connected pathway flow, ganglion dynamics and drop traffic. While sequences of oil-filling and water filling events of different nature cause ganglion dynamics (Rücker et al., 2015b) already at capillary numbers below Ca < 10 −6 , effects such as ganglion dynamics, where such events occur repeatedly at the same location (Lenormand et al., 1983;Datta et al., 2014a,b;Gao et al., 2017Gao et al., , 2019Gao et al., , 2020Spurin et al., 2019), are often observed close to the critical capillary number for capillary de-saturation (Berg and Cense, 2009) and likely caused by viscous mobilization and associated breakup, which can in the extreme case lead to "drop traffic" (Avraam and Payatakes, 1995a,b). Note that in the porous media literature, the term "intermittency" is frequently used instead of ganglion dynamics. However, outside the porous media community, this term is commonly associated with turbulent flow, which could lead to misconceptions as in porous media, and this behavior is often observed already in the capillary dominated regime. Periodic pressure fluctuations can also be caused by snap-off of non-wetting phase bubbles at the outlet although this is a mechanism observed more in spontaneous imbibition (Unsal et al., 2007a,b;Unsal et al., 2009;Bartels et al., 2019). While the exact pore scale mechanism can differ, the pressure signature is typically caused by capillary effects in a capillary-dominated regime . A recent pore scale simulation study that shows such fluctuations caused by pore-scale displacement events introduced the term "athermal" to characterize the energy scale of these fluctuations (Winkler et al., 2020).
The referenced literature suggests that these non-thermal fluctuations are a relatively universal phenomenon occurring in a very wide range of rock types and are observable on different scales. The key step that has been missing so far is a direct link between pore-scale fluctuations and Darcy scale behavior, which is the focus of this paper. This link between fluctuations observed in core flooding experiments at a length scale of several centimeters and their origin from pore scale displacement events is important because from both a very fundamental perspective and also for practical reasons. It is important to understand whether fluctuations are caused by pore scale displacement events, because such events are always to some extent irreversible and associated with the dissipation of energy (Morrow, 1970;Seth and Morrow, 2007;Berg et al., 2013). The big question is whether this energy dissipation by pore scale events is correctly captured when using time-averaged pressure and saturation data. At this moment, we do not fully understand the consequence of the averaging in case that fluctuations are not just noise but indeed caused by pore scale displacement events. Also, for the interpretation of such experiments by inverse modeling to determine the relative permeability (Berg S. et al., 2020), for a correct assessment of the associated uncertainty, it is important to understand whether fluctuations are noise or have another cause, e.g., a physical mechanism.
In this work we address this question of the link between fluctuations and pore-scale displacement events by analyzing the magnitude and statistics of Darcy scale "steady-state" fractional flow experiments for different rocks and fluid/fluid/solid wetting and comparison with pore-scale fractional flow experiments with synchrotron beamline-based fast x-ray computed tomography. We will discuss the reason why fluctuations involve hundreds or more pores. From a pore scale perspective alone, this would not be very obvious, because fluctuations are mainly associated with individual pore scale displacement events and larger fluctuations are interpreted from the pore scale perspective as ganglion dynamics. The dominance of fluctuations involving hundreds to thousands of individual pores becomes understandable as in the transition to the Darcy scale fluctuations can only persist if a permissible (Darcy-scale) fractional flow solution exists, for which we will present supporting evidence. That emphasizes the need to presenting both the pore scale and the Darcy scale perspective on the same phenomenon. Each experiment has significant complexity, which cannot be fully covered in this work, which focuses very much on the link between pore and Darcy scale. Therefore, we largely rely on thoroughly conducted experiments on each scale using dedicated experimental setups that are designed to rule out instrumental artifacts, i.e., avoid back pressure controllers and other potential source of fluctuations other than the ones caused by the multiphase flow in porous media. We make use of the best experimental data that are unfortunately not conducted on exactly the same sample. We use Fontainebleau sandstone for Darcy scale experiments and sintered glass of comparable pore size at the pore scale. The phenomena found in the sandstone sample and the sintered glass are comparable, which is also established through consistency with other studies in the literature, which suggest that phenomena are actually generic and observed in many different rock types. The experiments we present here have the highest level of confidence to rule out such artifacts, which is why we gave priority over fully consistent samples. Ultimately, the pressure signature but also associated fluctuations in connectivity (McClure et al., 2016) (which in the Darcy scale experiments is represented by electrical resistivity) (Liu Z. et al., 2018) provides the connection between Darcy scale fluctuations and pore scale displacement events.

METHODS AND MATERIALS
Two different set-ups were used in this study to relate pore-scale mechanisms with Darcy-scale responses. The different set-ups, sample preparation and analysis are described below.

Rocks and Fluids
For the Darcy scale experiments, a cylindrical Fontainebleau sandstone sample with a diameter of d = 2.54 cm and a length of L = 4.45 cm was used. Fontainebleau sandstone is a fine-grained, well-sorted cretaceous sandstone outcropping in South France, well-known for its regular grain framework, the almost monomineral composition up to 99.5% quartz and for the absence of clay in it (Jacquin, 1964;Al Saadi et al., 2017). The rock sample had a porosity of φ = 0.13 and a permeability of K abs = 5.32 × 10 −13 m 2 (539 mD). This permeability is particularly well-suited for special core analysis (SCAL) steady-state flow experiments, because it falls in the middle of the operating envelope of the setup, which has a maximum accessible permeability of about 2-3 Darcy. Before the experiment the rock sample was cleaned in a chloroform-methanol azeotrope in a Soxhlett extractor to complete water wetness. The sample was not aged with crude oil and only inert fluids were used, which means that the rock remained strongly water wet. In a second series of Darcy scale experiments addressing the effect of wettability through variation of fluids, 3 cylindrical rock samples of Bentheimer sandstone rock (Lin et al., 2019a,b) with 5 cm length and 3.8 cm diameter were used. The porosity ranged between 0.246 and 0.248 and the permeability between 2,577 and 2,642 mD. The samples were cleaned with by Soxhlett extraction and then saturated with brine for permeability measurements. The three rock samples were then de-saturated with 3 crude oils (basic properties are listed in Table 1 and analysis results of total acid number (TAN), total basic nitrogen (TBN) and SARA are listed in Table 2) and subsequently aged (Lin et al., 2019b) for 40 days. Based on an extensive analysis of the crude oils, the potential for wettability alteration is highest for crude 3, somewhat less for crude 2 and least for crude 1, i.e., it is expected that for crude 1 the Bentheimer rock remains more water wet, whereas for crude 3, the rock is expected to become more intermediate wet.
For the micro-CT (pore-scale) flow experiments, a sintered glass model purchased from Robuglass (Berg et al., 2014) was chosen because of its large pores that are particularly well-suited for imaging pore scale fluid distributions, i.e., much more suitable than the very narrow pores of Fontainebleau sandstone. The sample had a diameter of Ø = 4 mm and a length of L = 2 cm. The porosity determined from a segmented micro-CT image was φ = 0.36. The permeability of K abs = 1.38 × 10 −11 m 2 (14 D) was obtained from numerical flow simulation with the Navier-Stokes solver from GeoDICT (Math2Market) (Berg and van Wunnik, 2017). This permeability is too high for the SCAL steady-state experimental setup (the associated pressure drops are too small for the range of the pressure transducers). The sample is naturally (from the sintering process) clean and strongly water wet. The fluids used in both experiments was n-decane and brine doped for contrast with 14.3% CsCl in the pore-scale experiment and 5% in the Darcy scale experiment. Density ρ and viscosity µ of the fluids are listed in Table 1. The interfacial tension was around 40 mN/m. Decane and brine don't form any emulsions, which might hypothetically also cause pressure fluctuations (Cobos et al., 2009;Yu et al., 2019).

"Darcy Scale" Fractional Flow Experiments
The Darcy scale experiments were conducted in an experimental setup, which is routinely used for measuring relative permeability with the so-called "steady-state" method, where water and oil phases are co-injected at fractional flow: where q α are the fluxes and λ α = k r,α /µ α the mobility of phase α, with the relative permeability k r,α and viscosity µ α (α = w for wetting phase and α = n for the non-wetting phase). The co-injection at f w is realized by two dual-piston Quizzix pumps (Chandler, Metek) that operate with zero-dead volume switching valves for smooth and pulse-free continuous flow (Vindum Engineering). Fluids are re-circulated (after separation) without any back-pressure controller that could cause artifacts such as pressure oscillations. The core is enclosed in an X-ray transparent Hassler-type core holder made from carbon fiber epoxy composite that allows in-situ saturation monitoring by a self-built X-ray transmission system that measures saturation at 22 positions along the core (Cense et al., 2014). During the experiment, the pressures drop dP w is monitored at a 1 min time interval and saturation S w and resistivity at a time interval of 8 min 30 s. The conceptual overview of the setup is displayed in Figure 1A.
At the inlet, the oil and water supplied by 1/16 inch tubing are first dispersed by a set of diffusor plates with two rings of concentric circular grooves followed by a metal mesh (both at same cross section as the core) before entering the core. A similar end piece construction is used at the outlet. In the Darcy-scale experiments oil and brine are co-injected at varying flow rates and fractional flow f w with 2 dual-piston Quizzix pumps that operate with zero-dead volume switching valves for smooth and pulse-free continuous flow. Fluids are re-circulated (after separation) without any back-pressure controller. Saturation is measured in-situ during the experiment using a linear X-ray scanner. The micro-CT setup (B) uses a core holder with 2 integrated micro-pumps  for water and oil co-injection at specified fractional flow f w . At the outlet side the fluids are collected at ambient pressure, i.e., no back-pressure control system is used. The porous sample is wrapped with paper (for a hydrophilic boundary condition) and heat-shrunk into a poly-carbonate cylinder. The setup is self-contained as it is rotating during the micro-CT scans, i.e., no external flow lines are used. Power is supplied via slip rings and the pressure sensor readings are transmitted with a wireless system.
Starting with primary drainage, the fractional flow f w was systematically changed from 100% brine (f w = 1) to 100% oil (f w = 0), followed by 1st imbibition cycle, where the fractional flow f w was systematically changed from 100% oil (f w = 0) to 100% brine (f w = 1) in 10 saturation steps each. The total injection rate of q w + q o was thereby kept constantly at 3 ml/min. That corresponds to a capillary number of Ca = 10 −5 , which is still smaller than the onset of capillary de-saturation in imbibition and 1-2 orders below the onset of capillary de-saturation in drainage (Berg and Cense, 2009) but also high enough to suppress any capillary end-effect (Huang and Honarpour, 1998).
At each saturation step f w , the saturation and phase pressures were recorded when, based on a constant pressure with a variation of the respective average being <1%, steady state was reached. The pressure readings were obtained with differential pressure transducers (Rosemount, type 3051 with an accuracy of ±0.15%) (Kokkedee et al., 1996). In addition, conductivity and saturation profiles at 30 positions along the rock sample were measured. The saturation profiles were obtained by Xray transmission measurements calibrated to both 100% brine and oil saturation. The X-ray transmission is measured using a home-built linear X-ray scanner with a conventional X-ray tube and semiconductor-based detector, which makes individual point measurements (after an 8 min 30 s average to achieve sufficient signal-to-noise ratio). For details on the equipment, sample preparation, measurement procedures and analysis were published in previous studies (Kokkedee et al., 1996;Berg et al., 2016;Berg S. et al., 2020). For the particular case of the strongly water-wet Fontainebleau sandstone sample, the saturation profiles along the sample S w (x) (not shown) are flat as expected for a homogeneous sample and there is no noticeable capillary end-effect. The absence of saturation gradients allows us to directly compute the relative permeability-saturation relationships using the two-phase extension of Darcy's law for the phase fluxes where k r,α is the relative permeability of phase α, K is the absolute permeability of the porous medium, µ α is the viscosity and p α the pressure of phase α over sample length L. Resistivity is measured by a 4-point impedance measurement technique using (blackened platinum) electrodes embedded into the end pieces of the core holder.

Micro-CT Beamline Experiments (PS)
The experiments were performed in 2012 at the TOMCAT beamline of the Swiss-Light Source Synchrotron at the Paul-Scherrer Institute, Villigen, Switzerland. The flow cell was specially designed for fractional flow experiments with two remotely controlled micro piston pumps (based on linear motors from Physics Instruments with a modified gear box for low-rate injection) enabling pulse-free continuous, low rate injection while performing 3D imaging Rücker et al., 2015a). The pumps are built-into the core holder itself, which rotates continuously during the fast tomography in order to avoid any displacements caused by bending of flow lines. Flow lines of 0.8 mm diameter directly enter the bottom of the core, which has-because of lack of space-no extra diffusor plate. Since we only monitor the middle section of the core, the first section of the core acts as diffusor. The pressure of water and oil phases at injection side was monitored with 2 pressure sensors (Keller 2 MI), which transmitted the pressure data in real-time using a wireless method enabling a constant data acquisition in the field of view during the experiment. The conceptual overview of the setup is displayed in Figure 1B.
One full 3D image with a voxel size of 2.2 µm was recorded in 1 min. The reconstructed micro-CT images were further processed and segmented with the software package Avizo (ThermoFisher). Details on the flow cell, sample preparation, the experimental settings and image processing were described previously (Armstrong et al., 2014a;Rücker et al., 2015a). Further analysis was conducted with MatLab (MathWorks).
Even though at a resolution of 2.2 µm process such as film flow are not resolved anymore, such sub-resolution physics is effectively captured through the effective displacement events leading to filling of pores or changes in non-wetting phase connectivity, which is captured at this resolution.
The flow sequence included 3 fractional flow steps f w (f w = 0.8, f w = 0.5 and f w = 0.2), where at each f w the total injection rate q tot = q w + q o was increased step-wise from 3 µl/min to 30 µl/min and 300 µl/min. The associated capillary numbers, which will be shown later in more detail in Figure 13, are in the capillary-dominated regime below the critical onset for viscous mobilization (Oughanem et al., 2015).
The permeability of the segmented pore space (in the imaged field of view) K abs was computed by using a Stokes flow simulation conducted with GeoDICT (Math2Market). In addition, for selected time sequences of 2-phase flow experiments, the permeability of water K w and oil phases K o was computed in a similar fashion. The respective connected pathway relative permeability for water k r,w and oil k r,o was then obtained by dividing by the absolute permeability, i.e., k r,w/o = K w,o /K abs .

RESULTS AND DISCUSSION
One of the main findings of the Darcy scale experiment is that the raw data contains fluctuations which-while varying in amplitude significantly from experiment to experiment-have often a systematic trend with fractional flow showing often maximum amplitudes at intermediate fractional flows.
In section Fluctuations in "Darcy scale" Fractional Flow Experiments, a detailed analysis of the noise structure in Darcyscale steady-state core flooding experiments is presented. The focus is less the resulting relative permeability, which is also shown for completeness, but rather the raw data itself. We make a clear distinction between instrumental noise with Gaussian distribution (noise as a random process would result in a Gaussian probability distribution, fully described by its mean and standard deviation) and fluctuations with much larger amplitude than the instrumental noise for which we can largely rule out any instrumental artifacts. In section Dependency of the Fluctuations on Wettability the effect of wettability on the observed pressure fluctuation is discussed.
Finally, in section Energy Scale of Fluctuations, the energy scale of the large fluctuations observed in the Darcy scale flow experiments is analyzed and identified as a capillary energy scale, which is orders of magnitude larger than the thermal energies, i.e., fluctuations are non-thermal and of pore-scale origin. In section Time Scale of Fluctuations the time scale of the fluctuations is put into perspective with other characteristic times scales encountered in multiphase flow.
In Flow Experiments. This is another piece of evidence that the significant fluctuations observed in SCAL experiments are defacto caused by pore scale events. Using the insight of the insitu 3D imaging, in section Pore Scale: Breakup and Coalescence Processes -Ganglion Dynamics the underlying fundamental pore-scale mechanism of breakup and coalescence that constitute ganglion dynamics and cause such events are discussed. In section Identification of Pore Scale Flow Regimes by Collective Cluster Dynamics and the Pressure/Saturation Signature of Ganglion Dynamics the effective flow regimes are characterized on the basis of saturation and macroscopic capillary number. The pore scale events discussed in section Pressure and Saturation Fluctuations in Micro-CT Fractional Flow Experiments suggest that ganglion dynamics events are initiated by an imbibition cycle followed by a drainage cycle.

Differentiation of Capillary Fluctuations From Gaussian Noise
Fractional flow experiments, often also referred as the steady-state method, are routinely used to determine relative permeability from cylindrical rock samples. Two immiscible fluids, such as oil and water are co-injected at fractional flow f w FIGURE 2 | Water pressure drop dP w , water saturation S w and electrical conductivity during a SCAL fractional flow experiment to determine relative permeability on a centimeter-sized Fontainebleau sandstone sample. Overview of the water pressure drop for all fractional flows f w (A) and fluctuations (with base line subtracted) for Frontiers in Water | www.frontiersin.org FIGURE 2 | selected fractional flow vs. time for pressure (blue) and saturation (red), fluctuation histogram and associated power law spectrum (B-E). While the fractional flows f w ≥ 0.75, the fluctuation histogram follows a Gaussian shape (dotted light blue line in the middle in B-E) suggesting instrumental noise as cause. At 0.67 ≥ f w ≥ 0.02, the noise amplitude is significantly larger than instrumental noise and the histogram is clearly non-Gaussian and in most cases also non-symmetric and even bi-modal (D,E) suggesting other causes than instrumental noise. For some fractional flows, e.g., f w = 0.1 the pressure fluctuations are almost perfectly periodic (D) with an oscillation period of t ≈ 21 min, which is also reflected in a peak in the Fourier spectrum (marked with a blue arrow). Similarly, for S w also perfectly periodic behavior can be observed, although at different frequencies (marked with red arrow), which could be potentially caused by the different sampling rates of dP (1 min intervals) and S w (8 min 30 s intervals).
(Equation 1) following a sequence of fractional flow steps that starts at f w = 1 (only brine) to f w = 0 (only oil) for a drainage experiment and vice versa for imbibition. In Figure 2 the raw data of the experiment, i.e., the pressure drop p, saturation S w , and electrical conductivity are plotted for a drainage experiment in the Fontainebleau sandstone sample. Figure 2A shows an overview of the whole experiment, where f w , pressure drop dP w , saturation S w , and resistivity are plotted for a constant flow rate drainage experiment as a function of injected pore volumes (bottom axis) and time (top axis). Pressure drop dP w , saturation S w , and resistivity show fluctuations that are very small for high fractional flows but increase for intermediate fractional flows significantly. For a more detailed view and analysis for each fractional flow, the fluctuations were isolated by subtracting the base line (using a moving average window). In the left panel of Figures 2B-E for 4 selected fractional flows, the isolated fluctuations for pressure drop dP w and saturation S w are displayed.
In the middle panel of Figures 2B-E, the fluctuation histogram is shown. For fractional flows f w ≥ 0.75 (which is not a hard cut-off but rather based on the observations in the fractional flow sequence that has been established over several years as standard protocol), the histogram is symmetric and in good approximation of Gaussian shape, see also  For these cases, the fluctuations are highly periodic, which is also reflected in the peak in the Fourier transform of the fluctuation time signal in the right panel of Figure 2. The fluctuation spectra for saturation (and electrical conductivity data in Figure 2A) differ from the pressure fluctuations, which could be a consequence of the different sampling rates. The pressure is sampled at 1 min time intervals, which means that fluctuations at a 2 min time interval are resolved. Saturation and resistivity are only monitored at an 8 min 30 s time interval. Faster fluctuations are not explicitly resolved.
When performing time averages of the detected signal, a smoothening operation commonly applied during the analysis of multiphase flow experiments in porous media, the fluctuations, which are initially clearly non-Gaussian and even bi-modal as in Figure 2E, become increasingly symmetric as shown in Figure 4.
Here, a rolling average over a window of 2, 4 and 16 individual data points was used. After sufficient averaging, the histogram ( Figure 4B) becomes perfectly Gaussian (which is a consequence of the central limit theorem) and the information this signal may have originally contained get lost, which has implications for the interpretation of the energy dissipation (McClure et al., 2020a,b).
For f w < 0.67, the pressure fluctuations δ p are much larger than what would be expected from the behavior of the pure noise from instrumentation, which has been established based on amplitude and Gaussian shape of the histogram in Figure 3. The pure instrumental noise can be established from f w = 1.0, which corresponds to the injection of 100% brine into a fully brinesaturated rock without any oil, representing single-phase flow conditions. Nevertheless, the instrumental noise of the pressure transducer is proportional to the magnitude of the pressure. With δ p ∝ p we can extrapolate a noise trend from the first 6 f w , which is illustrated by the dotted line in Figure 5A. In Figure 5B a similar approach was applied to saturation S w (Cense et al., 2014). This way the fluctuations can be separated from the noise. Fluctuations falling into the noise regime (below the dotted lines in Figures 5A,B) can be either instrument noise or pore-scale events. But for the events with standard deviation significantly above that line, instrumentation noise can be ruled out entirely.

Differentiation of Capillary Fluctuations From Experimental Artifacts
As Figures 2-5 illustrate the strong pressure fluctuations do not occur for all fractional flows but primarily for fractional flows 0.5 < f w ≤ 0.02 as shown in Figure 5C, where the presumably both fluid phases are mobile and ganglion dynamics (Rücker et al., 2015b) is possible.
In order to characterize the periodic aspect of the fluctuations, the three largest peaks in the Fourier spectrum are displayed in Figures 5C,D. Note that in order to separate dominant peaks with significant amplitude from insignificant peaks the dot size represents peak amplitude and color its prominence (ratio of peak amplitude to that of the next peak). The same representation is used in Figures 5C,D. Based on the magnitude of the periodic fraction one can divide the fluctuations into a periodic and non-periodic dominated regime as indicated in Figures 5C,D.
While Figure 5C suggests that f w is a control parameter for the magnitude of fluctuations in general, i.e., including periodic and non-periodic ones and also the magnitude of the periodic oscillations (PSD) initially scales with f w , the frequency of the dominant peaks do not scale with f w as shown in Figure 5D.
As already established in section Differentiation of Capillary Fluctuations From Gaussian Noise, the significant fluctuations FIGURE 3 | Pressure drop dp and saturation S w fluctuation (S w fluctuations were magnified by factor 5 for better visibility) from Figure 2 for fractional flows 1.0 ≥ f w ≥ 0.75 showing an almost perfectly Gaussian behavior (dotted lines represent a Gaussian fit to the respective data). For dp the correlation coefficient R 2 is 0.94-0.96. For S w the correlation coefficient ranges between 0.41 and 0.69, which is likely because of the lower signal-to-noise of the S w data.  Figure 2E, and rolling-averages over a window of 2, 4, and 16 data points (A). While the histogram for the original data (blue) set is bi-modal and clearly distinct from a Gaussian (dotted line), the histograms of the averages become increasingly symmetric and can be fitted very well with Gaussians (B). The Fourier spectrum (C) shows that smaller frequencies remain unaffected up to the frequency associated with the timer period of the averaging window.
occurring at intermediate fractional flows are not caused by instrumental noise. We can also largely rule out other artifacts from the experimental setup. These would manifest in a scaling of the periodicity with f o , which is not observed. Also, common effects known from traditional experimental setups for fractional flow, can be ruled out due to the design chosen in the outlined work.
In traditional set-ups one likely cause of fluctuations is the back-pressure control system. Membrane-based back-pressure regulators may oscillate when two immiscible fluid phases pass through simultaneously, or pumps with a servo control system for back-pressure control could potentially cause oscillations as well. However, in this particular setup the fluids were internally circulated, and no back-pressure control system was used, see Figure 1A. The two Quizzix feed pumps do use a servo control system for maintaining a constant flow rate. But that servo control system is only reacting on feedback of the pump piston movement, which is independent of the flow dynamics within the rock sample. Furthermore, if this set-up would induce fluctuations, then these would be visible at all fractional flows including f w = 1.0 and f w = 0.0. But, at these fractional flows no oscillations were observed.
The frequency of the dominant Fourier peaks does not show any obvious systematic, e.g., linear dependency on oil fractional flow f o = 1 − f w . If-hypothetically speaking-the periodic oscillations were an instrumental artifact one would expect a dependency with the rate of oil injection, i.e., q o = f o · q tot at constant q tot . Doubling f o would then double the oscillation frequency. That was not observed. The dark-gray dotted line in Figure 5D representing such scaling with f o = 1 − f w made to pass (by choice of a suitable offset) the dominant peak in the Fourier spectrum at f w = 0.1 (Figure 2D) does not pass the dominant peaks at f w = 0.05 and f w = 0.02 ( Figure 2E). A more general linear model (light gray dotted line) does not describe the data either. It may pass through some of the dominant peaks but miss other dominant peaks of the Fourier spectrum. That clearly rules out any instrumental artifacts related to the regulation of oil fractional flow, or also injection instabilities into the rock, as they would all be expected to scale with oil injection rate, i.e., f o .

Analysis of Capillary Fluctuations at the Darcy-Scale
In order to gain deeper insight into the cause for the oscillations, which have also been observed in several independent pore scale experiments (Spurin et al., 2020;Menke et al., 2021), in Figure 6 we show the saturation profiles for f w = 0.1, which has the longest oscillation period and the sampling rate of saturation would be sufficient to resolve both pressure and saturation. A closeup of the pressure data from Figure 2D for the same time interval of saturation profiles is shown in Figure 6A. While the pressure data goes through 5 oscillations, the saturation shows only 2 full oscillations. The respective saturation profiles for the same time interval are shown in Figure 6B. For clarity, the saturation profiles have been vertically shifted. The wave structure is better visible in the space-time contour plot shown in Figure 6C.
The saturation profiles show a peak that travels from inlet to outlet at constant velocity. The dashed black lines in Figure 6B represent the wave, which is better visible in the space-time plot in Figure 6C. Even though Figure 6C may appear at first glance as a standing wave, a water saturation peak (red) moving showing pressure data in blue and saturation data (sampled at much lower frequency) in red. Respective saturation profiles (B,E, which have been vertically shifted for clarity) exhibit a water saturation peak, which moves from inlet to outlet at constant velocity, forming a traveling wave, which repeats one the peak reaches the outlet. The wave structure in (B,E) has been emphasized by adding the dashed black lines. The wave structure is better visible in the space-time contour plots in (C,F). For f w = 0.10 the wave velocity is v Sw = 0.75 mm/ min and for f w = 0.05 the wave velocity is between 1.7 and 2.1 mm/min.
horizontally from left to right at different times t (vertical coordinate). The fact that in the space-time plot in Figure 6C diagonal structures at constant angle and constant spacing are found suggests that the water saturation peak moves at constant velocity, which ultimately reveals a traveling wave structure, which through its repetitive pattern is emphasized in the spacetime representation. A similar behavior is observed for f w = 0.05 shown in Figures 6D-F. In the experiment the behavior was reproducible. In total, 3 drainage and 3 imbibition cycles were performed. The experiment started with the 1st drainage cycle arriving at connate water saturation of ∼4%. It was followed by the 1st imbibition cycle arriving at a residual oil saturation of ∼48%, which is typical for a strongly water-wet rock sample such as the non-aged Fontainebleau sandstone. Then the rock was cleaned in-situ by flushing with iso-propanol without removing the sample. Then the 2nd drainage and 2nd imbibition cycle were performed, followed by in-situ cleaning and finally conducting the 3rd drainage and 3rd imbibition cycle. The respective relative permeability computed from the raw data (such as the 1st drainage cycle from Figure 2) are displayed in Figure 7. For all 3 cycles the relative permeability data overlaps indicating excellent repeatability of the experiment. We observe overall larger noise amplitudes in the drainage cycles than for the imbibition cycles (not shown), which suggests that hysteresis may play a role.
By comparing the speed of the traveling wave from the experiment with the prediction of fractional flow theory we can test whether the wave propagation is a simple fractional flow phenomenon. According to fractional flow theory, v D = x D /t D = df w /dS w , which is shown as an inset in Figure 8B.
x D = x/L is the dimension distance (L is the length of the core) and t D = PV injected /PV the dimensionless time (PV =pore volume). For f w = 0.10 the wave velocity in the experiment is v = 0.76 mm/min, which corresponds to a dimensionless velocity v D = 0.46. For f w = 0.05 the wave velocity v is between 1.7 and 2.1 mm/min, which corresponds to dimensionless velocity v D = 2.1 − 2.6.
In Figure 8B, we show the fractional flow curves for the saturation range of interest, and in Figure 8C, the respective v D = df w /dS w . For f w = 0.05, the experimentally measured v D is more compatible with the imbibition curve, whereas for f w = 0.10 the wave velocity is more compatible with the drainage curve. At first glance that appears somewhat inconsistent. However, when taking a closer look at the fractional flow construction involving two fractional flow curves, which is analogous to a moving an oil bank known, e.g., from polymer flooding (Lake, 1984), we notice a somewhat subtle detail, which is an inflection point around f * w ≈ 0.12 above which the water saturation of the bank (wave peak) increases compared with the initial saturation S w,i for the respective flow. An example is shown in Figure 8A for f w = 0.15. If we consider the experimental uncertainties in terms of saturation measurements and instead of using dotted blue lines in Figures 8D,E as initial saturations but the lowest saturation observed, it is very well possible that in the experiment the inflection point is at 0.07 < f * w < 0.10. That would imply that for f w = 0.10 we have indeed a water bank pushed by oil (which is drainage) and for f w = 0.05 we would have an oil peak pushed by water (imbibition). From the experiment, it is difficult to distinguish between the two FIGURE 7 | Relative permeability for water (blue) and oil (red) for the SCAL experiment with the strongly water-wet Fontainebleau sandstone sample for the 1st, 2nd, and 3rd cycles for drainage (A-C) and imbibition (D-F). Each cycle consisted of drainage, followed by imbibition, and sub-sequent in-situ cleaning of the sample with iso-propanol inside the setup without unmounting/remounting the sample. There is excellent repeatability between individual experiments. The dotted lines show a Corey fit of the data. The fractional flow functions f w (S w ) for the first drainage and imbibition cycle are shown in (C,F), respectively. situations based on saturation profiles because the phenomenon is repetitive, and it is not fully clear, which phase is displacing and which phase is displaced. Within the experimental uncertainties, it is very well possible that the f w = 0.10 experiment is above f * w , whereas the f w = 0.05 experiment is below the inflection point f * w . That would explain why the experimentally observed wave velocities follow the prediction from fractional flow theory v D = df w /dS w for drainage at f w = 0.10 and imbibition for the f w = 0.05 shown in Figure 8C.
The correct modeling of the traveling water peak involves hysteresis in the relative permeability-saturation functions. For the situation encountered in the experiment, the relative permeability would actually reflect a scanning curve, which needs to be considered in respective analysis (Plohr et al., 2001). The wave for each fractional flow consists of a water peak moving from inlet to outlet, which has water displacing oil at the front and oil displacing water at the back. The respective fractional flow constructions are shown in (B), where the flow line is constructed from the intersection of the drainage fractional flow curve at respective f w step. The fractional flow solution is the intersection between the blue line with the imbibition fractional flow curve (Lake, 1984). Note that depending on the exact value of the initial saturations (D,E) there is an inflection point 0.07 < f * w < 0.12 above which the water saturation from the fractional flow solution increases compared to the initial saturation (A), while below f * w it decreases (B). That would imply a change from oil pushing a water bank (drainage) at f w > f * w to water pushing an oil bank (imbibition) for f w < f * w . That could explain why the experimentally observed wave velocities follow the prediction from fractional flow theory v D = df w /dS w for drainage at higher fractional flows and the imbibition line for smaller fractional flows. That transition is within the experimental uncertainty between f w = 0.10 and f w = 0.05. A saturation profile with moving banks has an overall higher total mobility than a flat saturation profile for the saturation range S w,c < S w ≤ 0.4 − 0.5 (F). The two experiments at f w = 0.10 and f w = 0.05 are added as points, where for the computation of mobility, respective bank saturations and lengths have been taken from (D) and (E), respectively.
Such traveling waves with non-monotonic saturation profiles are known to form under specific conditions for instance for infiltration of water into dry soil (Schneider et al., 2018;Kmec et al., 2019;Mitra et al., 2020) and can be theoretically described with a range of models ranging from generalized Darcy-scale formalisms with dynamic capillary pressure (Schneider et al., 2018) to hybrid models combining Darcy physics and pore scale aspects (Kmec et al., 2019). Mitra et al. (2020) provides a mathematical model, where for specific choices of parameters, the solution forms spirals in the S w − p c space. The hysteresis model by Plohr et al. (2001) demonstrates that for specific choices of parameter values "loop", i.e., recurring solutions can be found, which implies stable orbits in phase space (Strogatz, 1994). These consist of a sequence of water and oil bank, where the respective final saturation and fractional flow serves as an initial condition for the next wave sequence. In that way, the sequence becomes recurring. Although Corli and Fan (2018) raise questions around non-uniqueness, the principle does provide a reasonable working hypothesis about the origin of the periodic waves observed in the experiments.
Whether periodic wave structures are observed depends on the initial condition S w,i for each fractional flow f w , the fractional flow f w and the underlying relative permeability functions, and potentially also the domain size L. It cannot be excluded that the large but non-periodic oscillations in Figure 2C are in the end related to the same fractional flow effect, but cannot form a closed loop (Plohr et al., 2001). More generally, phenomena at the Darcy scale are only visible if a permissible fractional flow solution exists, otherwise they would be dampened and a periodic trigger from e.g., initial or boundary conditions would decay.

Physical Interpretation of the Wave Characteristics of the Observed Fluctuations
While mathematically the solutions are permissible by fractional flow theory, the question is now why from a physical perspective a bank solution would form. A solution with moving banks would be plausible if, for the same total pressure drop, more total flux was transported than for a flat saturation profile. Whether that is the case can be assessed by computing total mobility for a flat saturation profile and compare it with the respective total mobility for the situation with moving banks. For the flat saturation profile the total mobility of water and oil phases, weighted with the respective fractional flow using drainage relative permeability (because overall it is a drainage experiment). Note that the definition in Equation (3) is different from standard fractional flow theory and in that sense represents a pseudo-mobility. The weighting with the respective fractional flow is done in order to estimate the total flux based on the contributions of drainage and imbibition waves from segregated flow. For the situation with a moving bank, we assume based on the experimental observation that the saturation profile is composed of one water and one oil bank and compute the total mobility each bank separately. In absence of a rigorous fractional flow solution, which would first require a very detailed treatment for hysteresis, we simply estimate the top and bottom saturations by adding a typical saturation fluctuation amplitude of δS w = ±0.1 or δS w = ±0.2. For the water bank we use imbibition relative permeability and for the oil bank drainage relative permeability. We then weight the contribution of water and oil bank by their respective length in the saturation profile, i.e., the water bank is about 10-20% in terms length of the total saturation profile l w = 0.1 − 0.2 (see Figures 6C,F), i.e., λ banks = l w λ water bank + (1 − l w )λ oil bank . In order to assess whether the bank solution has higher mobility we divide λ banks /λ flat , which is plotted in Figure 8F. We see that depending on the magnitude of δS w the bank solution has a higher mobility for saturations below S w < 0.4 − 0.5, corresponding to approximately f w < 0.4, which is also the range in which the large fluctuations in pressure and saturation are observed (see Figure 5C). Note that the fractional flow picture does not capture the effects of capillarity, which can improve the microscopic sweep (Khorshidian et al., 2021). The remaining question is what actually triggers the fluctuations in general and waves in particular. In the following sections, we will establish from pore scale fractional flow experiments conducted in a similar steady-state manner as the core flooding experiments in Figures 2-5 that periodic variations of pressure and saturation already occur on the oil-cluster scale and ultimately caused by pore scale displacement events (Berg et al., , 2014Spurin et al., 2020;Menke et al., 2021). The pressure and connectivity signature and 3D imaging of respective pore scale fluid distributions observed in the fractional flow micro-CT experiment suggest that pore scale displacement events are the underlying cause of fluctuations that are still visible in experiments commonly assumed at Darcy scale. They also likely provide the trigger for the waves from Figures 6C,F.

Dependency of the Fluctuations on Wettability
In order to investigate the dependency on rock wettability, a systematic study was conducted using twin samples of Bentheimer outcrop rock taken from the same block and selected for homogeneity and similar properties in terms of porosity and permeability. The rock was then de-saturated and aged with 3 different crude oils ("crude 1, " "crude 2, " "crude 3") with properties listed in Tables 1, 2. In Figure 9 the raw data of the respective steady-state imbibition experiments (analogous to Figure 2A) are shown.
The pressure fluctuations are largest for the sample aged with crude 1 and become much smaller for Bentheimer samples aged with crude 2 and 3 ( Figure 9A). Based on the trends of residual oil saturation S o,r and water endpoint relative permeability k r,w (S o,r ) shown in Figure 9B the sample aged with crude 1 is the most water-wet case, whereas crude 2 and 3 become increasingly intermediate-wet. This confirms the experience in the SCAL community that significant fluctuations are mainly observed for water-wet cases.
When normalizing the standard deviations of the pressure fluctuations by the mean pressure along the lines of the learnings from Figure 5 we find that over-proportionally-large pressure FIGURE 9 | Dependency of the fluctuations on the wetting state of the rock. Pressure fluctuations are largest for the Bentheimer rock sample aged with crude 1 and become much smaller for Bentheimer samples aged with crude 2 and 3 (A). The trends for residual oil saturation S o,r and water endpoint relative permeability k r,w (S o,r ) suggest that the wettability for the sample aged with crude 1 is more water-wet, whereas samples aged with crude 2 and crude 3 become increasingly intermediate-wet (B). The standard deviation of the pressure fluctuations normalized by the mean pressure for each fractional flow f w is largest for crude 1 (which is the most water-wet case) and becomes significantly less for crude 2 and 3 (C).
fluctuations are only observed for the crude 1 case, i.e., the most water-wet sample. However, the trend toward crude 2 and 3 cases and respective wetting states seems highly non-linear and also not entirely consistent. For instance, in the crude 3 case we do observe fluctuations in electrical conductivity, which suggests that in intermediate cases fluctuations at the pore scale could still be present but due to the much smaller capillary pressure (Lin et al., 2019a) show a much smaller signature in the pressure but still are subject to ganglion dynamics. This is in line with porescale studies, in which fluid redistributions were still happening in mixed-wet scenarios, though much slower pace compared to water-wet situations . At intermediate wetting states, the connectivity of the brine phase, which controls the electrical conductivity (Liu Z. et al., 2018), is perhaps also more sensitive to ganglion dynamics, because water-films may be overall less connected. We begin with a number of model calculations for typical pore body and throat sizes in sandstone rock of d = 25 µm diameter for pore throats and d = 60 µm for pore bodies . Assuming furthermore an interfacial tension of σ = 25 mN/m for a typical crude oil-brine system (water-decane with highest purity n-decane would be ≈ 50 mN/m but polar components in crude oil can reduce interfacial tension below 1 mN/m, therefore 25 mN/m is a good average value) we can estimate capillary pressure

Energy Scale of Fluctuations
(with contact angle θ ≈ 0 • for a water-wet system and radius r = d/2) for a water-wet system as p c = 0.017 bar for pore bodies and p c = 0.04 bar for pore throats, which are also typical values for the pressure fluctuations in the steady-state experiment in Figures 2-5 (the maximum fluctuations can be a factor 3-4 larger but for most f w the range is very similar) (Berg C. F. et al., 2020). This suggests that the pressure fluctuations are likely caused by events on the capillary energy scale, which are typically pore scale events. Therefore, pore scale displacement events related to capillary instabilities (Singh et al., 2017Primkulov et al., 2019) are the most likely cause for the significant fluctuations.
The associated pressure-volume work of displacing the wetting-phase by non-wetting phase fluid (Morrow, 1970).
In a pore body ( V ≈ 1.1 · 10 −13 m 3 for a spherical pore) is 4.5 · 10 −10 J and for a pore throat ( V ≈ 1.2 · 10 −14 m 3 ) 4.9 · 10 −11 J. This can be expressed in terms of thermal energy units k B T, where k B = 1.38 · 10 −23 J/K is the Boltzmann constant and T = 300K (≈ 26.7 • C) the temperature arriving at 1.1 · 10 11 k B T and 1.2 · 10 10 k B T for pore bodies and throats, respectively, which is also summarized in Table 3. Since the associated energy of the pressure fluctuations is many orders of magnitude larger than k B T these fluctuations can be classified as non-thermal (Winkler et al., 2020).
To put this into perspective, we can also estimate the kinetic energy that is associated with laminar single-phase flow at a typical injection rate 1 foot/day , which is typical for a water flood in oil reservoirs far away from the well bore, which translates for a porosity φ = 0.25 into a linear flow velocity of v = 1.4·10 −5 ms −1 (Berg and van Wunnik, 2017).
The kinetic energy is for an aqueous phase fluid with density ρ = 1000 kg m −3 around 1.1 · 10 −20 J ≈ 2.63 k B T, which is roughly on the thermal energy scale. The linear flow velocity during a Haines jump, which can be computed numerically with a phase-field type of method and validated with experimental data, is around v = 1.0 · 10 −2 ms −1 (Armstrong et al., 2015). That translates into a kinetic energy of 5.7 · 10 −15 J ≈ 1.4 · 10 6 k B T, which is again several orders of magnitude larger than the thermal energy scale.
Note that for pressure-volume work and kinetic energy of pore-scale displacement events there is always a range, which depends on the non-local nature of the Haines jumps, pore size, wetting conditions etc. The important point to make here is that the kinetic energy is orders of magnitude larger than the thermal energy scale, which is consistent with the classification as non-thermal fluctuations. Since the energy scale of these fluctuations is many orders of magnitude above the thermal energy scale, increasing temperature cannot trigger pore scale displacement events (unless phase changes or disintegration of material occur).
We observe pressure fluctuation amplitudes in the range of capillary pressure representative, typical pore throats of sandstone rock and indicative of pore scale displacement events as origin (Figures 2-5). However, the magnitude of the saturation change in relation to the volume of an individual pore suggests that many individual pores are involved in the phenomenon. The pore volume of the Fontainebleau rock sample is V pore,tot = 2.85 cm 3 , whereas the pore volume of a single pore of d = 60 µm diameter is V pore,single = 1.1 · 10 −7 cm 3 , which means that a saturation change of 5% needs to involve thousands of individual pores. The associated energy scale can be estimated as p c · V pore,tot · δdS w with p c ≈ δdP w = 0.15 bar resulting in 0.043 J or 10 19 k B T, which is significantly larger than the energy scale of displacements involving individual pores.
That means that multiple individual pores are involved in these fluctuations. From a pore-scale perspective this is not too surprising as experiments have already shown that displacement events are not necessarily always isolated but can be cooperative causing avalanche-like cascading sequences of individual displacements adding up to an overall significant event, that involves many pores and lasts over a time period of between 0.5 and 2 s (Armstrong et al., 2014a). What is new is that the fluctuations measured at the Darcy-scale are periodic, which means that the respective pore-scale displacement events would also need to be repetitive under fractional flow conditions. In order to address that question, in the following sections time sequences of 3D imaging are presented, which give further insight into the pore-scale origin of the fluctuations.

Time Scale of Fluctuations
When comparing experiments of different scales, both, the magnitude of time and length differences need to be considered. The time series and frequency analysis of the pressure fluctuations shown in Figures 2, 4C, 5D suggest a frequency range of fluctuations between 10 −3 − 10 −2 Hz corresponding to a time period of T = 1/f = 1 − 20 minutes, which also reflects the dominant peaks in Figure 5D. The time scale of the saturation fluctuations is between 5 and 60 min, which is also about the time scale of the traveling waves from Figure 6.
In Figure 10, we put the time scale of pressure and saturation fluctuations and traveling waves into perspective with characteristic time scales encountered in multiphase flow in porous media. Displacement events in individual pores such as Haines jumps and snap-off are on the time scale around 1 ms (Armstrong et al., 2015). Cascading displacement events involving multiple pores show a relaxation time scale around one second (Armstrong et al., 2014a). The time scale of pressure fluctuations observed in this work, in particular the periodic oscillations in Figures 2-4, are on the order of 1-20 min, which is also the time scale of periodic pressure fluctuations reported in recent studies (Spurin et al., 2020;Menke et al., 2021) see also the supporting information of (Spurin et al., 2020). The time scale of the saturation fluctuations is 30-60 min, which is similar as the characteristic time of traveling waves from Figure 6, which is still faster than the slow relaxation behavior observed in (Schlüter et al., 2017). Figure 10 suggests that in terms of time scale the pressure and saturation fluctuations observed in the SCAL experiments are a phenomenon ranging right between pore scale and Darcy scale flow physics (Khorshidian et al., 2021). Imaging at the Darcy length scale may not give further insights into the underlying mechanisms and pore scale imaging is required. Based on Figure 10, we can expect a time scale of the order of minutes, which means that imaging at the time scale of few seconds will be sufficient to resolve the phenomenon. 3 | Pressure-volume work (p c V) for model scenarios of immiscible displacement and kinetic energy for 1 ft/day flow in units of Joule (J) and thermal energy (k B T).

Diameter (µm) Volume (m 3 ) Pressure volume work (J) Kinetic E (J) Energy (k B T)
p c V pore throat 25 1.2 · 10 −14 m 3 4.9 · 10 −11 1.2 · 10 10 p c V pore body 60 1.1 · 10 −13 m 3 4.5 · 10 −10 1.1 · 10 11 Kinetic energy laminar flow 60 1.1 · 10 −13 m 3 1.1 · 10 −20 2.63 Kinetic energy Haines jump 60 1.1 · 10 −13 m 3 5.7 · 10 −15 1.4 · 10 6 FIGURE 10 | Time scale of the fluctuations observed in this work (red) put into perspective with other characteristic times scales encountered in porous media. Displacement events in individual pores such as Haines jumps and snap-off are on the time scale around 1 ms (Armstrong et al., 2015). Cascading displacement events involving multiple pores show a relaxation time scale around one second (Armstrong et al., 2014a). The time scale of pressure fluctuations observed in this work and related studies (Spurin et al., 2020;Menke et al., 2021) are on the order of 1-20 min and the saturation fluctuations and traveling saturation waves from Figure 6 are on the order of 30-60 min. That is still faster than the slow relaxation behavior observed in Schlüter et al. (2017). Darcy flow, e.g., at typical field rates of 1 ft/day over a characteristic length scale such as a capillary dispersion zone of 3-5 centimeters length occurs at 7-10 h.

Pressure and Saturation Fluctuations in Micro-CT Fractional Flow Experiments
In the steady-state micro-CT flow experiment conducted at a much smaller sample (4 mm diameter, 2 cm length) also fluctuations in pressure and saturation are observed as shown in the example in Figure 11. These events are re-occurring between periods of meta-stable and predominantly connected pathway flow (Rücker et al., 2015b). Associated changes in pressure in Figure 11A and saturation in Figure 11B have similar magnitude as the fluctuation amplitudes in the SCAL steadystate experiment conducted at a much larger sample shown in Figures 2-5. The pressure fluctuation event in Figure 11 of 10-20 mbar amplitude is similar to the fluctuation amplitude in Figure 2C. The monitored fluctuations at the pore-scale are significantly larger than instrumental noise and are not caused by instrumental artifacts, which is evident from the larger scale analysis, the chosen experimental set-up as well as the additional information revealed by the imaging of the interior of the rock during flow (Figure 11). Note that-as discussed in section Energy Scale of Fluctuations-pressure fluctuations are mainly caused by capillary effects. Therefore, even though the micro-CT sample has a significantly larger permeability, but pore throat size distribution is actually similar to that of sandstone rock of comparable porosity and permeability as in the SCAL experiments (see Figure 1 in Berg et al., 2014), the magnitude of pressure fluctuations is similar.
These observed fluctuations and associated pore-scale events are related to change in fluid connectivity as the Euler characteristic (Herring et al., 2013), a mathematical measure for connectivity displayed in Figure 11G, and the connected pathway relative permeability computed from the fluid phases imaged by micro-CT displayed in Figure 11D shows. Even though the sample is significantly smaller than in the SCAL experiment, the associated saturation change is much bigger than the volume of individual pores, which implies that these events involve multiple individual pores. The event statistics, which was obtained by comparing 3D images of water and oil phases of subsequent time steps (Rücker et al., 2015b) displayed in Figure 11C show that the events observed in pressure and saturation consist of hundreds of individual pore level displacements.
That is consistent with cooperative displacements of advancing fluids fronts (Jung et al., 2016;Singh et al., 2017;Primkulov et al., 2019) and the observation in Armstrong et al. (2014a) with avalanche-like cascades of individual pore filling events operating overall at the cluster scale. The repetitive sequence of oil and water filling events are indicative of waves passing through the field of view, although the domain size in micro-CT experiments does not permit to resolve the spatial extent of the wave but rather see repetitive behavior in the average S w . The event statistics in Figure 11C shows that the macroscopically-visible events start with sequences of water filling events, which are then followed by sequences of capillary-driven oil filling events, which may provide the trigger for the Darcy scale observations, where traveling waves of water saturation peaks were observed ( Figure 6B) and permissible by fractional flow solutions (Figure 8), which operate in the viscous limit. This draws a direct connection between Darcy-scale observable flow regimes (which have to be permissible by fractional flow solution, otherwise they would be damped and decay) and cooperative pore-scale FIGURE 11 | Pressure drop (A), water saturation (B), event statistics (C) connected pathway relative permeability k r of connected water and oil phases (D), capillary pressure p c obtained from the mean curvature (E), interfacial area, (F) and the Euler characteristic for the oil phase (G) for the experiment at f w = 0.5 and v = 4 · 10 −4 m/s. The fluctuation events with a clear pressure and saturation signature are caused by a sequence of water and oil filling events. They are causing changes in the connectivity of water and oil phases as visible in the connected pathway relative permeability, which was computed by Stokes flow simulations of the imaged water and oil phases in relation to the (analogous computed) permeability of the respective pore space. The capillary pressure largely follows the saturation and connectivity behavior, which all typically lag behind water or oil filing events, respectively. Note that capillary pressure shows a clear signature also for the 2nd filling event, which is much less visible in saturation, interfacial area (Joekar-Niasar and Hassanizadeh, 2011) and Euler characteristic (Herring et al., 2013). displacement behavior. It furthermore suggests that analogous to ganglion dynamics (Rücker et al., 2015b) the dynamics is initiated by snap-off which is a capillary instability, where capillary energy is converted into kinetic energy, which can then drive event cascades. The pore-scale fluid configuration eventually reverts back to the original situation in terms of saturation, pressure drop and connectivity of each phase as suggested by Figures 11A,B,D. There is a range of possible reasons from the bubble instability (Unsal et al., 2007a,b;Unsal et al., 2009), which can in principle also occur in large pores inside of rock, relaxation time scales (Armstrong et al., 2014a) being longer than the re-supply of injected liquid at fractional flow f w , questions about the energy efficiency of connected pathways (Spurin et al., 2021) and inertial effects during displacement events causing ganglion dynamics (Rücker et al., 2015b), which can all explain why established flow paths can become unstable and disconnection processes are followed by re-connection processes. The capillary pressure determined from the liquid-liquid curvature (see, e.g., Lin et al., 2019b) shown in Figure 11D largely follows the saturation and connectivity behavior, which typically lag behind water or oil filing events, respectively. The magnitude of the capillary pressure changes is about 10 mbar, which is significant and suggests that capillarity is involved in the process. The system can reduce its capillary energy by water filling events, but because of the flux and fractional flow boundary conditions, oil-filling events temporarily increase the capillary energy again. This suggests that the condition is not only viscous-unstable as the fractional flow assessment suggests, but also capillary unstable, based on the competition between minimizing potential energy (where capillary energy is a significant fraction) and fractional flow boundary conditions. Note that pore filling events in Figure 11C are computed from difference images, which is the reason why they may appear somewhat earlier than the capillary pressure response in Figure 11E, which are computed from only one image.
FIGURE 12 | Breakup of oil clusters (A) and coalescence (B) processes at the pore scale. Breakup occurs often by simultaneous snap-off processes in a cluster thereby preserving topological loop structures (see insets in A). Clusters can also grow by coalescence, which is overall more a frontal propagation, whereby connectivity is maintained, and formation of loops is largely avoided. Both breakup and coalescence events impact cluster volume, i.e., saturation and cluster length, i.e., the macroscopic capillary number, which is proportional to the viscous pressure drop over the cluster length.
The question is now how we can use 3D pore scale imaging to understand more about the origin of these repeating events with cascading pore-level event structure.

Pore Scale: Breakup and Coalescence Processes-Ganglion Dynamics
In order to characterize the processes at the pore scale, we first establish a general insight into the evolution of oil clusters under fractional flow conditions, which is more complex than the ganglion dynamics during imbibition alone (Rücker et al., 2015b), because both wetting and non-wetting phases are permanently co-injected. One of the questions is how that coinjection at fractional flow f w is translating into the flow at pore scale, i.e., what the associated flow regimes are.
At the pore scale the fractional flow can be divided into time periods of connected pathway flow and time intervals with topological changes in the fluid configuration. These can be decomposed into a more imbibition-dominated regime with breakup of oil clusters by sequences of snap-off and piston-like displacement processes as shown in Figure 12A, and a more drainage-dominated regime with oil cluster growth by Hainesjumps and coalescence as shown in Figure 12B. Note that as clusters breakup /shrink or grow, also their length changes, which FIGURE 13 | Fractional flows f w (0.8, 0.5, and 0.2) and cluster-based capillary number Ca cluster from Equation (7) for the micro-CT flow experiments conducted at flow rates of 3, 30, and 300 µl/min (A). For all flow rates and fractional flows Ca cluster ≪ 1, which means that the flow regime is capillary dominated. The bounding line for Ca macro = 1, which is the transition to the viscous-dominated regime is plotted as red line in (B) in the flow velocity (left axis) and injection rate (right axis) vs. cluster length (C) space. Also, the largest cluster falls below this line even at the highest flow rate, which means that the flow regime is capillary dominated. Figure 11 show long time periods without any change but also time periods with discrete events. Based on the associated dynamics at the pore scale, respective pore scale flow regimes can be categorized into predominantly connected pathway (A) and ganglion dynamics (B). A "phase diagram" can be developed in analogy to Figure 12 but now for the whole field of view instead of individual clusters only, where connected pathway flow is characterized by a fixed point (C) and ganglion dynamics as trajectories (D) in the saturation-macroscopic capillary number space.
where L cluster is the length of the respective oil cluster, µ w the wetting phase viscosity, v w the wetting phase Darcy velocity, σ ow the wetting/non-wetting phase interfacial tension and θ the contact angle, and p c the capillary pressure at breakthrough.
In that way, for each oil cluster an individual capillary number can be assigned, which is proportional to the cluster length L cluster and increases or decrease according to how the cluster grows by coalescence or shrinks by breakup. Figure 12 is based on actual data from the micro-CT flow experiment conduced in the sintered glass sample. Conceptually similar behavior is expected also in other porous media but may change in terms of quantity, because of e.g., different scaling of permeability over capillary pressure.
Note that these examples characterize local flow regimes. These local regimes can also dominate the macroscopic flow regime in a macroscopic drainage or imbibition experiment, but they can also co-exist, or alternate in time, in fractional flow.

Identification of Pore Scale Flow Regimes by Collective Cluster Dynamics and the Pressure/Saturation Signature of Ganglion Dynamics
The micro-CT flow experiments involved three different fractional flows f w = 0.8, 0.5, 0.2 and three different injection rates of q tot = 3, 30, 300 µl/ min. Using the largest from the imaged cluster size distribution from Rücker et al. (2015a) these flow rates are then translated into the maximum cluster-based capillary number (Equation 7). The overview of all experiments is shown in Figure 13A in the fractional flow-macroscopic capillary number space.
All measurements obtained fall significantly below the Ca macro = 1 bounding line, which is the transition to the viscousdominated regime as shown in Figure 13B. This means that in every respect the flow regime is capillary dominated below the onset for viscous mobilization (Oughanem et al., 2015).
Nevertheless, at specific saturations and higher flow rates we observe ganglion dynamic behavior as shown in Figure 11. In Figure 14 we show the respective pore scale fluid distributions. In Figure 14A shows an example for a connected pathway flow period without any significant fluctuations. In Figure 14B shows the pore scale fluid distributions for the fluctuating event from Figure 11. Figures 14C,D are an extension of the concepts from Figure 12 but now for the whole oil saturation in the imaged field of view and not only individual clusters, the capillary number is calculated using the a volume weighted average cluster length. The connected pathway flow regime is characterized by a single point in Figure 14C, whereas the fluctuating regime is characterized by a trajectory in Figure 14D. Based on that observation, the saturation-macroscopic capillary number space displayed in Figures 14C,D can be understood as "phase diagram", where connected pathway flow is characterized by a fixed point and ganglion dynamics as trajectories.
The ganglion dynamic events from Figure 11 are characterized by a diagonal movement between two discrete states in the saturation-capillary number space as shown in Figure 14C. This repeated switch between two states is very similar as the pressure and saturation fluctuations between an maximum and a minimum saturation level, and a maximum and minimum pressure drop, observed in Darcy-scale SCAL experiments in Figures 2-5. From the change in saturation but also the associated 3D visualization of the respective pore scale fluid distribution shown in Figure 14B it is clear that such ganglion dynamics events involve cooperative displacement in a large number of individual pores.

SUMMARY AND CONCLUSION
In this work we establish that significant fluctuations in pressure and saturation well visible in multiphase flow experiments on centimeter-sized samples are not instrumental noise (or artifacts) but caused by cooperative pore scale displacement events. There is consistent evidence ranging from special core analysis (SCAL) core floods on samples few centimeter length and diameter and micro-CT flow experiments on few mm sized samples. Fluctuations are observed in several independent SCAL studies, show non-Gaussian statistics. Displacement events imaged by beamline based micro-CT flow experiments have similar characteristics in terms of amplitude and saturation change as the SCAL experiments.
Until now, the common view has been that porous media, e.g., rock samples of few centimeter in size are above the REV scale and therefore, pore-scale events should have averaged out. It is, however, not too surprising that pore-scale events do not average out at sample sizes of few centimeters since non-wetting phase cluster sizes can reach into the several mm to few centimeter range (Georgiadis et al., 2013;Armstrong et al., 2014b) and that the pressure signature of respective ganglion dynamics (Rücker et al., 2015b) remains visible. In order to average out, if at all, sample sizes that could enclose several free non-wetting phase clusters would be required.
By focusing on the specific case of steady-state experiments where wetting and non-wetting phases are co-injected at varying fractional flow, we can exclude transient effects. Even though in many cases of such fractional flow experiments, after injection of hundreds of pore volumes, pressure and saturation eventually stabilize, there are many other reported cases, where significant fluctuations in pressure and saturation persist in a wide range of different rock types, experimental conditions and scales investigated, which suggests that the phenomenon is somewhat universal and applicable to a wide range of rock types (Datta et al., 2014a,b;Masalmeh et al., 2014;Reynolds and Krevor, 2015;Sorop et al., 2015;Gao et al., 2017Gao et al., , 2019Gao et al., , 2020Lin Q. et al., 2018;Alcorn et al., 2019;Clennell et al., 2019;Lin et al., 2019a;Spurin et al., 2019;Wang and Masalmeh, 2019). Based on capillary number and pore-scale mechanisms involved the cause of the fluctuations is of capillary nature. The magnitude of the fluctuations depends strongly on the wetting state and may be noticed only in more water-wet cases, which does not necessarily mean that ganglion dynamics is absent in mixed-wet cases . It is still uncertain whether bi-continuous oil-water configurations and associated connected pathway flow regimes can be established for specific pore space morphologies and the spatial wettability distributions (Lin et al., 2019a).
In the analyzed experiments, the amplitude of the pressure fluctuations is at a capillary energy scale, which suggests that pore-scale events are the most likely cause. The volume associated with the saturation fluctuations are many orders of magnitude larger than individual pores suggesting that the pore-scale displacement events are cooperative and involve many pores, which is also reflected in the hundreds of individual pore scale displacements observed in the changes of the pore scale fluid distributions imaged in the micro-CT flow experiments.
The root cause for fluctuations are discontinuities associated with pore scale displacement events. However, the involvement of hundreds to thousands of individual pores stretches the limits of what is explainable at pore scale just by cooperative displacements of advancing fluids fronts (Jung et al., 2016;Singh et al., 2017;Primkulov et al., 2019) and avalanche-like cascades of individual pore filling events (Armstrong et al., 2014a). Cooperative displacements with locally either more imbibition and drainage dominated behavior with associated breakup and coalescence processes is ultimately reminiscent of what we would describe as ganglion dynamics (Rücker et al., 2015b) although the necessity for ganglion dynamics has not been fully comprehended. Associated flow regimes can be identified in the capillary number-saturation space where during breakup and coalescence, oil clusters follow specific trajectories. In the transition to the Darcy scale, e.g., when averaged over all oil clusters to characterize the average macroscopic behavior, connected pathway flow without significant saturation changes is a fixed point and ganglion dynamics form trajectories. The observation of fluctuating flow where saturation, pressure drop, and fluid connectivity change repetitively between two different states is somewhat reminiscent of ganglion dynamics, which operates on a capillary energy scale many orders of magnitude above thermal fluctuations suggesting that these fluctuations are non-thermal.
Ganglion dynamics events can indeed involve many individual pores. However, a more statistical distribution would be expected, where filling events with hundreds to thousands of individual pores would be relatively rare, because the event size distribution follows a power law (Berg et al., 2014). Therefore, from the pore scale perspective alone it is not clear why these large events would become the dominant feature, in a sometimes perfectly repetitive oscillatory manner.
The Darcy scale experiments provide the missing piece of the puzzle. At the centimeter scale oscillatory regimes with traveling waves are observed. The characteristic feature is a saturation peak traveling at constant velocity from inlet to outlet. Once the outlet is reached, a saturation peak emerges at the inlet repeating the cycle. Such waves are compatible with fractional flow principles (based on Darcy scale physics) and can be explained by a "loop" fractional flow solution (Plohr et al., 2001). Darcy scale phenomena are only observable over larger length and time scales due to permissibility by fractional flow solutions (otherwise they would be damped and decay exponentially), meaning that they enter a viscous-dominated flow regime. The trigger mechanism for such waves, however, are provided by capillary-driven ganglion dynamics events observed on the cluster scale. In a respective micro-CT flow experiments imaged by fast X-ray computed micro-tomography such recurring events were observed, which are initiated by an increase in waterfilling events. That is reminiscent to the ganglion dynamics observed in imbibition at low capillary numbers (Rücker et al., 2015b), where ganglion movement is only initiated after a sequence of snap-offs disconnecting an oil cluster, which then becomes mobile through a sequence of disconnection and coalescence processes, at a capillary energy scale, which upon snap-off is converted into kinetic energy. This may provide the trigger for Darcy-scale "loop" solutions, which consist of a moving water bank followed by an oil bank in a periodic "loop" between bounding drainage and imbibition fractional flow curves, potentially forming closed orbits in the phase space (Strogatz, 1994).
A possible physical explanation for the occurrence of fluctuations are indeed the underlying traveling waves which can transport more flux than a homogeneous saturation profile. In other words, the Darcy scale solution becomes unstable against fluctuations, because in this way, at the same pressure drop, more flux is transported, i.e., energy dissipation is minimized.
The phenomenon of the clearly visible fluctuations in core flooding experiments at the centimeter length scale involve both capillary driven dynamics at the cluster scale with pore scale origin, and aspects of fractional flow solutions in the viscous limit. In other words, these phenomena are bridging between pore and Darcy scale physics, and the complexity arises from the fact that between the two scales the physical description concept changes.
There are now several potential implications, ranging from the question on the rate of energy dissipation associated with the fluctuations (McClure et al., 2020a,b) to non-uniqueness of the "loop" fractional flow solutions of the periodic oscillations (Plohr et al., 2001;Corli and Fan, 2018) and implications for the relative permeability-saturation functions obtained from such experiments. Ultimately, we have to also ask ourselves from a practical application perspective how we can use the details of fluctuations and traveling waves revealed in this much more dedicated analysis to enrich the interpretation of the experiment, acknowledging that by averaging the information would be largely lost. For instance, fluctuations could be used as a wettability indicator. These questions are good starting points for future work on the subject.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

AUTHOR CONTRIBUTIONS
SB, MR, AG, RA, HO, FE, and MK conceptualized the research. SB, AG, HO, MR, FE, and MK planned the research. AG, MR, SB, RA, HO, NB, HvdL, and LS performed the research. MR and SB performed the data analysis and wrote the manuscript. All authors contributed to the article and approved the submitted version.