Onsager-Symmetry Obeyed in Athermal Mesoscopic Systems: Two-Phase Flow in Porous Media

We compute the fluid flow time-correlation functions of incompressible, immiscible two-phase flow in porous media using a 2D network model. Given a properly chosen representative elementary volume, the flow rate distributions are Gaussian and the integrals of time correlation functions of the flows are found to converge to a finite value. The integrated cross-correlations become symmetric, obeying Onsager's reciprocal relations. These findings support the proposal of a non-equilibrium thermodynamic description for two-phase flow in porous media.


I. INTRODUCTION
Athermal fluctuations occur in a number of phenomena in nature, important to biology, chemistry and physics [1][2][3]. Currently, an active effort is taking place to better understand the statistical physics of such systems and its use is realized for a growing number of research areas [1,[3][4][5][6][7][8]. A particular example is granular materials, which constituents are macroscopic. In the absence of an external driving force the material will stay in its current configuration, sharing some properties with non-ergodic systems [3]. However, when a granular material is exposed to an external force, a great number of states may be visited resulting in solid-or fluid-like behavior as a response to that force.
One area that seems to have not been analyzed in those terms, are flows driven through porous media. Such flows are important for numerous geological and technical processes, say in oil production, CO 2 sequestration, water transport in aquifers, or heterogeneous catalysis. An important class of porous media flows is the simultaneous flow of two immiscible fluids. In such a system, clusters of the two fluid phases, travelling through the porous media, are constantly forced to split and recombine. Thus, the fluid configuration in the pore space changes, leading to fluctuations in the flow rate of each phase (fractional flow rate), as well as in the total flow rate. These fluctuations are of a mechanical nature, different, but analogous to thermal fluctuations on the molecular level. The fluctuations appear on a mesoscopic scale much larger than the molecular scale of statistical thermodynamics, yet the mesoscopic scale which is defined by the pore sizes of the medium is very small compared to the overall system. In the most extreme cases the pores can be in the nanosize regime, while the system of interest, for instance in chalk oil reservoirs [9], has geological dimensions.
Our long-term aim is to find a non-equilibrium thermodynamic description for such flow systems. The art is then to define a suitable representative elementary volume (REV), where the essential assumption of local equilibrium, as expressed by the ergodic hypothesis, and microscopic reversibility holds. The statistical foundation of the theory was spelled out a long time ago [10]. Experimental [11] and computational [12,13] evidence, exist now, that the ergodic hypothesis can be expected to hold for immiscible two-phase flow in two-dimensional porous media of a minimum number of links.
Here, the aim is to move one step forward, and examine the idea of time-reversal invariance or microscopic reversibility of fluctuations [10,14]. Thus, our interest is the time-correlation functions of the flows. On the molecular scale, thermal fluctuations have correlation functions that are connected to transport coefficients. This is formulated in the Green-Kubo relations, which are frequently employed in molecular dynamic simulations. The method goes back to Onsager's regression hypothesis [15,16], which says that the decay of molecular fluctuations are governed by the same laws as the relaxation of macroscopic non-equilibrium disturbances. For the Onsager reciprocal relations [15,16] to apply, microscopic reversibility must hold. The idea of the present work is to apply Green-Kubo-like relations to the fluctuations in a REV of a network model. The Green-Kubo relations for the molecular level apply to global equilibrium. In the present approach we will use similar expressions, but for fluctuations on the mesoscale, thereby extending or expanding the Green-Kubo-scheme. We shall see that the system models a time reversal invariant process and that arXiv:2001.04332v1 [physics.flu-dyn] 13 Jan 2020 the integrated flux correlations satisfy Onsager symmetry. A similar approach to fluctuations in hydrodynamic dispersion was taken by Flekky et al. [17].
In the present study, we model the immiscible twophase flow through a porous media using a hexagonal lattice, where the links represent pore throats and have a distribution in link radii. The steady flow of two incompressible and immiscible fluids is driven by a constant pressure difference across the network, leading to a steady state with fluctuating fluid flow. The flow properties fluctuate around well-defined averages and system is in a non-equilibrium steady-state on the network level of description. A correlation of the two flows is thus unavoidable. But what is the nature of this correlation?
The answer will have an impact on how we may proceed with a theoretical description of the flows.
If one considers a steady state of the immiscible twophase flow, as we will in our model, we shall see that the fluctuations become Gaussian around a steady mean.
Hence, the concept of the REV at steady state is highly relevant and important for how we build a theory that can help us understand transport in porous media.

II. MODEL
The transport of the two immiscible fluids through a two-dimensional porous medium is represented by a dynamical pore network model [18,19]. This model has been in development over two decades and has a record of explaining experimental and theoretical results in steady and transient two-phase flow in porous media [11][12][13][18][19][20][21][22][23]. In this model the two fluids are separated by interfaces and are flowing in a network of links which are connected at nodes. The network has a honeycomb structure as illustrated in Figure 1, with equally long links and a distribution of link radii. The radii are drawn from a uniform distribution in the interval 0.1L to 0.4L, where L is the length of the links. The flow rate q ij inside a link connecting nodes j and i is given by: Here p j and p i are the pressures at nodes j and i, c ij is the capillary pressure, and g ij is the conductivity of the link. The links have an hourglass shape, thus the capillary pressure is a function of the interface positions, z ij : Here γ is the surface tension, r ij is the radius of link ij, and The effect of the χ-function is to introduce zones of length βr ij at each end of the links where the pressure discontinuity of any interface is zero. The conductivity of the link, g ij , contains a geometrical factor and the effective viscosity of the link: Here, r ij is the radius of the link and the viscosity is defined as were started with a random distribution of the two liquid phases, and were propagated at least 300000 time steps to allow for the system to reach steady state. Statistics for the time correlation functions were collected for 9.7 million time steps. The length of a link in the network was set to 1mm. We report results for each set of parameters as averages of at least 30 runs using different starting configurations of the two fluids. Volume flow rates and velocities refer to network averages. The properties of the steady state flow are determined by the pressure drop across the system, ∆P , the total wetting saturation of the network, the surface tension, γ, and the viscosity of the two fluids. We have ensured that the same steady state flows averages are obtained from different initial distributions of the two fluids including an initial configurations where the two phases are completely separated (i.e. each phase is in a single connected cluster).
Furthermore, it has been tested that a different link radii configuration, drawn from the same uniform distribution does not change the steady flow averages nor the appearance of the time-correlation functions computed in this study.
We investigated time correlation functions and their long-time convergence for two choices of the parameter set µ w , µ n and γ, which are viscosities of the wetting and the nonwetting phase and the surface tension, respectively.
Here, the capillary number is defined as Ca = vµ n /γ, with v being the seepage velocity. The case was chosen to represent flow where the two fluids are interchangeable with respect to their viscous dissipation.
Case (B) had µ w = 5µ n (µ w = 10 −3 Pa·s) and γ = 0. This case is typical at high flow rates where the contribution from γ essentially can be neglected, and the capillary number goes to infinity. It can be considered as a limiting case, chosen to elucidate the behavior when viscous forces dominate and the surface tension is negligible. Data was collected for wetting phase saturation S w =0.25, 0.5 and 0.75. Here, S w is the volume fraction of the wetting phase of the total pore volume in the network.

III. RESULTS AND DISCUSSION
We report first that the fluctuations in flow velocities are Gaussian when a suitable representative volume (REV) is chosen. We proceed to give the structure of the time correlation functions for the REV. The results for what we will call from now the Green-Kubo coefficients for the network follow from this.

A. Fluctuations
In case (A) the resistance is determined by the positions of the interfaces in the links only, as the two phases have the same viscosity. In case (B), the resistance to flow in link i is inversely proportional to the effective viscosity. The positions of the interfaces are then irrelevant as there is no surface tension and hence, no capillary pressure. A typical example of fluctuations in the total volume flow Q for case (A), with a pressure gradient (∆P/∆x of 100 kPa/m and surface tension 30 mN/m, is shown in Figure 2. By plotting the statistical frequency of the flow rate or the seepage velocity, we obtain a Gaussian distribution. This is shown in Figure 3 where we plot the statistical frequency of the fluid velocity (counts) on a logarithmic scale vs. sgn In such a plot a Gaussian distribution appears as a triangle, and this behavior is very well followed by the data. There is only a very small asymmetry in the distribution which is to be expected as the fluid velocity cannot be less than zero. A regular plot of the distributions is presented for the seepage velocity, and the velocities of the wetting and non-wetting phases in Figure 4. The distributions are shown for two network sizes of case (A), with 30×20 links and one twice the size of the former, 60×40 links. The shown distributions are normalized with the area, and the variance of the larger network is 1/2 the width of the smaller network. So, in spite of the apparent noise seen in Figure 2, one obtains the distributions in Figure 4, which has its analogue in a molecular picture, basic to thermodynamic equilibrium properties. This allows us to proceed with the next step, and construct time correlation functions for the meso-level.

B. Network Size and Representative Volume
Ideally, the system size of the simulation is sufficiently large and representative of the statistical ensemble. In this the entropy and other thermodynamic properties are proportional to the system size, (i.e. they are extensive [24]). With the Gaussian nature one may expect that the inverse variance 1/σ 2 of the fluctuations is proportional to the system size, i.e the area A. This relation is plotted in Figure 5 for network models of dimension 30×20, 60×40 and 120×60 links. It shows that this requirement is well met by a system with low Ca, but systems with higher Ca may be more susceptible to possible size effects. The size of the REV will be system dependent, see Savani et al. [13]. But in the present cases (A) and (B), a REV can be defined, for a range of Ca, different for the different cases. The results for the REV complies with the meso-level analog we are seeking.

C. Time Correlation Functions
With a well-defined REV, and with Gaussian fluctuations established, we can proceed to define the time correlation functions C RS for the fluctuating quantities R and S at the meso-level: where the brackets · · · indicate ensemble averages. The fluctuation from the mean, δR, is defined as and δR = 0.
(8) Figure 6 shows the time correlation functions of the total flow rate Q for different choices of the surface tension. After a rapid decay on a short timescale (below 10 −3 s) there is a slower, logarithmic decay which is more pronounced for larger values of γ (between 1 ms and 100 ms). These two regimes are followed by a slow long time decay. The rapid decay appears on the time scale that  correspond to the time necessary to evolve the flow by one average link volume. As shown in Figure 6, the decay is somewhat faster for smaller surface tensions as the total flow velocity is higher. The regime of the logarithmic decay is within the time of evolving the flow by the total volume of the network, and is more pronounced for higher surface tensions and thus higher capillary pressures in the pores. Hence, the decay corresponds to parts of the flow that is slow moving or frustrated. These are the regimes of interest here. They contain the relative movements of the two flows in terms of their mutual displacement.
It is interesting to note some similarities with time correlation functions of glass [25] or yield-stress fluids [26]. In these cases, the autocorrelation functions, like the selfscattering function, can be fitted to a function of the form: We attempted a fit of F (t) to the autocorrelation functions of the total flow, see Figure 7. Satisfying fits could be obtained with the exception that the local minimum and maximum at around 10 −3 s and in some cases the flat top (at times < 10 −4 s) are not well described. Fit parameters for the different choices of γ are summarized in Table I.

D. Convergence and Symmetry
The Green-Kubo method employs integrals of suitable time-correlation functions C RS (as defined in equation 6) to compute coefficients, L RS : The convergence of the integral over the time correlation function of the total flow is shown in Figure 8. As in molecular dynamics, where the Green-Kubo method is normally used, the convergence is slow and statistics has to be collected over long time-scales and/or multiple trajectories to achieve convergence of the integral when τ is approaching infinity. We computed the integrals for the autocorrelation and crosscorrelation functions of the wetting and nonwetting phases, with the indexes i, j = n, w referring to the nonwetting and wetting phase, respectively. The results are listed in Table II for case (A), where the surface tension is 30 mN/m and the fluids have the same viscosity, for three different choices of ∆P , the pressure difference across the network. Table IV summarizes results for infinite capillary numbers (case B), where the surface tension is zero but the fluids have different viscosities, for different choices of the saturation. In all cases (A) and (B) we find that the integral cross correlations obey the Onsager reciprocal relation, Λ ij = Λ ji , within the statistical error in the simulations. This result is new for a meso-level description, like the one used here, and is encouraging for the overall aim; to create a non-equilibrium thermodynamic description on for the macroscopic level. The finding applies to a well defined REV, for which we have a Gaussian distribution of fluctuations, analogous to the corresponding distribution on the molecular level.
It is interesting that the cross coefficients are all negative. This makes sense for network flow where one component cannot advance faster (on average) than the mean flow, unless the other component advances slower (on average).
The Λ ij s for case (B), where the surface tension is zero, show an extreme limit property, because the crosscorrelation functions obey Λ ww Λ nn −Λ wn Λ nw ≈ 0. A singular matrix of coefficients is the essence of complete coupling of the two fluids' flows; they are linearly dependent. For all choices of saturation Λ ww = ζ 2 Λ nw , where ζ is a constant [14]. On the other hand, for case (A), where the surface tension differs from zero, a deviation of this dependency is observed. The linear dependence of the fluxes in case (B), can thus be associated with lack of capillary forces. This can be understood in the following way: in case (B) the variation of mobility in a given link is a function of the saturation in the link only. However, if the mobility in one link is increased it has to decrease elsewhere. In contrast, for case (A) the variations in link mobility depend on the interface position and a change in the link mobility can take place without affecting the mobilities of other links.
The value of ζ for case (B) can be deduced by looking at the coefficients in Table IV. Within the accuracy, we find ζ 2 ≈ 20 or ζ = 4.5 ± 0.5 for all S w . The value is close to the ratio of fluid viscosities, which will describe the dissipation. Table IV show a dependence on the saturation, decreasing in value as the saturation is increasing. Here, the wetting fluid is the more viscous fluid and the increase in saturation reduces the effective permeability. The coefficients show also a dependence on the pressure drops across the network (see Table II), increasing with higher pressure drops. This is consistent with a higher effective permeability at higher pressures. In fact the dependence of the volumetric flow is non-linear as can be seen in Table III. A non-linear dependence of a flow rate on the pressure difference is a well-known phenomena in immiscible two-phase flow [27].

All coefficients in
To investigate the origin of the Onsager symmetry in more detail we examined the cross correlations in Figure 9. For molecular systems, one can find transport coefficients using the Green-Kubo method, see [28] and Onsager reciprocal relations apply given time reversal invariance. We have found that time reversal invariance applies also on the mesoscopic level, as formulated by: C AB (τ ) = C BA (τ ). This equality is pictured in Figure 9. There is agreement, except for very small times where the contribution to the integral is negligible. Moreover for case (B), the deviation from symmetry at small times is attributable to the finite system size. It is vanishing for larger network sizes. This is shown in the upper three panels of Figure 9.

IV. CONCLUSIONS
Our investigation of time correlation functions has revealed interesting parallels between the time correlation functions of two immiscible fluids in a porous media, those observed for glass and stress-yield fluids, and those for molecular fluctuations. A network with incompressible fluids has been used as model for the porous medium, but the findings should not be restricted to this. We have been able for the first time to find Onsager symmetry in athermal fluctuations on the meso-level. The symmetry of the coefficients implies time reversal invariance or microscopic reversibility of fluctuations also on the meso-level, in agreement with recent experimental [11] and computational evidence [12,13]. Time reversal invariance is here understood as C AB (τ ) = C BA (τ ), holding for all time scales except very short times.
We found that the structure of the time correlation functions depends on the surface tension. Integrals over auto and crosscorrelation functions of a REV, were found to converge and the integrals of the crosscorrelation functions essentially obeyed Onsager's symmetry. The coefficients obtained in this manner may have a relation to  porous media permeabilities. Further research of time correlation functions to compute transport properties of immiscible-two phase flow is therefore encouraged.