History-independence of steady-state in simultaneous two-phase flow through porous media

It is well known that the transient behavior during drainage or imbibition in multiphase flow in porous media strongly depends on the history and initial condition of the system. However, when the steady-state regime is reached and both drainage and imbibition take place at the pore level, the influence of the evolution history and initial preparation is an open question. Here, we present an extensive experimental and numerical work investigating the history dependence of simultaneous steady-state two-phase flow through porous media. Our experimental system consists of a Hele-Shaw cell filled with glass beads which we model numerically by a network of disordered pores transporting two immiscible fluids. From the measurements of global pressure evolution, histogram of saturation and cluster-size distributions, we find that when both phases are flowing through the porous medium, the steady state does not depend on the initial preparation of the system or on the way it has been reached.


I. INTRODUCTION
Understanding the physical mechanisms underlying multiphase flow in porous media is crucial to a wide variety of industrial and environmental problems, such as oil recovery, CO 2 transport and storage or ground water management. At the macroscopic scale, the flow of immiscible fluids in a porous material is described by specifying the relations between global quantities such as flow rate, pressure gradient or fluid saturation. At the pore scale, this flow is governed by the competition between capillary, viscous and gravitational forces. Understanding the link between the two levels of description requires to relate the position and shape of the interface(s) between the two phases to the values of the macroscopic variables. From an experimental point of view, two-dimensional model porous media providing direct pore-scale visualization of the flow structures are ideal tools to study multiphase flow mechanisms and their relations with global quantities in controlled, laboratory-scale situations. Over the past decades, micro channel networks etched in transparent plates [1], or prepared using molding techniques [2] and porous Hele-Shaw cells consisting of a layer of beads between two parallel plates [3] have become classical tools, to which improvements have been constantly proposed [4,5]. Experimental observations have been explained through extensive numerical simulations based on network models [6][7][8][9][10] and lattice Boltzmann methods [11][12][13][14][15][16][17], statistical models [18][19][20] and differential equations [21]. Most of the research in this area has been focused upon transient phenomena, i.e. drainage or imbibition -arising when one phase displaces the other in a porous medium. The relation between macroscopic flow parameters, fluid morphology and stability of the interface between the two phases has been thoroughly observed [3,4] and successfully modeled [18,19].
In this article, we deal with steady-state flow, where many questions are yet to be answered. In the steadystate regime, two phases are injected simultaneously into the porous medium and one observes that one or both are fragmented and transported in the form of clusters of various sizes, forming a complex flow pattern with multiple interfaces. After a characteristic time, the system reaches a steady-state in which the macroscopic flow variables fluctuate around constant values. The usual distinction between drainage and imbibition is irrelevant to describe steady-state two-phase flow, in which both processes occur simultaneously. New approaches are thus needed to understand this regime. In an effort to bring new insight, experimental and numerical studies have investigated the relations between macroscopic flow variables, and different models have been proposed to relate them to porescale flow mechanisms : Payatakes and co-workers carried out detailed experimental, numerical and theoretical studies of steady-state two-phase flow, with emphasis on the determination of relative permeabilities [1,9,10,[22][23][24][25]. The steady-state characteristics of macroscopic flow properties have also been investigated numerically by Knudsen et al. [26]. A power-law relation between pressure and steady-state flow rate has been observed experimentally by Tallakstad et al. in a 2D system [27,28], and by Rassi et al. in a 3D system [29]. Very recently, the relation between the steady-state flow rate and pressure drop has been derived analytically for two-phase flow through single capillaries [30] and through porous media [27,28,31] and also supported by extensive numerical simulation. Distributions of non-wetting or wetting clusters in the steady state have also been studied experimentally by Tallakstad et al. [27,28] and numerically by Ramstad and Hansen [32], and critical exponents were measured. It is worth noting though, that the comparison of experimental and numerical results is not always straightforward, due to the different boundary conditions used in the two cases. Yet from the point of view of sta-tistical physics, the existence of a genuine steady state is very promising to build a thermodynamic-like theoretical description of the system. In that context, it is crucial to determine whether the steady state is independent on the history of the process, or in other words, whether it is a real state in a thermodynamic sense [33]. It is well known that when drainage and imbibition occur successively, the relative permeabilities become history dependent and the pressure-saturation curves display an hysteresis [34,35], the underlying pore-scale mechanisms of which are known. Besides in the well-known magnetic and elastic systems, such history dependence and hysteresis has also been observed in different flow processes like hydrodynamic heat flow [36] and particle flow through random media [37]. However, in the case of two-phase flow through porous media, it is not trivial to predict whether such an hysteresis will come into play in the steady-state situation when drainage and imbibition occur simultaneously. It was proposed that a thermodynamic-like description for simultaneous twophase flow in porous media can be sketched [33] if the flow was history independent.
Here, we present an extensive experimental and numerical study in order to investigate the historyindependence of the steady state. Our experimental system consists of a porous Hele-shaw cell in which we simultaneously inject air and a viscous water-glycerol solution. We then compare steady-states obtained for a given flow rate with different initial conditions. We model the system by a network of disordered pores transporting two immiscible fluids. From pressure measurements and analysis of the statistical properties of the flow patterns, we observe no history-dependence for the steady state.

II. EXPERIMENTAL SETUP
We use a two-dimensional, transparent, porous Heleshaw cell. This experimental setup, shown on Figures 1 and 2, has been described in details in [27] and we recall its main features here. The porous medium consists of a random mono layer of glass beads, 1 mm in diameter, spread between the sticky sides of two sheets of contact paper. Its lateral boundaries are sealed with silicon glue. Attached on top of this layer, a Plexiglas plate with etched flow channels allows injection and evacuation of fluids into and from the porous matrix. A pressure cushion (see [27] for details) placed below the porous medium, and a thick glass plate on top prevent the system from bending when the pressure increases as fluids are injected. Clamps placed all around the setup maintain all the layers together. This way, we obtain a porous medium of constant thickness a = 1 mm, length L = 85 cm, width W = 42 cm in which the beads remain immobile. The porosity φ and absolute permeability κ 0 of the medium are found experimentally to be φ = 0.63 and κ 0 = 1.95 · 10 5 cm 2 [27].
We use the same fluid pair as Tallakstad et al. [27],  Figure 1. Sketch of the experimental setup : the 2D porous matrix consists of a disordered mono layer of glass beads spread between two sheets of contact paper. The boundaries are sealed with silicon glue. The upper part of the system consists of a Plexiglas plate with drilled inlet and outlet flow channels. A pressure cushion and a thick glass plate placed below and above the porous matrix ensure the overall rigidity of the system and maintain its thickness constant. Clamps maintain all the layers together. A lightbox illuminates the system from below and a digital camera is placed above to record images of the flow structure.
namely air as the non-wetting phase and a viscous waterglycerol solution (15-85% in mass) as the wetting phase. The latter is dyed in black with 0.1% Negrosine to obtain a good contrast on the experimental images (see Figure 3). As pointed out by Tallakstad et al., the use of air as one of the phases introduces more complexity in the two-phase flow problem for two reasons : first, it yields a high viscosity contrast with the glycerol solution, and second, its compressibility gives rise to rapid bursts or avalanches [27]. However, from an experimental point of view, it has the huge advantage of allowing us to reuse the same porous model for all experiments. Indeed, it can easily be flushed out, making it possible to obtain reproducible initial conditions. Therefore, we find it fully convenient for the present study. As shown on Figure 2, the two phases are injected simultaneously, using the same syringe pump, from 15 syringes (7 of air and 8 of water-glycerol solution), each connected to one of the 15 inlet nodes of the porous model. In the following, we note Q the total flow rate, while Q w = (8/15)Q and Q nw = (7/15)Q denote the wetting and non-wetting flow rates, respectively. To account for small temperature variations due to the heat released by the lightbox (see Figure 1), we monitor the temperature of the wetting phase at the outlet of the model. The viscosity µ w of the wetting phase is deduced accordingly using the empirical formula given in [38]. In the series of experiments presented here, the measured temperatures are in the range 24.4 − 29.3 • Celsius, giving 0.083 > µ w > 0.062 Pa.s. The viscosity of air being µ nw ≈ 1.9 × 10 −5 Pa.s, the viscosity ratio M = µ nw /µ w of the order of 10 −4 in all experiments.
With this setup, the total flow rate Q and the fractional flow F w = Q w /Q are controlled flow variables. The volumes of wetting and non-wetting fluids present in the porous matrix, V w and V nw , are free to vary with time. Thus, the saturations S w = V w /V and S nw = V nw /V , where V is the total pore volume, are free to fluctuate.  The flow is characterized by the capillary number : where µ w is the wetting phase viscosity, Q w is the total wetting fluid flow rate, γ ≈ 6.4 · 10 −2 N/m is the interfacial tension between the two phases [27] and A = W aφ is the cross-section of the porous matrix. In the present experiments, we have explored the range 3.33 10 −6 Ca 1.13 10 −4 . The highest experimental value of Ca is set by the maximum pressure that the porous model can hold. However, as we will see in Section IV, we have also explored higher values of Ca in numerical simulations, namely 1.92 × 10 −5 Ca 7.0 × 10 −2 .
Our analysis of steady-state flow and its history dependence relies on two kinds of information : measurements of the pressure inside the model, and pictures of the flow pattern. We measure the pressure P (t) in the wetting phase as a function of time t using flow-through pressure sensors (SensorTechnics 26PC0100G6G) placed at three different points of the model as indicated on Figure 2. The porous model is lit from below using a lightbox, and images of the flow structure are recorded regularly using a Nikon D200 digital reflex camera giving 2592 × 3872 pixels images with a spatial resolution of 8 × 8 pixels per mm 2 . Before further processing, images are cropped to remove boundaries, leaving us with 2400 × 3800 pixels on the final images. Figure 3(a) shows an example of steadystate image. The area of the imaged zone (see Figure 2) is large enough to contain many air clusters of various sizes. As illustrated by Figure 3(b), the high image resolution makes it possible to distinguish glass beads, air clusters, and viscous liquid. Each of these phases gives rise to a peak on the gray-scale image histogram (see Figure 3(c)). The heights of these peaks contain information about the proportion of wetting and non-wetting fluids in the system, thus about the saturations S w and S nw .

A. Principle of the experiments
We investigate the history-dependence of steady-state flow by comparing steady states obtained at the same flow rate but with different initial conditions. For this, we have performed experiments in which the flow rate is modified twice, as illustrated by Figure 4(a). The porous model is initially filled with the wetting phase only. Then, both phases are injected simultaneously at a fixed flow rate Q 1 . Once the system has reached a steady state (ss 1 ), we abruptly change the flow rate to a different value Q 2 and wait until a new steady state is established (ss 2 ). Finally, we change the flow rate back to its initial value Q 1 and let the system evolve towards a third steady state (ss 3 ). We have used different couples {Q 1 , Q 2 }, listed in Table I, to obtain different magnitudes and signs of ∆Q = Q 1 − Q 2 . We compare ss 1 and ss 3 using three criteria : the average pressure drop across the system, the fluid saturations, and the size distribution of the air clusters. From the pressures measured at the inlet, middle, and outlet of the system, we compute the pressure drops ∆P i (t) = P inlet (t)−P outlet (t) and ∆P m (t) = P middle (t)− P outlet (t) (see Figure 2).   Table  I) : (a) Temporal variation of the imposed flow rate (1/15) .Q for one syringe. -(b) Ratio of the total injected fluid volume to the available pore volume, computed from the imposed flow rate using Eq. 2 -(c) Corresponding pressure drop measurements. Plateaus characterize the different steady-states (ss). The sudden drop around t = 75h is an experimental artifact due to the refilling of the syringes, as explained in the text. Steady-states obtained before and after refilling are distinguished using the subscripts a and b. drops stabilize and fluctuate around constant average values. This behavior is identical to what has been observed previously [27] and defines steady state 1 (ss 1 ). The fluctuations of ∆P in steady state reflect the dynamics of the non-wetting clusters -transport, mergings and snap-offs -as the system explores different configurations [27,28]. Immediately after we change the flow rate to the value Q 2 > Q 1 , both pressure drops display a very rapid in-crease, and a new steady-state (ss 2 ), characterized by higher values of ∆P , is quickly established. When we change the flow rate back to Q 1 , we observe again a rapid variation of the pressure drops towards a third plateau defining steady-state 3 (ss 3 ).
Because experiments are performed at slow flow rates, and to obtain enough statistics on the measurements, they must typically run for several days. Indeed, the quality of the statistics is determined by the ratio N pv (t) of the total volume of fluids injected into the system to the total pore volume of the porous matrix, namely : where Q denotes the total flow rate, t denotes time and W , L, a and φ are the width, length, thickness and porosity of the porous matrix, respectively. Figure 4(b) shows N pv , calculated from the imposed flow rates as a function of time, for a typical experiment. In all the experiments presented here, the duration of the steady states correspond to the injection of 0.3 to 1.4 pore volumes in the system. As a consequence, it is necessary to refill the syringes one or several times in the course of an experiment. This is systematically performed using the following protocol : the outlet and inlets of the model are closed and the syringe pump is stopped immediately. The refilling process takes approximately 20 minutes. Because the vents used to close the model are not perfectly air-proof, the pressure in the system relaxes towards atmospheric pressure during this process, explaining the sudden drop of ∆P observed on Figure 4(c). However, by looking at the pictures recorded throughout the process, we have checked that this does not affect the flow structure, which remains immobile over the duration of the refilling procedure. Furthermore, we observe that once the syringe pump is restarted, ∆P retrieves its initial value after a delay originating from the compressibility of air. Thus this refilling procedure does not affect the results of the experiments. In the following, when necessary, we distinguish steady states obtained before and after refiling using the subscripts a and b, respectively (see Figure  4(c)).
B. History-independence of pressure drops Figure 5 shows the temporal evolution of the pressure drops ∆P i and ∆P m for the 6 experiments listed in Table  I. In all cases, we observe a behavior similar to what has been described in the previous section, namely rapid variations of ∆P upon changes of flow rate between plateaus characterizing steady states. Most importantly, we observe that whatever the value and sign of ∆Q, the pressure drops are similar in ss 1 and ss 3 . This suggests that the steady-state pressure drop ∆P ss only depends on the imposed flow rate, and not on the history of the system. Note also that for a given flow rate, ∆P ss values are reproducible from one experiment to the next, regardless of how steady state has been reached.  Table I.
To obtain a quantitative indication of this historyindependence, we compute average steady-state pressure drop values ∆P ss over the duration of each steady-state, thus over periods of time corresponding to the injection of 0.3 to 1.4 pore volumes in the system, as already indicated earlier. For each experiment, we compare the variation of ∆P ss between ss 1 and ss 3 to the fluctuations of ∆P ss within these two states. For this we compute the ratio : where · · · represents a temporal average over the duration of a steady-state and σ(∆P ss 1 ) ≈ σ(∆P ss 3 ) are the standard deviations of ∆P ss in ss 1 and ss 3 . The values of χ computed from the inlet and middle pressure drops for the different experiments are reported in Table I. As can be observed, almost all of these values are smaller than one. Slightly higher values are observed in Exp. 3 and Exp.4. However, our experimental temperature data suggests that temperature-induced viscosity fluctuations most likely explain this fact. Indeed, the largest temperature variations both between ss 1 and ss 3 and within ss 1 or ss 3 occur for these two experiments. There-fore, we find that within the precision of our measurements, the steady-state pressure drop values are historyindependent.
C. History-independence of saturation and non-wetting cluster size distributions Pressure measurements suggest that the observed steady-state only depends on the imposed flow rate. However, these measurements do not give us detailed spatial information. Thus, we now analyze the images of steady-state flow patterns. As explained earlier (see Figure 3), the flow structure consists of clusters of the nonwetting phase, air, surrounded by the wetting viscous water-glycerol solution. As air clusters are transported through the porous medium, they are fragmented or merged, giving many different realizations of the flow pattern. However, within steady-state, the saturation and the size distribution of air clusters remain constant on average [27]. Thus, similarly to what we have done for pressure drops, we compute average steady-state saturations and cluster size distributions and compare their properties in ss 1 and ss 3 . To obtain a good statistics on these quantities, we have chosen the frame rates so that successive images sample different realizations of the flow pattern. All averages are computed over series of images, typically 100, spanning a time range corresponding to the injection of 1/5 of the pore volume at the minimum.
As explained in Section II, the gray-scale histograms of the raw images give us a direct indication of the proportion of the two phases in the system, thus of saturation (see Figure 3). Figure 6 shows average steady state histograms for the 6 experiments listed in Table I. It is clear on all these graphs that the histograms corresponding to ss 1 and ss 3 can be distinguished from those corresponding to ss 2 : indeed, as expected, the saturation depends on the flow rate, as reflected by different air and liquid "peak" heights on the histograms. However, histograms are similar in ss 1 and ss 3 , suggesting that the steady-state saturation is also history-independent. It is important to note that histograms are directly obtained from raw images without prior processing, which excludes eventual artifacts due to image processing. However, they do not give us any information about the spatial repartition of the two phases in the system.
To obtain this information, we process the images to identify air clusters and compute their size distributions. Image processing is performed using ImageJ [39]. Raw images are thresholded to obtain binary (black and white) images on which we run a standard particle analysis algorithm 1 to identify air clusters and measure their sizes n. From steady-state images, we compute the 1. We use the "Particles4" ImageJ plugin written by G. Landini (see http ://www.dentistry.bham.ac.uk/landinig)  Table I. Note that no image processing has been applied before obtaining these histograms. Data for steady-states ss1, ss2 and ss3 are represented in black, blue and red, respectively. Different symbols refer to averages performed over different series of 100 images, namely : •, · for images in ss1a, +, × in ss 1b , , ♦ in ss2 △, ▽ in ss3a and ⊳, ⊲ in ss 3b .
normalized probability density functions of n, i.e, nonwetting cluster size distributions p(n) , where · · · represents an average over a series of ≈ 100 images. Figure 7 shows the distributions p(n) computed for the 6 experiments listed in Table I. These distributions typically display a power-law-like behavior with a cutoff at large cluster sizes [27,28]. As mentioned by Tallakstad et al. [27], the obtained distribution is affected by threshold values, which must thus be carefully chosen using visual inspection. Here, we focus on the variations of the distribution with the history of the system. Therefore, the most important requirement is that the image processing procedure is used consistently throughout one experiment. To avoid possible bias due to variations of illumination in the room, the experimental setup is isolated behind a dark curtain. The camera exposure time and aperture are the same for all experiments, and we use the same thresholding parameters, carefully chosen by visual inspection, for all experiments. This allows us to compare images obtained in ss 1 , ss 2 and ss 3 for a given experiment and from one experiment to the other in a meaningful way. Similarly to what we observed for the  Figure 7. Average cluster size distributions p(n) computed from steady-state images. We use the same symbols and colors as on Figure 6 histograms, it is possible to distinguish the ss 1 and ss 3 distributions from those corresponding to ss 2 (see Figure  7). This is coherent with the results of previous studies indicating that distributions are shifting towards higher clusters sizes when increasing the flow rate [27,28]. However, the ss 1 and ss 3 distributions are similar, meaning that the steady-state non-wetting cluster size distributions are history-independent. We have checked that whereas varying the threshold values affects the distributions, typically by shifting then towards lower or higher cluster sizes, it does not modify the results in terms of history independence.
The experimental boundary conditions imposed that the controlled flow variables were the total flow rate and the fractional flows. In the next Section, we turn to numerical simulations to further investigate the historydependence of the steady-state for different boundary conditions, as well as higher Ca values and different viscosity ratios M .

IV. THE NETWORK MODEL
The two-dimensional experimental porous medium is modeled by a network of tubes orientated at 45 • with respect to the overall flow direction. The tubes (or links) Overall flow direction Figure 8. Illustration of the network formed by tubes that are connected to each other through nodes where the dashed lines are intersect. One single tube is colored gray. Spherical glass beads makes the tubes as hour-glass shaped.
intersect at the vertices (or nodes) of the network with coordination number 4. The nodes are considered to have no volume, so the tubes consist of both the pore and throat volumes. The network is illustrated in Figure 8. The disorder due to the random position of the glass beads in the experiment is introduced in the model by choosing the radius r of each tube randomly from a uniform distribution of random numbers in the range [0.1l, 0.4l], where l is the length of a tube. In order to incorporate the shape of the pores in between spherical glass beads, each tube is considered as hour-glass shaped which introduces the capillary effect in the system. The network transports two immiscible fluids, one of which is more wetting than the other with respect to the pore walls. The fluids are separated by menisci and we obtain the capillary pressure p c at a meniscus inside the hour-glass shaped tubes from a modified form of Young-Laplace law [40,41], where x is the position of the meniscus and γ is the interfacial tension between the fluids. θ is the wetting angle and we consider perfectly wetting conditions, i.e. θ is either 0 • or 180 • . The flow is driven by maintaining a constant total flow rate Q throughout the network which introduces a global pressure drop. The local flow rate q in a tube with a pressure difference ∆p between its two ends follows the Washburn equation of capillary flow [40,42] where k = r 2 /8 is the permeability for cylindrical tubes. Any other cross-sectional shape of the tubes will lead only to an overall geometrical factor. Here a is the crosssectional area of the tube and µ eff (s nw ) is the volume average of the viscosities of the two phases present inside the tube. Hence, it is a function of the local non-wetting saturation s nw in that tube. The sum over p c runs over all the menisci inside the tube. The property that the fluid flux through every node will be zero is used to obtain the local pressures at the nodes. The set consisting of one equation (5) per tube, together with the Kirchhoff equations balancing the in and out flow at each node are then solved using the Cholesky factorization or the conjugate gradient method. The system is then integrated in time using an explicit Euler scheme. Inside a tube all menisci move with a speed determined by q. When a meniscus reaches the end of a tube, new menisci are formed in the neighboring tubes. In each link, a total maximum number of menisci is allowed to form. When the number is exceeded, two nearest menisci are merged keeping the volume of each fluid conserved. Here, we have considered a maximum of 4 menisci in one tube (i.e., 2 non-wetting bubbles), as it is not very likely to form a lot of menisci in one pore as seen from the experimental observations. Further details of the model and how the menisci are moved can be found in [40,43]. In order to reach the steady state in the simulation, we considered two different approaches. The conventional way is to implement the bi-periodic (BP) boundary condition, where the links at the outlet row are connected to the inlet links, so that the network acquires a toroidal topology [26]. In this case the simulation can go forever, regardless of whether one fluid percolates or not. However, in order to keep the flow going, the global pressure gradient is maintained by considering an invisible cut through the system in terms of the pressure. Since the system is closed with this boundary condition, the individual fluid volumes remain constant throughout the simulation. The non-wetting saturation S nw = V nw /V is therefore an independent parameter here, along with the total flow rate Q, whereas the non wetting fractional flow F nw = Q nw /Q fluctuates over time.
Implementing the bi-periodic (BP) boundary condition is of course impossible in the experiments. As described before, in the experimental setup, two fluids are injected at one edge of the Hele-Shaw cell through a series of alternate inlets and the opposite edge is kept open. In this case, the flow rates of the two fluids can be controlled independently. Thus, the control parameters are the total flow rate Q and the fractional flow F nw , whereas the saturation S nw fluctuates. In order to have a close emulation of the experimental ensemble, we also implement the open boundary conditions (OB) in the simulations, where the individual flow rates at the inlet links are controlled. Therefore, in OB, the system is open in the direction of total flow while we consider periodic boundary conditions in the direction perpendicular to the overall flow.

V. SIMULATION RESULTS
Simulations are performed considering networks of 256 × 256 links for BP and 128 × 192 links for OB. In order to avoid any traces from the inlets in OB, only a 128 × 128 segment of the network towards the outlets is considered for the analysis (see Fig. 10). Each link has a length of 1 mm, which is equal to the bead diame-ter used in the experiments. Most of the simulations are performed for the viscosity ratio M = 1. Three different capillary numbers, Ca 1 = 1.92 × 10 −5 , Ca 2 = 9.15 × 10 −3 and Ca 3 = 2.88 × 10 −2 are considered for BP. For OB, the capillary numbers considered are Ca 1 = 3.2 × 10 −3 , Ca 2 = 3.2 × 10 −2 and Ca 3 = 7.0 × 10 −2 . For OB we choose a similar fractional flow F nw = 0.5 as in the experiments. In BP, we run the simulation for the saturation S nw = 0.74 which we find close to the critical point for the range of parameters we considered here [32]. We will report only one set of simulations for M = 10 −4 with BP for a system with 128 × 128 links and S nw = 0.5 with capillary numbers 5.06 × 10 −2 and 1.24 × 10 −1 , as the conjugate gradient solver converges very slowly for viscosity unmatched fluids, making the simulation very computationally expensive.
To investigate the history dependence, we use a procedure analogous to the experimental one (see Sec. III A) : the simulation is started from the initial condition with a capillary number setting a constant flow rate Q 1 . Once a steady state ss 1 is reached, the capillary number is altered to a different value setting another constant flow rate Q 2 , and this new flow rate is maintained until a different steady state ss 2 is reached. Finally, the capillary number is set back to the initial value and the system is allowed to evolve towards a third steady state ss 3 . Each simulation thus involves two different capillary numbers among Ca 1 , Ca 2 , and Ca 3 , therefore 6 different simulations have been performed for each boundary condition. The steady states are compared with three different criteria -the average pressure drop, the distribution of fluid saturation over the system, and the non-wetting cluster size distribution. All the measurements are averaged over 50 to 100 configurations in the steady-state and 5 to 10 different realizations of the network.
Fluid morphologies for a typical simulation in the three steady states ss 1 , ss 2 and ss 3 are shown on Figures 9 and 10 for BP and OB respectively. Here the simulation starts from Ca = 9.15 × 10 −3 to reach ss 1 , then it is altered to Ca = 2.88×10 −2 to reach ss 2 , and then it is again turned back to the initial Ca = 9.15 × 10 −3 to reach ss 3 . Figures  9 and 10 show the distribution of saturation over the network in these three steady-states in (a), (b) and (c) respectively, where the gray-scales from black to white correspond to s nw = 1 to 0 inside a link. For BP, it is not possible to see any difference in the gray-scale saturation distributions, as the system is closed and the total saturation is conserved. However, as we will see, identifying the clusters allows us to distinguish ss 1 and ss 3 from ss 2 (see Figs. 9 (d), (e) and (f)). In OB, the saturation distributions look similar in ss 1 and ss 3 whereas ss 2 shows more non-wetting saturation than the others. This is consistent with previous numerical studies showing that the variation of saturation with Ca depends on the viscosity ratio, fractional flow and other flow parameters [26]. We then identify the non-wetting clusters using Hoshen-Kopelman algorithm [44]. As every link can be occupied with both the fluids, a clip threshold in the link saturation is considered to identify the clusters [45]. If a neighboring link has a non-wetting saturation higher than the clip threshold, the link is then considered to belong to the same cluster. The clusters are shown in (d),  Figure 13 for OB. The two different capillary numbers for each simulation corresponding to ss 1 , ss 2 and ss 3 are indicated in the plots. In BP, we measure the global pressure drop ∆P over the whole system. In OB, we measure the pressure drops ∆P i at the inlet nodes and ∆P m at the middle of the system, with respect to the outlet where the pressures are averaged over all the nodes in the corresponding row, in the direction perpendicular to the flow. We observe very similar behaviors in the pressure curves as in the experiments. In the case of BP, the system is initialized by filling the tubes with the fluids randomly at the desired saturation S nw , which will be constant throughout the simulation. This random initialization has the advantage of decreasing the simulation time significantly, since the steady state is reached faster than with an initial condition in which the two fluids are completely segregated. Due to this initial random filling in BP, the global pressure ∆P starts from a higher value and decreases with time due to the formation of clusters. Subsequently, it reaches the steady state ss 1 and ∆P fluctuates around a constant average value as seen on Figures 11 and 12. In OB, the system is initialized by saturating the network completely with the wetting fluid and then the simulation is started by injecting two fluids simultaneously through a series of alternate inlets. The flow rates of individual fluids are controlled to obtain the required fractional flow F nw . Both drainage and imbibition therefore take place at the pore level creating new menisci, which increase the pressure drop with time as seen in Figure 13. Away from the inlets, the trace of the injection channels vanishes, and a steady-state ss 1 is attained, and ∆P i and ∆P m fluctuate around constant average values. Next, as soon as the capillary number is changed to a different value, a rapid change in the pressure drops is observed for both BP and OB, and the new steady-state ss 2 is reached, characterized by different constant values in the average pressure drops. When the capillary number is again altered to the initial value to reach the steady-state ss 3 , we find that the global pressure drops change back to the initial average value. Moreover, the global pressure drops corresponding to the same capillary numbers in different simulations have the same average value, no matter from which condition they have been reached. The global pressure estimates therefore completely support the experimental observations, i.e, that the steady state only depends on the imposed flow rate, and not on the initial condition. Now, in order to find a detailed microscopic information in this regard, we measure the distribution of link saturation over the system. This measurement provides us with similar information as that of experimental gray-scale image histograms, despite being computed slighly differently. More precisely, in  Figure 13. Time evolution of global pressure drops ∆Pi at the inlet and ∆Pm at the middle of the system for OB with M = 1. A rapid change in both the pressure drops can be observed as soon as Ca is altered, but they again settle back to the initial value when the flow rate is restored to the initial one.
in different simulations are identical. A minor difference in the histograms for ss 1 and ss 3 is observed only for Ca = 1.92 × 10 −5 in BP. This is due to the appearance, at low Ca and high saturation in BP, of percolating nonwetting flow channels yielding a high non-wetting fractional flow (F nw ≈ 0.99), as shown in the inset of Figure  14 (e), while the rest of the system is immobilized.
On the histograms, we observe two distinct peaks in BP for Ca = 9.15 × 10 −3 : one at s nw 0.8 corresponding to the links mostly filled with non-wetting fluid, and the other at s nw 0.4 corresponding to the links mostly filled with wetting fluid. Therefore links can be divided into two categories : either highly saturated with the non-wetting fluid or with the wetting fluid, rather than containing a mixing of two phases. This in turn indicates the presence of large clusters at this Ca, as already observed in the fluid morphology on Figures 9 (d) and (f). For the other two Ca values in BP, the histograms are changing towards having one peak, with a more flat shape. For M = 10 −4 in BP, and in OB, the histograms also display one peak, roughly centered but whose position shifts along the x-axis with Ca. This indicates that most of the links are filled with a mixture of both fluids. Therefore it is difficult to obtain large clusters in such conditions, as observed in the fluid morphology for OB ( Figure 10).
Finally we compute the non-wetting cluster size distri- butions at different steady states. Here, the size n of a cluster is defined by the total number of links that belong to the cluster. The probability p(n) to have a n-sized cluster is then defined as p(n) = N (n)/N tot , where N (n) is the number of n-sized clusters out of total N tot clusters identified. p(n) is averaged over different configurations in the steady state and different samples of the network. In Figures 17 and 19, p(n) is plotted in log-log scale for BP and OB respectively, for different simulations. For all the simulations, we find that the cluster size distributions are identical for ss 1 and ss 3 whereas they are different from ss 2 . Distributions for the same Ca for different simulations are also the same, showing that the steady-state cluster size distributions are history independent.

VI. CONCLUSIONS
In this article we have considered the question of history dependence in the steady-state two-phase flow in porous media and presented detailed experimental and numerical investigations in this context. Experimentally, a quasi two-dimensional laboratory model consisting a Hele-Shaw cell filled with glass beads is considered, through which two phases, a gas-liquid pair with a viscosity ratio 10 −4 flows at constant flow rate. The system is allowed to evolve to a steady-state where the global pressure drop fluctuates around a constant average value. Steady-states corresponding to the same control parameters (e.g. capillary number) are attained from different initial conditions. In order to characterize the complex flow patterns in the steady-state, the gray-scale histo-  only depend on the external parameters, but do not depend on the initial preparation of the system or the history of the process. We therefore conclude that, within the range of parameters explored in this study, the steady states in simultaneous flow of two phases through porous medium are history independent.