Viscous Crossflow in Layered Porous Media

Abstract

We examine the effect of viscous forces on the displacement of one fluid by a second, immiscible fluid along parallel layers of contrasting porosity, absolute permeability and relative permeability. Flow is characterized using five dimensionless numbers and the dimensionless storage efficiency, so results are directly applicable, regardless of scale, to geologic carbon storage. The storage efficiency is numerically equivalent to the recovery efficiency, applicable to hydrocarbon production. We quantify the shock-front velocities at the leading edge of the displacing phase using asymptotic flow solutions obtained in the limits of no crossflow and equilibrium crossflow. The shock-front velocities can be used to identify a fast layer and a slow layer, although in some cases the shock-front velocities are identical even though the layers have contrasting properties. Three crossflow regimes are identified and defined with respect to the fast and slow shock-front mobility ratios, using both theoretical predictions and confirmation from numerical flow simulations. Previous studies have identified only two crossflow regimes. Contrasts in porosity and relative permeability exert a significant influence on contrasts in the shock-front velocities and on storage efficiency, in addition to previously examined contrasts in absolute permeability. Previous studies concluded that the maximum storage efficiency is obtained for unit permeability ratio; this is true only if there are no contrasts in porosity and relative permeability. The impact of crossflow on storage efficiency depends on the mobility ratio evaluated across the fast shock-front and on the time at which the efficiency is measured.

1 Introduction

Alternations of parallel, continuous layers of different lithologic and physical properties are a ubiquitous type of geologic heterogeneity observed at many different length scales, including lamination (millimeter-thick layers), bedding (centimeter- to meter-thick layers) and laterally extensive genetic and stratigraphic units that may correspond to flow zones in groundwater aquifers and hydrocarbon reservoirs, typically several meters to tens of meters in thickness (e.g., Campbell 1967; Ringrose et al. 1993a, b; Jones et al. 1994, 1995; Koltermann and Gorelick 1996; Marsily et al. 1998; White and Barton 1999; Li and White 2003; Jackson et al. 2003; Deveugle et al. 2011). Understanding multiphase flow in layered porous media is therefore important for accurate prediction of many subsurface processes; examples include geologic carbon storage, migration of non-aqueous phase liquids (NAPLs) in contaminated aquifers and hydrocarbon production.

We examine immiscible, two-phase flow along homogeneous and parallel layers of contrasting petrophysical properties such as porosity and permeability. During the displacement of one phase by another, fluids may crossflow between adjacent layers, due to any combination of the viscous, capillary and gravitational forces which drive multiphase flow. This work investigates crossflow caused by viscous forces, which is commonly termed viscous crossflow (see Fig. 1a). We do not attempt to predict exactly the behavior of a given geologic reservoir or multiphase flow system; rather, our objectives are to predict where and how viscous crossflow affects the displacement of one fluid phase by another, as a function of a small number of key dimensionless numbers, in order to support mechanistic interpretations of more complex numerical model predictions, regardless of length scale (e.g., King and Mansfield 1999).

Previous laboratory experiments which examined the displacement of one fluid phase by another along parallel layers of contrasting grain size (e.g., Bertin et al. 1990; Dawe et al. 1992; Cinar et al. 2006; Alhamdan et al. 2012; Datta and Weitz 2013) include capillary effects, making the results inapplicable at larger scales at which viscous forces may be more important than capillary forces. Theoretical models and numerical simulations of flow that omit the contribution of capillary and buoyancy effects are thus required to examine the individual effect of viscous forces. Previous studies using theoretical models have investigated only flow along layers of contrasting absolute permeability, assuming layer properties are otherwise identical (Zapata and Lake 1981; Yortsos 1995), and piston-like displacements, in which the entire saturation change occurs at the shock-front (Zapata and Lake 1981 see Fig. 1b). However, in many geologic systems, layers are associated with contrasting porosity and relative permeability, in addition to the permeability contrasts investigated previously. Moreover, many displacement processes are not piston-like. Here we report a more general treatment that includes (i) non-piston-like displacements, in which some of the saturation change occurs over the shock-front, and some occurs over a rarefaction wave following the shock-front (see Fig. 1c), and (ii) layers of contrasting porosity and relative permeability, in addition to contrasting absolute permeability.

Fig. 1

a Schematic distribution of regions contacted by the displacing phase (represented in gray) during injection along a two-layered porous medium without significant capillary and buoyancy effects. As discussed later in this paper, the contrasting petrophysical properties of the layers may cause the shock-fronts in each layer to move with different velocities $U_1$ and $U_2$, and crossflow driven by viscous forces to occur between layers (represented as transverse arrows between the layers on the figure). Typical injected-phase saturation profiles along each layer (in the absence of crossflow) are reported for b piston-like displacements, in which the entire saturation change occurs at the shock-front, and c non-piston-like displacements, in which some of the saturation change occurs over the shock-front, and some occurs over a rarefaction wave following the shock-front

After first presenting the mathematical model (Sect. 2), we then present in Sect. 3 a comprehensive set of five dimensionless numbers characterizing immiscible, two-phase flow along layered porous media. We show in Sect. 4.1 that contrasts in porosity and relative permeability affect the shock-front propagation rates, in addition to previously examined contrasts in absolute permeability, and we identify flow regimes in which a fast shock-front and a slow shock-front can be defined. We rationalize crossflow behaviors for displacements which are not piston-like in Sect. 4.2, using a mobility ratio evaluated across the shock-front in each layer and find an additional crossflow regime, which has not been reported previously. We also present, in Sect. 4.3, complex crossflow patterns that have not been observed in previous studies. We report in Sect. 4.4 the sensitivity of storage efficiency, defined as the fraction of the moveable pore volume (MPV) occupied by the injected phase, on the governing dimensionless groups. These results directly support the interpretation of complex numerical subsurface models used to predict the amount of $\hbox {CO}_{2}$ that can be stored in the subsurface (e.g., Cavanagh and Ringrose 2011). The storage efficiency is numerically equivalent to the recovery efficiency, defined as the fraction of the MPV of the displaced phase that leaves the outflow face of the model, so the results are also applicable to hydrocarbon production (e.g., Christie and Blunt 2001).

2 Mathematical Model

We investigate two-phase, immiscible and isothermal flow through a two-layered porous medium in which the layers have contrasting petrophysical properties (Fig. 2). The model is a two-dimensional (2-D) symmetry element of an n-layered system in which alternating layers have the same contrasting properties. Assuming that the fluids and pore space are incompressible and that the pore space is completely filled with both fluids, and neglecting gravity and capillary forces, flow is described by the continuity equations

$$\begin{aligned}&\displaystyle \phi {\Delta }S\frac{\partial s}{\partial t}+\nabla \cdot {\varvec{q}}_{\varvec{i}} =0, \end{aligned}$$

(1)

$$\begin{aligned}&\displaystyle \nabla \cdot {\varvec{q}}_{\varvec{T}} =\nabla \cdot \left[ {{\varvec{q}}_{\varvec{i}} +{\varvec{q}}_{\varvec{d}} } \right] =0, \end{aligned}$$

(2)

and the multiphase Darcy’s law

$$\begin{aligned} {\varvec{q}}_{\varvec{{i}}} =-\frac{k_{r,i} \left( s \right) }{\mu _i }\mathop k\limits ^{=} \cdot {\varvec{\nabla }} {\varvec{P}} \end{aligned}$$

(3)

$$\begin{aligned} {\varvec{q}}_{\varvec{d}} =-\frac{k_{r,d} \left( s \right) }{\mu _d }\mathop k\limits ^{=} \cdot {\varvec{\nabla }} {\varvec{P}} \end{aligned}$$

(4)

where

$$\begin{aligned} s=\frac{S_i -S_{i,r} }{1-S_{i,r} -S_{d,r} } \end{aligned}$$

(5)

is the normalized injected phase saturation, which varies between 0 and 1, $S_i$ is the injected phase saturation, $S_{i,r}$ and $S_{d,r}$ are the injected and displaced phase residual saturations, ${\Delta }S=1-S_{i,r} -S_{d,r}$ is the moveable saturation, $\phi$ is the porosity, ${\varvec{q}}_{\varvec{i}}$ and ${\varvec{q}}_{\varvec{d}}$ are the injected and displaced phase volumetric fluxes per unit area, ${\varvec{q}}_{\varvec{T}} ={\varvec{q}}_{\varvec{i}} +{\varvec{q}}_{\varvec{d}}$ is the total volumetric fluid flux per unit area, $k_{r,i} (s)$ and $k_{r,d} \left( s \right)$ are the relative permeabilities of the injected and displaced phases, $\mu _i$ and $\mu _d$ are the viscosities of the injected and displaced phase, $\mathop k\limits ^{=}$ the absolute permeability tensor (assumed diagonal here) and P the fluid pressure. We further define the total fluid mobility $\lambda _T$ as

$$\begin{aligned} \lambda _T =\frac{k_{r,i} }{\mu _i }+\frac{k_{r,d} }{\mu _d}. \end{aligned}$$

(6)

The relative permeability curves are represented as functions of the normalized saturation s by the parametric forms

$$\begin{aligned}&\displaystyle k_{r,i} \left( s \right) =k_{r,i}^e s^{n_i }, \end{aligned}$$

(7)

$$\begin{aligned}&\displaystyle k_{r,d} \left( s \right) =k_{r,d}^e \left( {1-s} \right) ^{n_d }, \end{aligned}$$

(8)

where $k_{r,i}^e$ and $k_{r,d}^e$ are the end-point relative permeabilities, and $n_i$ and $n_d$ the Corey exponents of the injected and displaced phases, respectively. Each layer is internally homogeneous with identical length L and thickness H / 2 but, in contrast to previous studies, may differ in longitudinal absolute permeability $k_x$, porosity $\phi$, moveable saturation ${\Delta }S$, end-point relative permeabilities $k_{r,i}^e$ and $k_{r,d}^e$, and Corey exponents $n_i$ and $n_d$. We assume the layers have identical transverse permeability $k_z$ such that $k_z^1 =k_z^2 \le \hbox {min}\left( {k_x^1 ,k_x^2 } \right)$, where $k_z^1$ and $k_z^2$ are the transverse permeabilities of layers 1 and 2, respectively. Numerical experiments, not reported here, show that relaxing this constraint, to allow $k_z^1 \ne k_z^2$, but enforcing $k_z \le k_x$ in each layer as typically observed in geologic systems, has negligible impact on the results because crossflow is dominantly controlled by the lowest transverse permeability of the two layers.

Fig. 2

Schematic diagram of the two-layer model used in this work

Initially ($t=0)$, the pressure is uniform ($P=P_0 )$ and the normalized saturation is zero throughout the domain ($s=0)$. At the inlet face, the boundary conditions are a constant average volumetric influx per unit area $\bar{q}_{in}$, distributed to maintain a uniform but time-varying pressure ($P=P_{in} \left( t \right) )$. The other boundary conditions are a fixed pressure, equal to the initial pressure $P_0$, on the outlet (opposing) face, and no flow across the other faces. This choice of boundary conditions is consistent with conventional and physically reasonable borehole boundary conditions used in groundwater and oil reservoir models (Aziz and Settari 1979; Wu 2000).

3 Scaling Analysis

We report here a comprehensive set of five dimensionless numbers characterizing immiscible, two-phase flow along two-layered media. In the following treatment, we use the dimensionless distances $\hat{x}$ and $\hat{z}$, the dimensionless time $\hat{t}$ (which is equivalent to the number of moveable pore volumes of displacing phase injected), the dimensionless longitudinal and transverse volumetric fluid fluxes per unit area $\hat{q}_x$ and $\hat{q}_z$, the dimensionless shock-front velocity in the longitudinal direction $\hat{U}$, the normalized total mobility $\hat{\lambda }_T^j$ in layer $j=1,2$, and the dimensionless pressure $\hat{P}$, which are defined as

$$\begin{aligned}&\displaystyle {\hat{x}}=\frac{x}{L} \end{aligned}$$

(9)

$$\begin{aligned}&\displaystyle {\hat{z}}=\frac{z}{H} \end{aligned}$$

(10)

$$\begin{aligned}&\displaystyle \hat{t}=t\frac{\bar{q}_{in}}{\overline{\phi {\Delta }S}L} \end{aligned}$$

(11)

$$\begin{aligned}&\displaystyle {\hat{q}}_x =\frac{q_x }{\bar{q}_{in}} \end{aligned}$$

(12)

$$\begin{aligned}&\displaystyle \hat{q}_z =\frac{q_z }{\bar{q}_{in}}\frac{L}{H}\frac{\overline{k_z k_{r,d}^e }}{\overline{k_x k_{r,d}^e }} \end{aligned}$$

(13)

$$\begin{aligned}&\displaystyle \hat{U}=U\frac{\overline{\phi {\Delta }S}}{\bar{q}_{in} } \end{aligned}$$

(14)

$$\begin{aligned}&\displaystyle \hat{\lambda }_T^j =\frac{\lambda _T^j }{k_{r,d}^{e,j} /\mu _d },\quad j=1,2 \end{aligned}$$

(15)

$$\begin{aligned}&\displaystyle \hat{P}=\frac{\overline{k_x k_{r,d}^e }}{L{\bar{q}}_{in} \mu _d }\left( {P-P_0 } \right) \end{aligned}$$

(16)

where the barred quantities correspond to arithmetic averages, or the harmonic average in the case of transverse permeability. Note we apply the same scaling for the volumetric fluid fluxes per unit area to the total fluid flux ${\varvec{q}}_{\varvec{T}}$ and the injected phase fluid flux ${\varvec{q}}_{\varvec{i}}$.

We obtain a comprehensive set of five dimensionless numbers characterizing and providing rapid insights on flow by (i) identifying the required numbers from a dimensionless form of the flow equations and (ii) modifying the numbers obtained into an equivalent, more informative set of numbers directly applicable to characterize flow features such as the shock-front velocity ratio and crossflow behavior. The methodology used to derive the dimensionless numbers is reported in more detail in Appendix 1. The dimensionless numbers are summarized in Table 1 and briefly introduced below.

The effective aspect ratio is given by

$$\begin{aligned} R_L =\frac{L}{H}\sqrt{\frac{\overline{k_z k_{r,d}^e }}{\overline{k_x k_{r,d}^e }}} \end{aligned}$$

(17)

and quantifies the amount of crossflow between layers. When the effective aspect ratio is small ($R_L \ll 1)$, for example if the layers are separated by an impermeable barrier, then the layers are non-communicating (Lake 1989). Conversely, in the limit of large effective aspect ratios ($R_L \gg 1)$, the layers are perfectly communicating and crossflow occurs instantaneously such that transverse pressure gradients are negligible compared with longitudinal pressure gradients (Yortsos 1995). We refer to this limit as equilibrium crossflow; flow may be then treated as layer parallel (Zapata and Lake 1981; Lake 1989; Yortsos 1995). The ratio presented in Eq. (17) is similar to that presented previously by Zapata and Lake (1981), but here we account for contrasts in the relative permeability end-points between layers. The average of the longitudinal permeability to the injected phase is arithmetic while it is harmonic in the transverse direction.

The longitudinal permeability ratio is given by

$$\begin{aligned} \sigma _x =\frac{\left[ {k_x k_{r,d}^e } \right] _1 }{\left[ {k_x k_{r,d}^e } \right] _2 }, \end{aligned}$$

(18)

where the subscripts 1 and 2 denote the two layers, and scales the ratio of the influxes into each layer. Equation (18) is similar to the expression suggested by Goddin et al. (1966) except that we again account here for contrasts in the relative permeability end-points between layers. The storage ratio is given by

$$\begin{aligned} R_s =\frac{\left[ {\phi {\Delta }Ss_\mathrm{av} } \right] _1 }{\left[ {\phi {\Delta }Ss_\mathrm{av} } \right] _2 } \end{aligned}$$

(19)

and scales the ratio of moveable fluid volumes (accounting for the contribution of two-phase flow effects to the reduction of the moveable pore volume) in each layer. The storage ratio is new to this study as we include, for the first time, contrasts in porosity between layers, and also contrasts in the relative permeability curves which dictate the mobile saturation ${\Delta }S$ and the average normalized saturation behind the shock-front $s_\mathrm{av}$ that would be obtained without layer property contrasts (defined in Appendix 1).

The final two dimensionless numbers are the shock-front mobility ratio in each of the two layers 1 and 2. The shock-front mobility ratio describes the ratio of the total mobility across the shock-front (i.e., the total mobility calculated at the saturation values that bound the discontinuity that defines the shock) and is given by

$$\begin{aligned} M_f^1 =\frac{\lambda _T^1 \left( {s_f^1 } \right) }{\lambda _T^1 \left( {s_\infty ^1 } \right) } \end{aligned}$$

(20a)

$$\begin{aligned} M_f^2 =\frac{\lambda _T^2 \left( {s_f^2 } \right) }{\lambda _T^2 \left( {s_\infty ^2 } \right) } \end{aligned}$$

(20b)

in layers 1 and 2, respectively. Here, $s_f$ denotes the upper bound of saturation across the shock-front, and $s_\infty$ denotes the lower bound. We show later in Sect. 4 that the shock-front mobility ratios in each layer control the viscous crossflow behavior. Previous studies have used a single mobility ratio to characterize immiscible flow in layered porous media (e.g., Zapata and Lake 1981); here, we define a mobility ratio for each layer to account for relative permeability contrasts between layers.

In the immiscible two-phase flow studied here, the dimensionless numbers allow extrapolation of results from one system to another if the systems have relative permeability curves with identical shapes (Rapoport 1955), here expressed by the Corey exponents. In the case of piston-like displacements, which were studied previously (e.g., Zapata and Lake 1981), the relative permeability curves only influence the displacement via their end-point values because the end-point and shock-front mobility ratios are identical, so it is only required that the two systems have identical end-point mobility ratios. In this work, we consider displacements which are not piston-like, so the relative permeabilities of the two phases vary behind the shock-front, making exact scaling with two dimensionless numbers inaccessible. However, we show later, via asymptotic flow solutions, that the dimensionless numbers presented above (and summarized in Table 1) are sufficient to achieve an approximate scaling of flow suitable for our purpose.

Table 1 Governing dimensionless numbers and range of values explored

The effective aspect ratio $R_L$ typically varies over the range 0.1–100 in layered sedimentary systems (L / H is typically large, of order 10–100, while $k_z /k_x$ is typically small, of order 10$^{-4}$ – 1; see Table 1). The longitudinal permeability and storage ratio can vary over several orders of magnitude, due to variations in grain size ($\sigma _x$ and $R_s \gg 1)$ or sorting ($\sigma _x \gg 1$ and $R_s \ll 1)$. Several studies have explored correlations between porosity, absolute permeability, and relative permeability, which could be used to restrict the possible combinations of permeability ratio and storage ratio (e.g., Nelson 1994). However, measured data often show a poor correlation (e.g., Thompson et al. 1987) and here, for the sake of generality, we do not restrict the combinations of permeability ratio and storage ratio investigated, varying the permeability ratio over the range 1–10,000 and the storage ratio over the range 0.01–100. A suite of core-plug measurements taken along a single well from a North Sea field (Tjølsen et al. 1991) shows that plausible combinations of the permeability and storage ratios span one or two orders of magnitude (Fig. 3a). As we show later, there is a consistent change in crossflow behavior at a shock-front mobility ratio $M_f =1$, and we investigate values that span this threshold, ranging from 0.6 to 1.4. Plausible combinations of the fast and slow shock-front mobility ratios are shown in Fig. 3b.

Fig. 3

Possible combinations of a permeability/storage ratios and b shock-front mobility ratios calculated from core-plug measurements taken along a single well from a North Sea field (data from Tjølsen et al. 1991). Combinations of the shock-front mobility ratios are calculated using the relative permeability data from Tjølsen et al. (1991), assuming viscosity ratios $\mu _{_d} /\mu _i =1$, 10 and 100

We quantify the impact of the dimensionless numbers on flow characterized in terms of a dimensionless storage efficiency, defined as

$$\begin{aligned} E_s =\frac{\int \int s\phi {\Delta }S\hbox {d}x\hbox {d}z}{\int \int \phi {\Delta }S\hbox {d}x\hbox {d}z}. \end{aligned}$$

(21)

The storage efficiency measures how effectively the injected phase is retained within the model and is relevant when characterizing the geologic storage of carbon in subsurface reservoirs and the location of NAPLs in contaminated aquifers. The storage efficiency is also numerically equivalent to the recovery efficiency, which is a measure of how effectively the displaced phase is removed from the model and is relevant to hydrocarbon production. Quantifying the effect of the dimensionless numbers in terms of the storage/recovery efficiency (henceforth termed the storage efficiency) therefore yields results of broad interest. To further analyze the controls on storage efficiency, we also decompose storage efficiency following Lake (1989) into the product of sweep efficiency, defined as

$$\begin{aligned} E_{sw} =\frac{\int \int \phi {\Delta }S1_{s>0} \hbox {d}x\hbox {d}z}{\int \int \phi {\Delta }S\hbox {d}x\hbox {d}z}, \end{aligned}$$

(22)

which measures the fraction of the moveable pore volume contacted by the injectant, and displacement efficiency, defined as

$$\begin{aligned} E_d =\frac{\int \int s\phi {\Delta }S\;\hbox {d}x\hbox {d}z}{\int \int \phi {\Delta }S1_{s>0}\; \hbox {d}x\hbox {d}z}, \end{aligned}$$

(23)

which quantifies the fraction of the contacted moveable pore volume which has been displaced by the injected phase.

4 Results

4.1 Impact of Layer Property Contrasts on the Shock-Front Velocity

Material property contrasts in layered systems cause the shock-front to propagate with a different velocity in each layer. This is important, because contrasts in the shock-front velocity impact on viscous crossflow and also on the storage efficiency. For layers that only differ in their longitudinal permeabilities, it is well known that the shock-front propagates faster through the high permeability layers (Dykstra and Parsons 1950; Zapata and Lake 1981). However, the behavior is more complex when the layers also differ in their porosities and relative permeabilities. We analyze the impact of these material property contrasts by quantifying the shock-front velocities in the limits of no crossflow ($R_L \ll 1$; see Appendix 2) and equilibrium crossflow ($R_L \gg 1$; see Appendix 3). Figure 4a, b summarizes the flow regions in which a fast shock-front and a slow shock-front can be defined as a function of the permeability ratio ($\sigma _x$) and the storage ratio ($R_s )$. When the ratio $\sigma _x /R_s$ is sufficiently larger or smaller than one, it is possible to define without ambiguity a fast shock-front and a slow shock-front. When the ratio $\sigma _x /R_s$ is close to one, it is not necessarily possible to define a fast shock-front and a slow shock-front. In the equilibrium crossflow limit, the two shock-fronts move at equal velocities. In the no-crossflow limit, an initially slower shock-front may become faster than the initially faster shock-front because the total mobility changes differently in each layer (see Appendix 2 for further details).

Fig. 4

Ratio of shock-front velocities in layer 1 and 2 as function of the dimensionless groups in the limits of a no crossflow ($R_L =0$) and b equilibrium crossflow ($R_L \gg 1$). In the no-crossflow limit (a), fast and slow shock-fronts cannot be defined without ambiguity in the region between the two dashed lines. In this part of the parameter space, the relative shock-front velocities are sensitive to the specific shapes of the relative permeability curves

Fig. 5

Schematic of pressure profiles obtained in the fast (F) and slow (S) layers without crossflow ($R_L \ll 1$) for various shock-front mobility ratios $M_f^F$ and $M_f^S$. The pressure profiles are used to predict the displacing-phase crossflow shown in the adjacent schematics. Arrows indicate the displacing-phase crossflow direction across the interface separating the two layers; crosses indicate no displacing-phase crossflow

4.2 Prediction of Crossflow Regimes

Having identified fast and slow shock-fronts in the previous section, we now consider the crossflow regimes, predicting these by comparing the pressure profiles obtained in each layer with negligible crossflow ($R_L \ll 1)$. Displacing-phase crossflow is predicted to occur from high to low pressure regions, although it is important to note that the displacing phase can only crossflow out of the invaded part of a layer behind the shock-front. Here we consider displacements which are not piston-like, so pressure gradients may vary behind the shock-front due to spatial changes in total fluid mobility. We construct pressure profiles for such cases by assuming that the total mobility behind the shock-front may be approximated by the total mobility at the shock-front

$$\begin{aligned} \hat{\lambda }_T \left( s \right) \approx \hat{\lambda }_T \left( {s_f } \right) \end{aligned}$$

(24)

Under this approximation, the pressure profile along a given layer is parameterized by the shock-front mobility ratio $M_f$: pressure decreases linearly on each side of the shock-front, and the ratio of pressure gradients ahead of and behind the shock-front is equal to $M_f$ (see Appendix 2). For the sake of simplicity, we assume one of the shock-fronts consistently propagates faster than the other in the limit $R_L \ll 1$, and we refer to this shock-front as the fast (F) shock-front as opposed to the slow (S) shock-front. This is reasonable for $\sigma _x /R_s$ ratios that are sufficiently larger or smaller than one (Fig. 4a). Pressure profiles are expressed with respect to the fast and slow shock-front mobility ratios, $M_f^F$ and $M_f^S$, and are summarized in Fig. 5. We can interpret these to predict injectant crossflow by comparing the pressure and the location of the shock-front in each layer. We observe three distinct crossflow regimes. For $M_f^F =M_f^S <1$ and $M_f^F>1>M_f^S$ (the lower two quadrants in Fig. 5), crossflow of the injected phase occurs from the fast to the slow layer. For $M_f^F<1<M_f^S$ (the upper left quadrant in Fig. 5), crossflow occurs in the opposite direction, from the slow to the fast layer. This crossflow regime, introduced by the contrasts in relative permeability, has not been reported previously. For $M_f^F =M_f^S >1$ (the upper right quadrant in Fig. 5), significant crossflow occurs in both directions: from the slow layer to the fast layer behind the slow shock-front and from the fast to the slow layer ahead of the slow shock-front. As we show in the next section, this latter crossflow pattern gives rise to complex changes in saturation that have not been reported previously. Note that the crossflow regimes shown in Fig. 5 apply irrespective of the value of the permeability ($\sigma _x$) and storage ($R_s$) ratios, which only define the fast and slow layers (Fig. 4a). For piston-like displacements without relative permeability contrasts, the fast and slow shock-front mobility ratios are equal to the end-point mobility ratio $M_e =k_{r,i} \mu _d /k_{r,d} \mu _i$, and we recover the two crossflow regimes identified by Zapata and Lake (1981) for $M_e <1$ and $M_e >1$.

Fig. 6

Injectant crossflow behaviors as function of $M_f^F$ and $M_f^S$. The fast and slow shock-front mobility ratios are alternatively chosen as 0.6 and 1.2, with $R_L =1$, $\sigma _x =5$ and $R_s =1$. Dashed lines represent the shock-front positions obtained without crossflow

4.3 Impact of Moderate Crossflow on Non-piston-Like Displacements

We now examine crossflow for a moderate effective aspect ratio ($R_L =1$) to confirm analytical predictions from the previous section. Solutions of the multiphase flow Eqs. (1)–(4) were obtained using a commercial code that implements a finite-volume-finite-difference approach to discretize the governing equations (Eclipse 100). Flow was simulated using two-dimensional Cartesian grids with resolutions ranging from $200 \times 200$ up to $800 \times 800$ cells to demonstrate that the solutions were converged. Time stepping was fully implicit. The selected Corey exponents are $n_i =4$, $n_d =2$ in layer 1 and $n_i =2$, $n_d =3$ in layer 2; these values are typical for geologic porous media and yield non-piston-like displacements (e.g., Fig. 1c). The permeability ratio and the storage ratio are chosen as $\sigma _x =5$ and $R_s =1$. The shock-front travels faster in layer 1 ($\sigma _x /R_s =5)$, so we henceforth refer to layer 1 as the fast layer and layer 2 as the slow layer. Figure 6 shows the injected phase saturation obtained for values of the shock-front mobility ratios chosen to be either 0.6 or 1.2. Crossflow directions are obtained through examination of the transverse fluxes of the injected phase, measured along the interface between the two layers, and are reported as white arrows on the saturation maps. These numerical experiments confirm the crossflow directions predicted analytically in Sect. 4.2. We also observe the formation of viscous fingers in the fast layer (see plots on right of Fig. 6), which is expected when the mobility ratio across the shock-front is larger than one (Riaz and Tchelepi 2006). Viscous fingers in the slow layer were also observed at later times, when the mobility ratio across the slow shock-front is larger than one. Comparison of shock-front positions with and without crossflow (see dashed lines in Fig. 6) shows crossflow also influences the volumes of displacing phase injected into each layer, as crossflow influences the total fluid mobility in each layer. We find that viscous crossflow only reduces the distance between the two shock-fronts when $M_f^F <1$. The impact of crossflow on storage efficiency will be further explored in Sect. 4.4.

In most cases, crossflow occurs only from the fast layer to the slow layer or vice-versa; however, for $M_f^F =M_f^S >1$, crossflow occurs from the slow layer to the fast layer behind the slow shock-front and from the fast to the slow layer ahead of the slow shock-front (Figs. 6, 7). Similar crossflow behaviors were also reported by Zapata and Lake (1981). A point of zero crossflow can be defined between these two opposing crossflow directions, which travels along the interface between the two layers toward the outflow face (Fig. 7b). The rotational nature of the crossflow around this point leads to the development of complex saturation patterns that have not been observed previously (Fig. 7a). These complex patterns were also obtained without relative permeability contrasts (Corey exponents were chosen as $n_i =n_d =2$ in both layers), confirming that the crossflow regimes identified here are not specific to our initial choice of relative permeability curves or to the relative permeability contrast. To eliminate the possibility of the complex saturation pattern representing a numerical artifact, a grid sensitivity test was conducted with grid resolution varying from 100 x 100 up to 800 x 800 cells. Numerical solutions were also obtained using the control-volume finite-element method of Jackson et al. (2015). These sensitivity tests yielded the same crossflow behavior regardless of grid resolution or numerical method.

Fig. 7

a Injected phase saturation at $t=1$ MPVI for $M_f^F =M_f^S =1.4$, $R_L =1$, $\sigma _x =5$, $R_s =1$, and b dimensionless transverse displacing phase volumetric flux per unit area $\hat{q}_{i,z}$ across the interface between the two layers corresponding to (a). Positive fluxes represent fluxes from the slow layer to the fast layer

4.4 Storage/Recovery Efficiency as Function of the Dimensionless Numbers

We continue our analysis by exploring the storage/recovery efficiency as a function of the permeability ratio, the storage ratio, the shock-front mobility ratios and the effective aspect ratio (Table 1). The correlations reported here provide a useful basis for the interpretation of more complex numerical subsurface models used to quantitatively predict storage or recovery efficiency.

Fig. 8

Storage efficiency versus the permeability ratio and the storage ratio for a unit mobility ratio piston-like displacement

Variations in $\sigma _x$ and $R_s .$ To investigate the relationship between $\sigma _x$, $R_s$ and storage efficiency, we initially consider the special case of a piston-like displacement, with the additional constraint that the injected and displaced phases have equal mobilities, so there are no fluid mobility contrasts and $M_f^F =M_f^S =1$. This is ensured by choosing suitable values of fluid viscosity and end-point relative permeability. For this case, the shock-front velocities are constant in each layer (but vary between layers), the shock-front velocity ratio being equal to $\sigma _x /R_s$ (apply $\overline{\lambda _T }^{1}\left( {\hat{t}} \right) =\overline{\lambda _T }^{2}\left( {\hat{t}} \right) =1$ in Eq. (37)). Solutions for this type of displacement are independent of the value of the effective aspect ratio $R_L$ as there is no crossflow, and are also independent of the Corey exponents (i.e., the shape of the relative permeability curves) so long as these are chosen to yield a piston-like displacement. Given this, and the fixed value of $M_f$ in each layer, flow is described in terms of just two dimensionless parameters ($\sigma _x$ and $R_s )$.

These assumptions allow the storage efficiency to be expressed as a function of $\sigma _x$ and $R_s$ (see Appendix 4). They also allow storage efficiency to be described in terms of sweep efficiency (defined in Eq. 22), these being equivalent for piston-like displacements. Figure 8 plots the storage efficiency contours obtained at 1 MPVI. We find that storage efficiency is maximum when the shock-fronts travel at identical velocities, i.e., $\sigma _x /R_s =1$, and decreases with increasing differences in the shock-front velocity, i.e., when $\sigma _x /R_s$ increases or decreases away from one. While shock-front velocity differences indicate a reduced fraction of the rock volume contacted by the displacing phase, additional knowledge of the storage ratio is required to quantify their impact on the fraction of the moveable pore volume contacted by the displacing phase which, for piston-like displacements, corresponds to the sweep and storage efficiency. This additional dependency on the storage ratio explains the observed, nonlinear relationship between storage efficiency, the permeability ratio and the storage ratio. Previous studies concluded that the maximum storage efficiency (often expressed previously as the equivalent recovery efficiency) is obtained for unit permeability ratio ($\sigma _x =1)$ (Dykstra and Parsons 1950); we see here that this is true only if there are no storage contrasts. Moreover, the permeability ratio only has a significant impact on storage efficiency for values less than $\sigma _x \approx 10$; at higher permeability ratio, the storage efficiency is approximately constant irrespective of permeability ratio. In contrast, the storage ratio impacts on storage efficiency over most of the range investigated, and becomes more significant as the permeability ratio increases above 10. These results suggest that is important to consider contrasts in porosity and relative permeability, in addition to contrasts in permeability, when assessing storage efficiency in layered systems. They provide a benchmark against which we can compare non-piston-like displacements with varying fluid mobility, which are considered next.

We now consider non-piston-like displacements with fixed and indicative values of $R_L =1$ and $M_f^F =M_f^S =0.6$, so moderate amounts of crossflow occur. Solutions reported below (Fig. 9) and in the rest of the section are obtained via numerical simulations, as described in the previous section. Figure 9 shows that storage efficiency follows the trend observed in Fig. 8, which was obtained assuming piston-like displacements without crossflow. Additional numerical experiments, not reported here, confirm these trends are maintained regardless of the values of the other dimensionless groups. Without the piston-like assumption, the displacement efficiency (defined in Eq. 23) can also vary with the dimensionless parameters, in addition to the sweep efficiency. However, we find here that the displacement efficiency only weakly depends on the permeability ratio and the storage ratio.

Fig. 9

Storage efficiency at t=1 MPVI versus a, b the permeability ratio (with $R_s =0.1$ and 10, respectively) and c the storage ratio (with $\sigma _x =2$ ), for $R_L =1$, $M_f^1 =M_f^2 =0.6$

Variations in $M_f^1\,$ and $M_f^2$. We now vary $M_f^1$ and $M_f^2$ with fixed and indicative values of $R_L =1$, $\sigma _x =5$ and $R_s =1$ to yield moderate crossflow and shock-front velocity ratio, finding that the storage efficiency decreases as both the fast and slow shock-front mobility ratios increase (Fig. 10a). Additional numerical experiments, not reported here, confirm these trends are maintained regardless of the values of the other dimensionless groups. The reduction in storage efficiency is explained by the reduced displacement efficiency (defined in Eq. 23) with which the phase initially occupying the pore space is displaced by the injected phase. It is well known that displacement efficiency in 1D flow, quantified by the average saturation behind the shock-front (see also Appendix 1), decreases with increasing shock-front mobility ratio (Fig. 10b).

Fig. 10

a Storage efficiency at $t=1$ MPVI versus $M_f^1$ and $M_f^2$ (with $R_L =1$, $\sigma _x =5$ and $R_s =1$ ), and b average saturation behind the shock-front assuming unidimensional flow versus the shock-front mobility ratio for two sets of Corey exponents $\left( {n_i ,n_d } \right)$

Variations in $R_L$. We finish by varying $R_L$ to explore the impact of crossflow on storage efficiency for the various crossflow regimes presented in Fig. 6. Without loss of generality, we choose the permeability ratio and the storage ratio ($\sigma _x =5$ and $R_s =1$) such that crossflow occurs. The layer in which the shock-front travels faster is referred to as the fast layer. We find that the effect of changing the effective aspect ratio $R_L$ on storage efficiency depends on the dimensionless time (number of moveable pore volumes injected) at which efficiency is measured, and on the fast shock-front mobility ratio (Fig. 11). Regardless of the mobility ratios, storage efficiency at breakthrough (defined as the time at which the displacing phase reaches the outlet face of the model) increases with $R_L$ (Fig. 11).

At early times post-breakthrough, increased crossflow (i.e., increased $R_L )$ yields increased storage efficiency for $M_f^F <1$, but decreased storage efficiency for $M_f^F >1$ (see, for example, changes in storage efficiency with $R_L$ at $t=1$ MPVI in the left quadrants of Fig. 11, and at $t=3.2$ and 8 MPVI in the right quadrants of Fig. 11). At late times, storage efficiency becomes more weakly dependent on $R_L$ (see, for example, the negligible changes in storage efficiency at $t=20$ MPVI in Fig. 11). These results are consistent with earlier findings obtained by Zapata and Lake (1981) for piston-like displacements. However, Zapata and Lake (1981) found the influence of $R_L$ on storage efficiency at breakthrough time (note that they did not report their results in terms of storage efficiency, but they can be expressed in this way) to be more significant than reported here. This may be explained by their use of a unidimensional flow model to calculate storage efficiency in the equilibrium crossflow limit, which overestimates the impact of crossflow. Numerical solutions reported in Appendix 3 (Fig. 15) show the injected phase flowing ahead of the slow shock-front position predicted using an equilibrium crossflow model similar to the one used by Zapata and Lake (1981). This results in numerical solutions yielding earlier breakthrough than unidimensional flow models for large effective aspect ratios $R_L$, thereby predicting smaller changes in storage efficiency with $R_L$ than previously observed.

Fig. 11

Change in storage efficiency compared to efficiency obtained without crossflow ($R_L =0.01$) for various shock-front mobility ratios. Results are shown at breakthrough, at early times post-breakthrough ($t=1$, 3.2 or 8 MPVI) and late time ($t=20$ MPVI). In all cases, $\sigma _x =5$ and $R_s =1$

5 Discussion

The results reported here are applicable to immiscible, incompressible flow in layered porous media irrespective of material property contrasts, fluid property contrasts and length scale, so long as capillary and buoyancy effects are negligible. Previous scaling analysis has shown that the viscous limit explored here is typically observed at high flow rates, large length scales, low permeabilities, and with fluids having similar densities (e.g., Shook et al. 1992). Example applications for our model include plug-scale experiments in the laboratory (10’s cm scale), water flooding of hydrocarbon reservoirs (100’s m scale) and $\hbox {CO}_{2}$ storage in regional aquifers (km scale). The five dimensionless numbers (the two shock-front mobility ratios, the effective aspect ratio, the permeability ratio and the storage ratio) allow possible flow regimes to be assessed in the context of data-poor scenarios of geologic heterogeneity and the associated uncertainty in storage or recovery efficiency to be explored.

Our results show that crossflow is inevitable unless the effective aspect ratio is very small (i.e., for $R_L <0.1)$. In most geologic settings, the effective aspect ratio is large; layers tend to be long and thin (so L / H is large) and, although $k_v /k_h$ tends to be small, the effective aspect ratio depends on the square root of this property. Thus, taking $L/H\sim 10$ and $k_v /k_h \sim 10^{-4}$ as typical minimum values yields a minimum $R_L$ of $\sim$0.1. Viscous crossflow in layered geologic systems is therefore likely unless there is a continuous barrier to flow (such as a mudstone or cemented horizon) between layers.

Our results also show that knowledge of the end-point mobility ratio is not sufficient to evaluate crossflow. The shape of the relative permeability curves is also important, as these strongly influence the shock-front mobility ratio in each layer, and therefore the crossflow regime. For example, the end-point mobility ratio for $\hbox {CO}_{2}$ injection into deep saline aquifers is typically large: the $\hbox {CO}_{2}$-brine viscosity ratio is low at the pressure and temperature conditions relevant to geologic carbon storage (between 0.02 and 0.2), while the end-point relative permeability is larger for the (wetting) brine phase than for the (non-wetting) $\hbox {CO}_{2}$ phase. However, depending on the shape of the relative permeability curves, the shock-front mobility ratio may be smaller or larger than one (Fig. 12). Likewise, knowledge of the permeability ratio is not sufficient to predict crossflow. Both the permeability ratio, which captures the impact of longitudinal permeability contrasts, and the storage ratio, which captures the impact of porosity and end-point saturation contrasts, are important to assess storage efficiency.

Fig. 12

Shock-front versus end-point mobility ratio for drainage CO$_2$-brine relative permeability measurements at reservoir conditions from Bennion and Bachu (2005, 2006)

6 Conclusions

This work examined the effect of viscous forces on the displacement of one fluid by a second, immiscible fluid along parallel, continuous layers of contrasting porosity, permeability and relative permeability. The two-phase flow was characterized using five dimensionless numbers that are new to this study:

• 1. The permeability ratio $\sigma _x$, which captures contrasts in longitudinal (layer-parallel) permeability and relative permeability.

• 2. The storage ratio $R_s$, which captures contrasts in porosity and the end-point saturations of the relative permeability curves.

• 3. The effective aspect ratio $R_L$, which rescales the aspect ratio to account for anisotropic permeability and relative permeability.

• 4. Two shock-front mobility ratios defined in each layer, $M_f^1$, and $M_f^2$, which capture the fluid mobility contrast across the shock-front at the leading edge of the displacement in each layer.

The impact of the dimensionless numbers on flow was quantified in terms of a dimensionless storage efficiency, so results are directly applicable, regardless of scale, to geologic carbon storage. The storage efficiency is also numerically equivalent to the recovery efficiency, relevant to hydrocarbon production.

In addition to contrasts in longitudinal permeability, characterized by the permeability ratio ($\sigma _x )$, we show that contrasts in porosity and relative permeability, characterized by the storage ratio ($R_s )$, affect the velocity with which the shock-fronts move through each layer and the storage efficiency. The difference in shock-front velocities increases, and the storage efficiency decreases, as the ratio $\sigma _x /R_s$ deviates from one. When this ratio is one, the storage efficiency is maximal. Previous studies have concluded that the maximum storage efficiency (usually expressed as the equivalent recovery efficiency) in layered systems is obtained for unit permeability ratio ($\sigma _x =1)$; we show that this is true only if there are no storage contrasts ($R_s =1)$. Moreover, the permeability ratio only has a significant impact on storage efficiency when it is less than approximately 10; at higher permeability ratio, the storage efficiency is approximately constant irrespective of permeability ratio. In contrast, the storage ratio impacts on storage efficiency over most of the range investigated and becomes more significant as the permeability ratio increases above 10. It is therefore important to consider contrasts in porosity and relative permeability, in addition to contrasts in permeability, when assessing storage efficiency in layered systems.

When crossflow occurs, it is possible to use the shock-front velocity to identify a fast layer and a slow layer; the shock-front velocity is higher in the fast layer than the slow layer. In some cases, the shock-front velocities are identical in each layer even though the layers have contrasting porosity, absolute permeability and relative permeability. Crossflow can allow the displacing phase to move with a uniform shock-front through a heterogeneous, layered system; this is a highly counterintuitive result. Three crossflow regimes can be identified based on the shock-front mobility ratios in the fast ($M_f^F$) and slow ($M_f^S$ ) layers. When $M_f^F =M_f^S <1$ and $M_f^F>1>M_f^S$, then displacing-phase crossflow only occurs from the fast layer to the slow layer; conversely, when $M_f^F<1<M_f^S$, then displacing-phase crossflow only occurs from the slow layer to the fast layer. When $M_f^F =M_f^S >1$, we observe complex crossflow patterns that have not been reported previously, with displacing-phase crossflow from the slow layer to the fast layer behind the slow shock-front and from the fast layer to the slow layer ahead of the slow shock-front.

Regardless of the values of the other dimensionless numbers, storage efficiency decreases with increasing shock-front mobility ratio. The impact of crossflow, quantified by the effective aspect ratio ($R_L )$, on storage efficiency depends on the time at which the storage efficiency is measured. Regardless of the values of the other dimensionless numbers, storage efficiency at breakthrough increases with the effective aspect ratio. This increase is maintained at late times when $M_f^F <1$, but only until some finite post-breakthrough time when $M_f^F >1$.

Appendices

Appendix 1: Derivation of Governing Dimensionless Numbers

The set of dimensionless groups presented in Table 1, which is used to describe immiscible, two-phase flow along a two-layered porous medium, is obtained in two steps. A first set of dimensionless numbers required to comprehensively characterize flow is obtained from a dimensionless form of the flow equations, following the commonly termed ‘inspectional analysis’, which has been previously applied to homogeneous (Shook et al. 1992) and simple layered porous media (Zhou et al. 1997). Guided by asymptotic flow solutions discussed in this work, we replace some of the obtained numbers by equivalent numbers to provide deeper insight into key flow features such as the shock-front velocity ratio or the crossflow behavior.

Inspectional analysis of the flow equations. Before non-dimensionalizing the flow equations, we first express the governing flow Eqs. (1)–(4) within the so-called fractional flow formulation. This requires the injected phase volumetric flux per unit area to be expressed as a function of the total flux,

$$\begin{aligned} {\varvec{q}}_{\varvec{i}} =-\frac{k_{r,i} }{\mu _i }\mathop k\limits ^{=} \cdot {\varvec{\nabla }} {\varvec{P}}=\frac{\frac{k_{r,i} }{\mu _i }}{\lambda _T }\cdot \left[ {-\lambda _T \mathop k\limits ^{=} \cdot {\varvec{\nabla }} {\varvec{P}}} \right] =f_i q_T , \end{aligned}$$

(25)

where the dimensionless ratio of the injected to the total fluid mobility ($f_i$, which is commonly called ‘fractional flow’) controls the efficiency with which the injected phase displaces the phase initially in place (see a typical fractional flow curve and its impact on a two-phase displacement in Fig. 13).

Fig. 13

a Typical fractional flow curve as a function of saturation, and b corresponding saturation profile (commonly called ‘Buckley–Leverett profile’) during a displacement along a homogeneous porous medium. The shock-front saturation $s_f$ and average saturation behind the shock-front $s_\mathrm{av}$ obtained in a homogeneous porous medium can be determined graphically as shown on the figure from the construction of the Welge tangent (represented by the dashed line tangent to the fractional flow curve) (Welge 1952). The average saturation behind the shock-front also corresponds to the reciprocal of $\left[ {\mathbf{d}f_i /\mathbf{d}s} \right] \left( {s_f } \right)$, which quantifies the shock-front velocity increase related to the bypassing of the displaced phase by the displacing phase

Re-injecting the latter expression into the continuity equation (1), yields the ‘fractional flow’ (dimensional) formulation of the governing flow equations

$$\begin{aligned}&\displaystyle \phi {\Delta }S\frac{\partial s}{\partial t}+{\varvec{q}}_{\varvec{T}} \cdot {\varvec{\nabla }} {\varvec{f}}_{\varvec{i}} =0, \end{aligned}$$

(26)

$$\begin{aligned}&\displaystyle \nabla \cdot {\varvec{q}}_{\varvec{T}} =\nabla \cdot \left[ {-\lambda _{T}} {\mathop k\limits ^{=}} \cdot {\varvec{\nabla }} {\varvec{P}} \right] =0. \end{aligned}$$

(27)

Normalizing flow Eqs. (26), (27) using the dimensionless quantities defined in Eqs. (9)–(15) and the following dimensionless quantities, $\hat{k}_{x,d} =k_x k_{r,d}^e /\overline{k_x k_{r,d}^e }$, $\hat{k}_{z,d} =k_z k_{r,d}^e /\overline{k_z k_{r,d}^e }$ and $\hat{C}_s =\phi {\Delta }S/\overline{\phi {\Delta }S}$, we obtain the following dimensionless form governing flow equations

$$\begin{aligned}&\displaystyle \hat{C}_s \frac{\partial s}{\partial \hat{t}}+\hat{q}_{T,x} \frac{\partial f_i }{\partial \hat{x}}+R_L^2 \hat{q}_{T,z} \frac{\partial f_i }{\partial \hat{z}}=0, \end{aligned}$$

(28)

$$\begin{aligned}&\displaystyle \frac{\partial \hat{q}_{T,x} }{\partial \hat{x}}+R_L^2 \frac{\partial \hat{q}_{T,z} }{\partial \hat{z}}=\frac{\partial }{\partial \hat{x}}\left[ {-\hat{k}_{x,d} \hat{\lambda }_T \frac{\partial \hat{P}}{\partial \hat{x}}} \right] +R_L^2 \frac{\partial }{\partial \hat{z}}\left[ {-\hat{k}_{z,d} \hat{\lambda }_T \frac{\partial \hat{P}}{\partial \hat{z}}} \right] =0. \end{aligned}$$

(29)

From these dimensionless flow equations, we identify the governing dimensionless numbers as constant coefficients appearing in the partial differential equations or constant coefficients controlling functionals appearing in the differential equations. The effective aspect ratio $R_L$ directly appears as a constant coefficient within the coupled partial differential equations. The end-point mobility ratios $M_e =k_{r,i}^e \mu _d /k_{r,d}^e \mu _i$ and the Corey exponents $n_i$ and $n_d$, defined in each layer $j=1,2$, are identified as a consequence of our relative permeability parameterization from

$$\begin{aligned} f_i^j \left( s \right) =\frac{M_e^j s^{n_{i,j} }}{M_e^j s^{n_{i,j} }+\left( {1-s} \right) ^{n_{d,j} }}, \end{aligned}$$

(30)

and

$$\begin{aligned} \hat{\lambda }_T^j =M_e^j s^{n_{i,j} }+\left( {1-s} \right) ^{n_{d,j} }. \end{aligned}$$

(31)

The longitudinal and transverse permeability ratios, $\left( {k_x k_{r,d}^e } \right) _1\Big /\left( {k_x k_{r,d}^e } \right) _2$ and $\left( {k_z k_{r,d}^e } \right) _1\Big /\left( {k_z k_{r,d}^e } \right) _2$, as well as the ratio of moveable pore volumes in each layer, $\left( {\phi {\Delta }S} \right) _1 /\left( {\phi {\Delta }S} \right) _2$, are identified from the spatially dependent functionals $\hat{k}_{x,d}$, $\hat{k}_{z,d}$ and $\hat{C}_s$. Considering these functionals are constructed so their spatial averages (in the arithmetic sense for $\hat{k}_{x,d}$ and $\hat{C}_s$, but in the harmonic sense for $\hat{k}_{z,d} )$ is equal to one, the latter permeability ratios and storage ratio are sufficient to parameterize the functionals.

Contribution of asymptotic flow solutions to dimensionless number selection. Having identified the required numbers to comprehensively characterize flow, we choose to replace the end-point mobility ratios $M_e =k_{r,i}^e \mu _d /k_{r,d}^e \mu _i$ in each layer, and the ratio of moveable fluid volumes $\left( {\phi {\Delta }S} \right) _1 /\left( {\phi {\Delta }S} \right) _2$, by the shock-front mobility ratios in each layer defined in Eq. (20) and the storage ratio defined in Eq. (19). First, the shock-front mobility ratios in each layer, rather than the end-point mobility ratios, capture the leading-order crossflow behavior in layered systems, as discussed in Sect. 4.2. Because the shock-front mobility ratio is a strictly increasing function of the end-point mobility ratio for fixed Corey exponents, replacing the end-point mobility ratios by the shock-front mobility ratio does not affect the consistency of the proposed dimensionless groups, so we prefer to retain the latter as being more informative of flow behavior. Second, we replace the ratio of the moveable pore volumes $\left( {\phi {\Delta }S} \right) _1 /\left( {\phi {\Delta }S} \right) _2$ by the storage ratio defined in Eq. (19), which includes the average saturation $s_\mathrm{av}$ that would be obtained in each layer without layer property contrasts. This accounts for the contribution of two-phase flow effects to the reduction of the moveable pore volume and yields a compact mathematical expression of the shock-front velocity ratios calculated in the limit of no crossflow and equilibrium crossflow (see “Appendices 2 and 3”). Because the average saturation $s_\mathrm{av}$ defined in each layer is entirely determined by the end-point mobility ratio and the Corey exponents, this change does not affect the consistency of the proposed dimensionless groups.

Numerical validation of the dimensionless numbers. To validate the proposed dimensionless numbers, we present a numerical sensitivity test (using the numerical method described in Sect. 4.3) that varies several of the dimensional parameters but maintain the values of the dimensionless numbers constant (see Table 2). Storage efficiency histories and saturation distributions obtained were identical between the three cases (see storage efficiency histories reported in Table 3; observed differences are less than one-tenth of one percent). This confirms that the presented dimensionless results are reproducible independent of scale.

Table 2 Data used for the numerical sensitivity tests to validate the dimensionless numbers

Table 3 Dimensionless results for three cases described in Table 2

While the transverse permeability ratio $\sigma _z$ is formally required to scale the displacement, the ratio has no influence on storage efficiency histories and normalized saturation distributions (see for instance storage efficiency histories reported in Table 4; observed differences are less than one-tenth of one percent). We therefore omit the transverse permeability ratio from further mention in this text.

Table 4 Storage efficiency histories obtained for various transverse permeability ratio $\sigma _z$

Appendix 2: Limiting Case of No Crossflow

We calculate herein the limit of no crossflow the ratio of the shock-front velocities in each layer, then identify flow regions in which we can define a fast and a slow shock-front and finally predict the shape of the pressure profiles in each layer. We first confirm using numerical simulation (as described in Sect. 4.3) the occurrence of this zero crossflow regime for $R_L =0.01$ (e.g., Fig. 14). Saturation varies longitudinally in each layer, and the 1D saturation profiles along each layer are reported in Fig. 1c.

Fig. 14

Injected phase saturation at time $t = 0.46$ MPVI for a typical non-piston-like displacement through non-communicating layers (with $R_L =0.01$, $\sigma _x =5$, $R_s =1$ and $M_f^1 =M_f^2 =0.6$ ). The corresponding saturation as a function of distance along each layer is shown in Fig. 1c. The solution was obtained using numerical flow simulation, as described in Sect. 4.3

We begin with the calculation of the shock-front velocities. Integrating the total volumetric fluid flux per unit area along the model, using the incompressibility equation

$$\begin{aligned} \frac{\partial \hat{q}_{T,x} }{\partial \hat{x}}=0 \end{aligned}$$

(32)

obtained with $R_L =0$, yields the dimensionless pressure drop

$$\begin{aligned} {\Delta }\hat{P}\left( {\hat{t}} \right) =\frac{\hat{q}_{T,x} \left( {\hat{t}} \right) }{\hat{k}_{x,d} \overline{\lambda _T }\left( {\hat{t}} \right) }, \end{aligned}$$

(33)

where $\overline{\lambda _T }$ denotes the harmonic average of the dimensionless total mobility $\hat{\lambda }_T$ along the model. Equating the pressure drops along the two layers yields the ratio of inlet fluxes into layers 1 and 2,

$$\begin{aligned} \frac{\hat{q}_{T,x}^1 \left( {\hat{t} } \right) }{\hat{q}_{T,x}^2 \left( {\hat{t} } \right) }=\sigma _x \frac{\overline{\lambda _T }^{1}\left( {\hat{t} } \right) }{\overline{\lambda _T }^{2}\left( {\hat{t} } \right) }. \end{aligned}$$

(34)

The ratio of influxes is initially weighted by the end-point, displaced phase permeabilities $\left( {k_x k_\mathrm{rd}^e } \right)$ but varies with time in response to total mobility changes in each layer. This is a consequence of the uniform pressure boundary conditions. The hyperbolic saturation equation obtained in the no-crossflow limit,

$$\begin{aligned} \hat{C}_s \frac{\partial s}{\partial \hat{t} }+\hat{q}_{T,x} \left( {\hat{t} } \right) \frac{\hbox {d}f_i }{\hbox {d}s}\frac{\partial s}{\partial \hat{x} }=0, \end{aligned}$$

(35)

shows that the dimensionless shock-front velocity $\hat{U}_j$ in layer $j=1,2$, is related to the total flow rate as

$$\begin{aligned} \hat{U}_j =\frac{\left[ {\hbox {d}f_i /\hbox {d}s} \right] \left( {s_f } \right) }{\hat{C}_s }\hat{q}_{T,x}^j =\frac{1}{\hat{C}_s s_\mathrm{av} }\hat{q}_{T,x}^j , \end{aligned}$$

(36)

the latter equality making use of the relation $s_\mathrm{av} =1/f_i {\prime }\left( {s_f } \right)$. We obtain the shock-front velocity ratio as follows,

$$\begin{aligned} \frac{\hat{U}_1 \left( t \right) }{\hat{U}_2 \left( t \right) }=\frac{\sigma _x }{R_s }\frac{\overline{\lambda _T }^{1}\left( t \right) }{\overline{\lambda _T }^{2}\left( t \right) }. \end{aligned}$$

(37)

This ratio indicates that shock-front velocities change with time in response to mobility changes. When $\sigma _x /R_s >{\max }_{\mathrm{s}} \hat{\lambda }_T^2 \left( s \right) /{\min }_{\mathrm{s }} \hat{\lambda }_T^1 \left( s \right)$, the latter equation indicates that the shock-front always moves faster through layer 1. The converse is true for $\sigma _x /R_s <{\min }_{\mathrm{s }} \hat{\lambda }_T^2 \left( s \right) /{\max }_{\mathrm{s }} \hat{\lambda }_T^1 \left( s \right)$. Between these two cases, it is not necessarily possible to define without ambiguity a fast and a slow shock-front; for example, when the ratio $\overline{\lambda _T }^{1}/\overline{\lambda _T }^{2}$ alternatively takes values above and below one, and the ratio $\sigma _x /R_s$ is sufficiently close to one, such that an initially slower shock-front becomes faster than the initially faster shock-front.

We now predict the shape of pressure profiles along non-communicating layers. Along each layer, the longitudinal pressure gradient is calculated from the dimensionless Darcy’s law

$$\begin{aligned} \frac{\partial \hat{P}}{\partial \hat{x} }=-\frac{\hat{q}_{T,x} }{\hat{k}_{x,d} \hat{\lambda }_T }. \end{aligned}$$

(38)

Ahead of the shock-front, saturation remains uniform (at its initial value), so the total mobility $\hat{\lambda }_T$ does not vary and pressure decreases linearly. However, this is not necessarily the case behind the shock-front, where saturation varies if the displacement is not piston-like. To approximate pressure profiles and predict crossflow regimes for non-piston-like displacements, we assume the total mobility $\hat{\lambda }_T$ behind the shock-front equals total mobility at the shock-front $\hat{\lambda }_T{\left( {s_f } \right) }$, so pressure decreases linearly behind the shock-front. The pressure profile can therefore be described as linear by parts, with the ratio of pressure gradients ahead and behind the shock-front being equal to the shock-front mobility ratio. Uniform pressure boundary conditions at the inlet and outlet allow direct comparison of the pressure profiles obtained in the two layers as shown in Fig. 6. In the limit of piston-like displacements, the shock-front mobility ratio is equal to the end-point mobility ratio and the result becomes exact.

Appendix 3: Limiting Case of Equilibrium Crossflow

We calculate herein the limit of equilibrium crossflow the ratio of the shock-front velocities in the two layers and establish quantitative boundaries between flow regions in which one shock-front moves faster than the other. For large effective aspect ratios ($R_L \gg 1)$, a regular asymptotic expansion with respect to $R_L$ shows that

$$\begin{aligned} \frac{\partial \hat{P}}{\partial \hat{z}}=0 \end{aligned}$$

(39)

at leading order (transverse equilibrium; Yortsos 1995). Calculating the total fluid flux through a transverse cross section of the model,

$$\begin{aligned} \mathop \int \limits _0^1 \hat{q}_{T,x} \hbox {d}\hat{z}=-\frac{\partial \hat{P}}{\partial \hat{x} }\mathop \int \limits _0^1 \hat{k}_{x,d} \cdot \hat{\lambda }_T \hbox {d}\hat{z}, \end{aligned}$$

(40)

which is also equal to one by total fluid volume conservation, yields the dimensionless total flux in the longitudinal direction

$$\begin{aligned} \hat{q}_{T,x} =\frac{\hat{k}_{x,d} \cdot \hat{\lambda }_T }{\mathop \smallint \nolimits _0^1 \hat{k}_{x,d} \cdot \hat{\lambda }_T \hbox {d}\hat{z}}. \end{aligned}$$

(41)

In what follows, we assume there are no transverse variations in saturation within each layer (1D flow models within layers). We also consider, without loss of generality, that the shock-front moves faster through layer 1. The longitudinal total flux along this fast shock-front in layer 1 is given by

$$\begin{aligned} \hat{q}_{T,x}^1 =\frac{2\sigma _x M_f^1 }{M_f^1 \sigma _x +1}, \end{aligned}$$

(42)

and, using the incompressibility condition $\left( {\hat{q}_{T,x}^1 +\hat{q}_{T,x}^2 } \right) /2=1$, the longitudinal total flux along the slow shock-front in layer 2 by

$$\begin{aligned} \hat{q}_{T,x}^2 =\frac{2}{M_f^1 \sigma _x +1}. \end{aligned}$$

(43)

The shock-front velocities are obtained by dividing the total fluxes at the shock-front by $\hat{C}_s s_\mathrm{av}$. Therefore, we find the ratio of the shock-front velocities to be constant and given by

$$\begin{aligned} \frac{\hat{U}_1 }{\hat{U}_2 }=\frac{\sigma _x }{R_s }M_f^1 . \end{aligned}$$

(44)

Numerical solutions confirm the occurrence of the equilibrium crossflow for $R_L =100$ (see Fig. 15). These also indicate that the unidimensional flow model captures accurately the fast shock-front position but only provides an approximate position for the slow shock-front: The numerical solution shows injected phase flowing ahead of the slow front position.

Fig. 15

Injected phase saturation at time $t=0.5$ MPVI for typical displacements through perfectly communicating layers (with $R_L =100$, $\sigma _x =5$ and $R_s =1$ ) for a $M_f^1 =M_f^2 =0.6$ (where white dashed lines represent shock-front positions calculated analytically) and b $M_f^1 =0.4$ and $M_f^2 =1$ (where the unidimensional flow model predicts that shock-fronts move at equal velocities). The solutions were obtained using numerical flow simulation, as described in Sect. 6

In the flow region $\sigma _x /R_s >1$, the latter formula indicates that the fast shock-front is located in layer 1 if $\sigma _x /R_s >\hbox {max}(1,1/M_f^1 )$. Likewise, repeating the same calculation assuming that the fast shock-front is located in layer 2, we find in the flow region $\sigma _x /R_s <1$ that the fast shock-front is located in layer 2 if $\sigma _x /R_s <\hbox {min}(1,M_f^2 )$. In the flow region $\min \left( {1,M_f^2 } \right)<\sigma _x /R_s <\hbox {max}(1,1/M_f^1 )$, the two shock-fronts move at equal velocities.

Appendix 4: Unit Mobility Ratio Piston-Like Displacements

In this limiting case, there are no transverse contrasts in total fluid mobility, so no crossflow occurs between layers ($M_e =1$; Zapata and Lake 1981) and flow is described by the no-crossflow limit outlined in Appendix 2. Here we calculate storage efficiency at 1 MPVI as function of the dimensionless groups $\sigma _x$ and $R_s$. Without loss of generality, we consider the case $\sigma _x >R_s$ so layer 1 is the fast layer. Other cases can be treated in a similar fashion. At $t=1$ MPVI, the fast layer 1 is fully swept while only a fraction $\left( {{\Delta }\hat{x} } \right) _2 =\hat{U}_2 \le 1$ of the slow layer 2 is swept. The slow shock-front velocity is given by

$$\begin{aligned} \hat{U}_2 =\frac{2}{1+\sigma _x }\frac{1}{\left( {\hat{C}_s s_\mathrm{av} } \right) _2 }. \end{aligned}$$

(45)

The displacement being piston-like, we have $s_\mathrm{av} =1$, and the storage efficiency at 1 MPVI can be calculated as follows by

$$\begin{aligned} E_s =\frac{\int \int s\hat{C}_s d\hat{\hbox {x}} d\hat{\hbox {z}}}{\int \int \hat{C}_s d\hat{\hbox {x}}d\hat{\hbox {z}}}=\frac{1}{2}\left[ {\left( {\hat{C}_s } \right) _1 +\left( {{\Delta }\hat{x} \cdot \hat{C}_s } \right) _2 } \right] =\frac{R_s }{R_s +1}+\frac{1}{\sigma _x +1}. \end{aligned}$$

(46)