Empirical mode decomposition of multiphase flows in porous media: characteristic scales and speed of convergence

We apply a proper orthogonal decomposition (POD) to data stemming from numerical simulations of a fingering instability in a multiphase flow passing through obstacles in a porous medium, to study water injection processes in the production of hydrocarbon reservoirs. We show that the time evolution of a properly defined flow correlation length can be used to identify the onset of the fingering instability. Computation of characteristic lengths for each of the modes resulting from the POD provides further information on the dynamics of the system. Finally, using numerical simulations with different viscosity ratios, we show that the convergence of the POD depends non-trivially on whether the fingering instability develops or not. This result has implications on proposed methods to decrease the dimensionality of the problem by deriving reduced dynamical systems after truncating the system’s governing equations to a few POD modes.


Introduction
Water flooding is a major oil production method in the petroleum industry. The idea behind this process, that is, generally applied as a secondary phase when the original energy of the reservoir has been exhausted, is to use water injection wells to generate an immiscible piston-like oil bank and push it toward the producers. However, the flow of two immiscible fluids with different viscosities in a porous medium gives rise to the so-called viscous fingering instability. The early work by Buckley and Leverett (1942) was pioneer in analyzing fluid displacements in this case. Viscous fingering is an intrinsically nonlinear phenomenon that substantially affects the performance of the recovery. A direct measure of its performance, under the effect of nonlinear instabilities, is the arriving time of the fluid injected to the well, which is known as the water breakthrough. An early breakthrough time indicates the presence of instabilities and thus also a reduction in the estimated ultimate recovery. These nonlinear instabilities are affected by the reservoir permeability distribution, fluid viscosities, capillary effects, dispersion and gravity. Thus, even when the reservoir is perfectly homogeneous (which is never encountered in practical situations), instabilities can still develop.
After the work of Buckley and Leverett, many researchers advocated the need to study the problem of fluid displacement in porous media combining theoretical, numerical and experimental approaches. Historically, studies have been focused on the recovery efficiency and on the estimation of the breakthrough time. Several studies considered the conditions under which viscous fingering appears, using tools of linear instability analysis, such as those of Koval (1963), Heller (1966), Peters et al.(1984), Chikhliwala and Yortsos (1985) and Guzman and Fayers (1997). Another approach is to numerically simulate the system of equations that governs the dynamics of the system. In this way, solutions to the problem can be obtained for a given permeability field and geometry, thus extending theoretical studies to the full nonlinear case as well as to more realistic scenarios (Jha et al. 2011;Nicolaides et al. 2015;Riaz and Tchelepi 2016). With the growth of computing power, problems tackled using this approach have become larger both in space and Edited by Yan-Hua Sun time resolution. A pioneer work in this area was done by Peaceman and Rachford (1962), where the authors solved the system of equations using a finite differences scheme. More recent works (Juanes 2005) tackle the problem using multiscale methods. Also, experimental studies were able to obtain data of the fluid in the porous medium using tomographic techniques [see, e.g., the work by Siddiqui et al. (1996)]. As these experimental studies are hard to perform, direct numeric simulations are still important to study fluid displacement, as it gives access to datasets and information that may be hard to obtain otherwise.
Using data from the numerical simulations, it is possible to apply different reduction techniques to decrease the dimensionality of the problem. The main goals of these techniques are: (1) to identify dominant modes in the dynamics, and to extract relevant information from those modes, (2) to compress the datasets and (3) to build a basis of solutions that can be used to derive a reduced system of equations for the evolution of the most energetic modes, neglecting contributions from modes with lower energy. In general, these techniques have a wide spectrum of applications, that in the oil industry go from seismic data reconstruction (Zhao and Song 2012), prediction of instabilities in oil well perforations (Lin et al. 2018), to more recently, the extraction of principal components from measured data using machine learning techniques. A pioneer work in the application of these so-called empirical mode decompositions upon a flow in a porous medium was that of Gharbi et al. (1997). In this work, the authors used a Karhunen-Loeve decomposition with the goal of characterizing the fluid displacement in the medium. The study showed that using such a decomposition can be useful to identify coherent spatial structures in the flow and that their behavior can be correctly captured by a finite set of modes which define a reduced basis. The same decomposition, used in conjunction with machine learning techniques, was used to study the fingering process in porous media (Smaoui and Gharbi 2000). Other works, such as those of Walton et al. (2013), considered the applicability of radial basis functions. Another (directly related) empirical decomposition is given by the proper orthogonal decomposition (POD), introduced by Berkooz et al. (1993) to identify coherent structures in turbulent flow. After this work, many authors used this technique for a large variety of problems from the study of solar magnetic fields (Mininni et al. 2004) to the analysis of multiphase displacements in porous media (Ghasemi and Gildin 2015).
In this work, we solve numerically the equations for multiphase flow in a porous medium using the finite volume method, and we study viscous fingering as the flow passes through a fixed number of obstacles in a channel-like geometry. The obstacles were randomly placed in order to generate the onset of viscous instabilities. We consider two phases (oil and water, to mimic the conditions of water flooding in oil reservoirs), and we vary the ratio of viscosities of the two fluids as well as the number of obstacles. Using these data, we define several characteristic lengths based on the probability density function (PDF) of the flow correlation length. We then perform a POD of the data to obtain an expansion of the datasets into a basis of empirical modes. In practice, PODs are often used to construct a surrogate model to improve the efficiency in intensive computation cases, such as in uncertainty quantification and optimization. In this light, the results presented in this work focus on the first step toward the construction of such reduced-order models, showing conditions for the convergence and useful metrics for its characterization. The main objective of this work is thus to characterize the behavior of the empirical modes and its geometrical properties, and in particular, to study the convergence of a truncated series of modes to the actual data as the parameters in the simulations are varied. Our main findings are that: (1) Properly defined correlation lengths, over the entire data as well as over individual modes of the POD, can be used to characterize viscous fingering and to identify the time of its onset. (2) The convergence of the truncated series of POD modes varies non-trivially with the ratio of viscosities, with the case with smaller viscosity ratio (and thus with less viscous fingering) requiring more modes to capture the flow profiles at fixed error, while cases with larger viscosity ratios (and thus with small-scale structures arising from the fingering) require less modes. The results have implications on proposed methods to decrease the dimensionality of the problem by deriving reduced dynamical models using truncated series of empirical modes, as the determination of the number of modes required to approximate the solutions with fixed error depends on the speed of convergence of the decomposition.

Theory and governing equations
In this work, we solve the saturation and pressure field equations in the absence of gravity (Peaceman 2000), given by where ∅ is the porosity field; K is the absolute permeability field; S is the saturation field; k r is the relative permeability field; is the viscosity; p is the pressure; q is the mass injection rate; is the density; and subscripts w and o represent the water and oil phases, respectively. Note there is no term associated to compressibility, as there is no gaseous phase in the system. The two saturation fields are related by It is useful to introduce the phase mobility for oil and water as We can also define the total mobility t and the difference mobility d , respectively, as In order to rearrange Eqs. (1) and (2), we define the average pressure p avg and the capillary pressure p c as We can also define the velocity of each phase as And the total velocity as Finally, we can define the volumetric rates of injection per unit of volume Q w and Q o as Equation (3) allows us to reduce the variables to a single saturation profile S , which in our case will correspond to the water saturation field, i.e., S = S w . Using these new variables, Eqs. (1) and (2) can be written as where f w is the so-called ratio of the wetting phase to total mobility, or fractional flow, and is given by Defining the mobility ratio M as the fractional flow becomes To solve these equations, we need to specify functions for the relative permeabilities k ro and k rw . It is known that inhomogeneities reduce de-performance of oil recovery (Khataniar and Peters 1992) and that fingering instabilities depend on the permeability distribution (Giordano et al. 1995). Based on these studies, in the present study we consider the following dependencies for each of them Although this is an arbitrary parameterization, it is frequently used. From these permeabilities, we can rewrite f w as where m is the viscosity ratio In Moissis et al. (1993), it was shown that the viscosity ratio plays an important role in controlling viscous fingering. To solve numerically these equations, we use dimensionless quantities based on a characteristic length L 0 , a characteristic volumetric injection rate Q 0 and a characteristic density 0 . Dimensionless quantities are defined from dimensional quantities (denoted with a tilde) as 1 3 with these choices, dimensional quantities can be obtained by fixing values for L 0 , Q 0 and 0 , and by multiplying all dimensionless quantities in the next sections by their corresponding factor.

Simulations and domain geometry
We solved Eqs. (14) and (15) numerically. To this end, we implemented a code using an implicit-explicit scheme (IMPES). The effectiveness of these methods can be seen (Chen et al. 2004). Spatial discretization was done using a first-order finite volume method (FVM), while time integration was done using a first-order Euler method. The code was written in C and parallelized using OpenMP. To properly capture fine structures associated with the fingering instability, a spatial resolution of N x × N y = 1000 × 500 grid points was used (in a two-dimensional channel with 2:1 aspect ratio), and the time step was determined by the CFL condition.
As mentioned in the previous section, the geometry considered is that of a channel with no flux in the Y direction over the horizontal boundaries. In the vertical boundaries, we set the left border as the injection well, and the right border as the production well. In the left border, we set a constant rate of injected water. Thus, the water and oil that flows on the right border adjust to satisfy mass conservation. In all simulations, we used a homogeneous permeability field, except in a region near the injection area where obstacles are present, with the purpose of perturbing the flow and triggering the fingering process. At t = 0 , the value of S w is zero for the entire reservoir except on the left side boundary and its vicinity, where the saturation profile S w has a value S w = 1 in x = 0 for all values of y, and decays smoothly reaching S w = 0 in x = 0.05L 0 . The oil saturation field is thus S o = 1 − S w , and as a result, at the beginning of the simulation is close to 1 in most of the domain. It is worth noting that as we are using a homogeneous permeability field in the whole reservoir (with the exception of the obstacles close to the left boundary), our fingering will not be too strong, and we will not observe stronger permeability patterns that can typically appear in the inhomogeneous case [as an example, for studies of the dependence of the fingering pattern on lateral permeability, see Waggoner et al. (1992)].
As already mentioned, in order to characterize the scale of the fingering (and later, to compare POD decompositions with varying number of obstacles) we set a group of obstacles near the injection well. These obstacles have a circular section, with permeability field inside them equal to zero. The configuration gives a simplified conceptualization of heterogeneity in oil recovery problems, representing lowpermeable zones. For practical purposes, this method is similar to introducing a small perturbation in the saturation field. In previous works (Christie and Bond 1987;Araktingi and Orr 1993), it was shown that under certain conditions, such as in the presence of an uncorrelated permeability field, the final behavior of the flow is independent of the method used to excite the instabilities. Thus, we choose the position of the center of the obstacles in the X direction to be the same for all the obstacles, and we do not vary it between simulations. The distances between the centers of the obstacles in the Y direction are given by L 0 ∕(1 + obs) + i,i+1 , where obs is the number of obstacles, and i,i+1 is a random displacement between the obstacle i and the obstacle i + 1 . The random displacement i,i+1 has a normal distribution centered around zero, and thus has small values that can be positive or negative. The radius of these obstacles also follows a normal distribution with mean r mean = 0.05L 0 . In Fig. 1, we show the detail of the permeability maps used in all simulations. Figure 1a shows the whole computational domain for the simulations with five obstacles, while Fig. 1b-e only shows the permeability maps in the subregions with obstacles, respectively, for the simulations with 1, 2, 3 and 4 obstacles. Instabilities start soon after the injected phase goes through the region with the obstacles. In different simulations below, the parameters that we vary are the number of obstacles ( obs ) and the viscosity ratio m . For the viscosity ratio, we use either m = 5 , m = 10 or m = 15 . For each of these values of m , we then perform five simulations with 1, 2, 3, 4 or 5 obstacles. The incoming flow has a fixed value Q w = Q 0 , and the capillary pressure p c is zero. As mentioned above, the permeability K is equal to unity in the entire reservoir except in the obstacles. The porosity ∅ is also equal to unity in the entire reservoir. Finally, the viscosity of the water is set to w = 1 (in dimensionless units as discussed above). The data from the simulation are saved at regular intervals of Δt = 0.2∕Q 0 . A full list of the simulations and of the parameters we vary is shown in Table 1, where 'Name' indicates the label used for each simulation, m is the viscosity ratio, and 'obs' indicates the number of obstacles distributed randomly in each run. Note the simulations are separated in three groups, each with the same value of m and varying number of obstacles.

Proper orthogonal decomposition
The POD (Berkooz et al. 1993) has the main goal of projecting a high-dimensional dataset (usually obtained from numerical simulations, field observations or from experiments) into an optimal basis of orthogonal modes. The expansion is optimal in the sense that the modes are ordered in decreasing energy (where the energy is associated to the power contained in that mode), and with the fastest possible convergence. As the basis obtained from the POD is empirical and specific for a given dataset, the POD is useful to identify coherent structures and to find bifurcations in the dynamics of the problem. In particular, the POD has often been used as a first step to construct reduced systems of ordinary differential equations from a set of partial differential equations (PDEs), by doing a Galerking projection of the PDEs into the empirical basis truncated to a few modes.
Given a scalar field u(x, t) , we can decompose it into orthonormal modes using the POD. As a result, we obtain a basis of spatial modes u � i (x) , and a set of orthogonal temporal coefficients a i (t) , such that As already mentioned, the scalar field u(x, t) is usually obtained from experiments or simulations. In our case, the dataset will be the result of the numerical simulations described in Sect. 3.1. For the moment, we can think of the scalar field as a discrete set of arrays where u i (i = 1, 2, 3, … , N) are snapshots of the scalar field at fixed times, and N in the total number of time steps in the simulation. Given u(x, t) , we can build a linear operator where the overline denotes the complex conjugate. Assuming u is a compact operator, we can rewrite these relations as an eigenvalue problem for eigenfunctions k and k , both with eigenvalue a 2 k , given by The scalar field u(x, t) can then be expanded in terms of the solutions to Eq. (34) as Here, k are the amplitudes of each mode (with a 2 k associated to their energies), and k (x) and k (t) are the empirical orthogonal modes which spectrally decompose u(x, t) . The spatial modes k (x) are usually called the 'topos,' while the temporal eigenfunctions k (t) are known as 'cronos.' To apply this decomposition to experimental or numerical data (which are discrete in space and in time) of saturation, we can use the so-called snapshot method as described by Sirovich (1987), which in practice is a decomposition in the time domain. In this work, the saturation field S corresponds to the results of the numerical simulations described in Sect. 3.1, saved to N snapshots S i every dt intervals in time. Thus, Eq. (31) can be rewritten as , it is the saturation at fixed time and as a function of the two spatial coordinates). The base saturation is then defined as the average of the ensemble, We can use the base saturation to rewrite each of the snapshots S i as a mean field plus fluctuations, as We then define the correlation matrix between any two snapshots p and q as where the brackets denote the average over space. Using this matrix, the equations above are equivalent to solving the eigenvalue problem given by The eigenvectors of this equation are the temporal coefficients of the expansion b i = b i 1 , … , b i N , which are equivalent to the coefficients a i (t) defined above. The spatial modes can then be obtained as (36) S(x, y, t) = S 1 , S 2 , S 3 , ..., S N , The equivalent of Eq. (30), giving the reconstruction of the saturation field at any point and at any time, finally becomes where for brevity we replaced the index k (corresponding to the discrete time snapshot) by the time label t. The larger the number of modes we use in the reconstruction, the more accurate it will be.
Once the decomposition is obtained in this way, it can also be used to project the governing equations into a few modes (the modes that dominate the dynamics), allowing for a reduction in the dimensionality of the problem. In other words, the set of PDEs in Sect. 2 can be truncated to a finite number of ordinary differential equations. In the Appendix, we describe the procedure used for that purpose. Figure 2a shows a typical saturation profile for m = 5 and one obstacle at time t = 53∕Q 0 (run I). Water is injected from the left, pushing the oil phase to the right. As the water phase passes through the obstacles, it can develop a viscous fingering instability which, in all runs, appears for times larger than t ≈ 17∕Q 0 . A vertical cut across any fixed value of L x to the right of the obstacle then results in a saturation profile which is inhomogeneous and displays alternating maxima and minima.

Global evolution
With the aim of characterizing the development of viscous fingering along the flow direction, we can look at characteristic scales in the flow and at their time evolution. A usual definition for the characteristic scale of a flow is given Light colors indicate maximum value of saturation while dark colors indicate a minimum value. As an illustration, white and gray arrows along a vertical cut of S w (water cut) show regions with saturation above and below the average. b A vertical profile of the saturation along the vertical line indicated in a. Gray and black arrows correspond to the same regions as, respectively, gray and white arrows in a by correlation length L S corr , which here we compute in the Y direction. Simply put, the length L S corr is obtained from the correlation function of the saturation field S = S w , as a function of y and averaged over all values of x . Figure 3 shows this length for simulations with m = 10 and varying number of obstacles. As expected, the correlation length is dominated by the large-scale structures. While at early times the behavior of L S corr is the same for all runs (before the water phase passes through the obstacles), at later times L S corr depends on the number of obstacles, reaching a value that decreases with increasing number of obstacles. To have a better understanding of how the correlation length varies as a function of the number of obstacles, we average L S corr over time using a temporal range over which it remains approximately constant (see the squared region in Fig. 3). Figure 4 shows the values obtained from the average, ⟨L S corr ⟩ , as a function of the mean distance between obstacles, for m = 5 , m = 10 and m = 15 . For large distances between obstacles (i.e., for a few obstacles in the domain), ⟨L S corr ⟩ converges to the same value independently of m , confirming that ⟨L S corr ⟩ is dominated by the contribution of the large-scale structure in the flow associated to the distribution of the obstacles. However, for larger numbers of obstacles (or, equivalently, for smaller distances between obstacles) ⟨L S corr ⟩ also depends on the value of m (taking smaller values for larger values of m ) and thus appears to become sensitive to the contribution of the small-scale structures associated with the fingering instability.
This correlation length will be useful to identify the characteristic scale of each POD mode in Sect. 4.2. However, to study the fingering instability from global data, we need a different definition for a characteristic scale that is more sensitive to the small-scale correlations associated with the viscous fingers. We thus introduce now a 'crosslength' L S cross also in the Y direction and over the saturation field S . We will see this magnitude is more useful to quantify the size of small-scale structures, and therefore, to identify the onset of the instability. The length L S cross is defined as a cross-length, i.e., as the length we have to displace in the Y direction to cross a given value of the concentration. The detailed procedure is as follows: (1) At each time, we compute the mean concentration S . (2) For each value of x , we move across the Y direction and compute the distance between all points with S =S , as shown in Fig. 2b (i.e., we compute the distance between 'crossings' of S with its mean value). (3) The process is repeated for different times in the simulation. (4) Finally, we compute the PDF of all these distances. From the PDF, we can compute different moments (as, e.g., the mean cross-length, its median or its deviation).
The PDF of L S cross for different simulations with three obstacles and viscosity ratios m = 5 , m = 10 and m = 15 are shown in Fig. 5. The arrows in each figure indicate the separation distance between obstacles, as well as between the obstacles and the wall of the computational domain (note obstacles are distributed randomly, so the separation is not the same for all obstacles). The peaks in the squared region remain approximately the same as m is varied; they correspond to large values of L S cross , and they lay close to the distance between obstacles. Thus, we can conclude that this part of the PDF of L S cross is associated with large-scale geometrical features in our domain. However, for smaller values of L s cross (see the white region on the left of Fig. 5) the shape of the PDF changes as m is increased. In particular, for larger values of m , the PDFs display a broader peak around L S cross ≈ 0.05∕obs (where obs is the number of obstacles). As corr ⟩ as a function of the mean distance between obstacles L 0 /obs in the runs with larger values of m , we have more viscous fingering, we will see that these changes in the peak of the PDF are associated with the growth of small-scale structures in the flow, resulting from the fingering instability.
From the observation in Fig. 5 that a part of the PDF of L S cross is sensitive to the value of m , we can build a characteristic length associated with the instability as the median of L S cross , denoted as L S cross . Figure 6 shows this quantity as a function of time, first for simulations with fixed m = 10 and with different numbers of obstacles in Fig. 6a. The median of the cross-length decreases with the number of obstacles, just as the correlation length also does. However, a sudden drop can be also seen at intermediate times in all cases. Figure 6b shows a detail of the median for simulations with just one obstacle, with m = 5 , m = 10 and m = 15 . In the simulation with m = 5 , L S cross remains constant in time, once the flow of water passes through the obstacles. However, for all other values of m , L S cross shows again a sudden drop in its value at later times. Visual exploration of the simulations indicates that the time of the drop corresponds to the onset of the viscous fingering instability. This can be understood, as smallscale correlations appear in the data once the fingering starts to grow (see Fig. 7).
To further illustrate this effect, the squared regions in Fig. 6b define two times for each simulation: The left boundary of the squared region, corresponding to the beginning of the drop in L S cross , defines a time t 0 . The right boundary defines another time t 1 . At these two times, Fig. 7 shows the saturation field in the simulation with one obstacle and m = 10 at t 0 . It can be seen in this figure that structures associated to fingering appear at time t 0 .

POD analysis
Applying the POD method to the simulations, we can decompose them into orthogonal modes, study each mode separately and reconstruct the saturation field using a reduced number of modes (both topos and chronos) as The elements of the basis, the spatial functions i and the temporal coefficients b i , were obtained using the snapshot method from the correlation matrix C pq . Figure 8 shows a simulation with m = 10 for one obstacle, together with a reconstruction using 10, 20 and 30 modes, respectively. As expected, if we increase the number of modes, the reconstruction is more accurate. We need to set a criterion to estimate how accurate the reconstruction is for a certain number of modes. As an attempt to accomplish this, we calculate the energy spectrum as a function of the mode number (i.e., the energy contained per POD mode). Figure 8 shows this spectrum normalized for runs with m = 5 , m = 10 and m = 15 as a function of k (the mode number), and for simulations with three obstacles. We found that for different numbers of obstacles used, the results practically did not change. To make an estimation of the error made in the reconstruction as a function of k , we This energy quantifies the amount of total energy gathered in a reconstruction with k modes, while the residual is the energy lost, and is a usual measure to quantify the quality of the reconstruction in PODs. The result is shown in Fig. 9b.
The green line indicates, as a reference, a fixed value of energy of 0.98E 0 (where E 0 is the total energy in the original field). This means that, in this case, we need five modes to reconstruct (up to an error in the energy of 2%) the simulation with m = 15 . Note that the larger the viscosity ratio m , the larger is the number of modes we need to attain the same accuracy. The results shown in Fig. 9 were obtained for simulations with three obstacles; however, the results are similar regardless of the number of obstacles.
To explicitly show the dependence with m of the number of modes needed to reach 98% of E 0 , we plotted this number in Fig. 10 as a function of m . The decaying behavior as m grows is counter-intuitive. We know that when m increases, the structures associated to viscous fingering also grow, and that more small-scale features are present in the flow. Thus, a larger number of POD modes can be expected to be required to reconstruct the field with a prescribed error when m is increased, unlike what is observed in the figure.
The reason for this behavior is as follows: As m increases, the number and amplitude of the structures associated to the fingering indeed increase. However, when m increases, the gradient separating the water phase and the oil phase also becomes smoother. This permits a better reconstruction using less modes. Figure 11 shows a profile of the advancing front for simulations with one obstacle and with m = 5 , m = 10 and m = 15 , respectively. Each one of them shows the advancing front in the simulation, and the reconstruction of the same front performed with 8, 16 32 and 64 modes. We can see how as the number of modes increases, the reconstruction becomes a better approximation to the advancing front in the simulation. However, for the sharper front in the simulation with m = 5 , more modes are required to get a reasonable approximation: While 64 modes in the simulation with m = 15 are very close to the actual profile, the same number of modes displays strong fluctuations near the front in the m = 5 case.
In the same way as we calculated the correlation length over the saturation field L S corr , we can now calculate the  Figure 12 shows the correlation length over the Y direction as a function of the mode number k (i.e., for the kth topo mode), for different numbers of obstacles and for m = 15 . We can see that this correlation length is also ordered as a function of the number of obstacles in a descending way. Between k = 50 and k = 80 , and specially for low number of obstacles, the function is compatible with a power law, which suggests that for a range of scales the viscous fingering process may become independent of the scale and selfsimilar. A certain scale independence of the process can be expected from the dynamic equations in Sect. 2, which when we made dimensionless using the length L 0 and all other units in Eqs. (23) to (29) are not explicitly dependent on the length L 0 nor on any dimensionless number based on L 0 . Note that as m increases, the maximum of the corresponding curve (i.e., the maximum correlation length per mode) becomes smaller, in both cases (i.e., for both obs = 1 and 5). This behavior indicates that the spatial representation of the viscous fingering seen in the simulations also emerges in the individual topos associated to these runs, and in particular, in their correlation lengths. In Fig. 13, it can also be seen that the maximum correlation takes place in topos with smaller k as m increases. This effect is illustrated in better detail in Fig. 14, which shows the mode number k corresponding to the maximum of L t corr as a function of m in Fig. 14a, as well as the average value of L t corr as a function of the distance between obstacles in Fig. 14b.
In Sect. 4.1, we also introduced the cross-length over the saturation field S, as a way to quantify the onset of the fingering process. We can now do the same analysis over each spatial mode, by computing a cross-length L t cross (k) for each topo (each labeled by the index k). The length L t cross (k) is calculated following the same procedure as before for the L S cross length, only for each topo instead of for the total concentration S . This will allow us to identify which individual spatial modes capture the fingering process, and which ones are associated to the large-scale (and smoother) flow. Figure 15 shows this length as a function of the kth mode for the different simulations. In Fig. 15a, we present the value of the cross-length for multiple simulations with m = 10 and with different numbers of obstacles. The dashed lines demark the values for L S cross obtained in Sect. 4.1 for the entire saturation field S before the fingering process is started (and for each simulation). In other words, it indicates the characteristic cross-length of the bulk (large-scale) flow. Note the first three modes in all cases have cross-lengths similar to those obtained for the entire saturation field. This indicates the first spatial modes of the POD are associated with the large-scale flow. For k = 4 and larger, L t cross (k) quickly drops to values that are comparable to those seen in Fig. 6 (for L S cross ) after the fingering starts. Thus, these modes have cross-lengths comparable to those of the viscous fingers and are the modes that capture the dynamics of this process. Figure 15b shows how this changes as we change m from 5 to 15, for only one obstacle. The behavior of L t cross (k) is much alike to the one already described, but the number of modes associated with the large-scale flow decreases as m increases.

Conclusions
We presented numerical simulations of the fingering instability in multiphase flow in porous media, varying the ratio of viscosities between the oil and water phases, as well as the number of obstacles used to trigger the instability. Our main goal was to characterize the process of viscous fingering using global correlation lengths, as well as an empirical mode decomposition. To this end, we studied the evolution of the saturation correlation length L S corr and introduced a new definition for a characteristic length based on the distance between crossings of the saturation with its mean value and the cross-length L S cross . We also performed a POD decomposition of the simulations and studied its convergence as well as the correlation length L t corr and cross-length L t cross of each topo. We showed that, when computed on the entire (nondecomposed) saturation field, the correlation length L S corr is dominated by the number of obstacles in the flow, while the median of the cross-length, L S cross is a better indicator of the onset of the fingering instability. However, when applied to each individual mode of the POD, the correlation length L t corr also becomes sensitive to the growth of small-scale features associated to the instability. Moreover, the cross-length of each topo L t cross can also be used to distinguish modes associated with large-and small-scale features in the flow.
We showed that when attempting a reconstruction of the saturation field using a finite number of POD modes, the convergence is non-trivially dependent of the value of the viscosity ratio m . While in all cases a small number of modes suffice to obtain the saturation to a good approximation (with errors smaller than 2% or 5%), more modes are required for small values of m (which do not display strong viscous fingering) than in the cases with larger values of m (which display fingering, and thus, small-scale features). This result is associated with the sharpness of the boundary between water and oil in the former case, which requires more modes to capture the stronger gradients in the saturation field. Also, this result is important when reduced dynamical systems for multiphase flow are derived from empirical modes, e.g., by truncating the governing equations to a few POD modes, as the number of modes required to properly approximate the solutions depends on the speed of convergence of the series, and thus on the viscosity ratio.  Fig. 6) for the complete saturation field S before the appearance of the fingering (i.e., the cross-length of the large-scale flow). b Length L t cross (k) as a function of the mode number, but now for simulations with m = 5 , m = 10 and m = 15 and for only one obstacle. The dashed line indicates the length L S cross measured for the entire saturation field before the fingering process started whose parameters can be different from the former one. To this end, we must do a Galerkin projection of the full set of differential equations in Sect. 2.
To illustrate the procedure, we consider Eq. (14) (other equations can be projected following the same steps), and in particular we provide details for the projection of one term in this equation. The projection for this term (as well as for any other term in the PDEs) can be done using the orthogonality of the basis obtained from the POD. Starting from Eq. (14), we multiply this equation by the mode k to obtain and we integrate over space Looking only at the first term, given by we can expand S using the POD modes as Here M < N is the number of modes used in the truncation, i are the spatial modes obtained from the POD, and i are new coefficients which are unknown, and for which we want to derive a set of ODEs which will prescribe their evolution. From this evolution, the concentration can be reobtained at any moment by using Eq. (46). We can also expand p avg as using the same spatial modes i obtained for S. As before, the coefficients i are unknown. Now we can rewrite Eq. (45) using Eqs. (46) and (47)  By rearranging the terms in Eq. (48), we get where A k i , B k ij and C k ijl are given by the following integrals The A k i , B k ij and C k ijl coefficients can be calculated numerically from the POD spatial modes and stored. Repeating this procedure to the other terms in Eq. (44), a set of algebraic equations for the coefficients i and i is finally obtained. Further repeating this procedure for Eq. (15) results in another equation for i and i , but with a term dependent on ̇i (i.e., its time derivative), which follows from the time derivative of S. The resulting new system of ODEs is a truncation of the full system of PDEs that can be solved for a finite number of modes M.