Upscaled Dynamic Relative Permeability for Unstable CO2 Flow in Stratified Porous Media

In geologic porous media, layering is ubiquitous on all length scales ranging from sub-mm compositional laminations to km-thick formations. The current study presents a general framework for estimating the two-phase viscous limit dynamic relative permeability (kr)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(k_{r} )$$\end{document} for non-communicating stratified sedimentary rocks. The approach is a two-stage upscaling taking into account the influence of nested heterogeneity at intra- and inter-layer scales. Sub-layer heterogeneities, promoting microflow instabilities in individual layers, are modelled by tuning input relative permeability parameters based on the viscosity contrast between fluids. Macroscopic flow instability at the domain scale is tackled by a new analytical solution for the saturation distribution in stratified media considering the complexities of frontal advance theory. For each layer, three flow stages are considered, each with its unique time-dependent behaviour. The overall solution is presented as a set of equations derived by overlapping of the solutions in different layers. The upscaled dynamic kr\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$k_{r}$$\end{document} is figured out by history matching. The practicality of the method is demonstrated in application to measurements collected from the Otway International Research Facility, Victoria, Australia. The results show the significant impact of layer separation on the upscaled endpoint saturations in addition to required modifications on the current kr\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$k_{r}$$\end{document} formulas. A comparison with conventionally modelled steady-state sweep highlights the inadequacy of steady-state upscaling for viscous limit flows. The applicability of the method to estimate kr\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$k_{r}$$\end{document} from the unsteady-state lab experiments conducted on heterogeneous cores is also discussed.

1 3 k r Relative permeability l Flood front distance from injection boundary (L) M Mobility ratio N Number of layers p, p in and p out Pressure, inlet pressure, and outlet pressure (ML − 1 T − 2) q inj Injection rate (L 3 T − 1) s and s r Saturation and residual saturation t Time (T) w Width (L) x and X Position measured from inlet boundary and total model length (L)

Introduction
Multiphase flow upscaling is considered as one of the most challenging steps in the model building and reservoir simulation workflow. The existence of more than one flowing phase requires the upscaling of saturation functions, i.e., capillary pressure and relative permeability (e.g., Fadili and Ababou 2004;Pickup and Serbie 1996). Relative permeability upscaling has been an active research area since the 1970s and a large number of approaches have been devised for this purpose (e.g., Pickup and Sorbie 1996;Darman et al. 2000;Virnovsky et al. 2004;Lohne et al. 2006). Most recent upscaling efforts focus on sub-core scale heterogeneity and X-ray tomography (e.g., Perrin and Benson 2010). Comparatively, little work is conducted on continuum-scale heterogeneity (e.g., Sedaghat et al. 2017). One of the traditional approaches that is still attractive is classifying continuumscale heterogeneity in terms of combinations of rock types at the sub-core scale (Boon et al. 2022).  studied laminated composites while Lohne et al. (2006) modelled fine-scale heterogeneity as a checkerboard structure. The idea of utilizing upscaled dynamic relative permeability to describe unsteadystate flow behaviour at the reservoir scale is not new (Hawthorne 1975;Hearn 1971) and has been widely investigated in the literature (e.g. Kossack et al. 1990). Yet, in this century, steady-state upscaling studies appear to dominate over dynamic ones. Arguably, steady-state methods are more popular since they are simpler to calculate saturation distributions in heterogeneous porous media for both, the capillary limit (CL, cf. Corey and Rathjens 1956;Ekrann et al. 1996) and the viscous limit (VL, cf. Ekrann and Aasen 2000). Unsteady-state upscaling is riddled with the complexity of space-time-dependent behaviour (Potashev 2014). These challenges persist even at the VL limit (see Debbabi et al. 2017). Furthermore, preparation of the upscaled data for coarse-scale simulation requires extensive sensitivity analyses that may have to be repeated at consecutive length scales (Pickup and Stephen 2000). Moreover, the derived dynamic relative permeability curves may exhibit non-standard behaviour like non-monotonicity. The reader is referred to Barker and Dupouy (1996) for a quantitative comparison between widely used techniques to estimate upscaled dynamic relative permeability and Barker and Thibeau (1997) for a qualitative one.
The vast majority of recently published upscaling studies rely on local steady-state upscaling (c.f. Durlofsky 2005), at least for the first stage, estimating saturation functions at their limits (e.g. . Others use a specific pressure gradient that is deemed relevant for the production of a particular reservoir (Lohne et al. 2006). Yet, recent experiments conducted by Zhang et al. (2017) cast doubt on the applicability of upscaled steady relative permeability curves to unsteady-state flow situations. Zhang concluded that co-current co-injection-induced flow, considered by steady-state upscaling, is a different displacement mechanism than the sweep of one phase by another in a reservoir. This assessment also matches the observations on micro-displacement efficiency by Shetty et al. (2014) who compared water alternating gas injection (WAG) with simultaneous water and gas (SWAG) injection.
Arguably, laminations and layering constitute the most commonly formed heterogeneity of sedimentary rocks across multiple length scales. Using a medical CT scanner, among others, Zhang et al. (2014Zhang et al. ( , 2017 demonstrate how strongly the distribution of immiscible fluids is affected by this structure causing flow channeling along the layers. This was recognized early, motivating the mathematical descriptions that can capture this behaviour. Dykstra and Parsons (1950) presented the first model of non-communicating layered media. Hearn (1971) generated the first upscaled dynamic relative permeability curves for stratified porous media under the assumption of complete crossflow between the layers (vertical equilibrium). All following analytical approaches (e.g., Pande et al. 1987) ignored frontal advance theory (Buckley and Leverett 1942), assuming that saturation changes only at the saturation front. El-Khatib (2001) was the first to consider the effect of the rarefaction wave from frontal advance theory in his study of non-communicating stratified media. However, his solution is limited to the case where all layers have the same relative permeability and initial saturation.
In this article, we analyse the unsteady-state relative permeability of a stack of noncommunicating porous layers, considering the variation of flow properties among them. CO 2 geo-sequestration in a layered non-fractured aquifer with an injection well (line drive) is considered as an application for the results of our analysis. The objectives of this study are four fold. First, we want to develop an analytical solution to the noncommunicating parallel layer water flooding problem, revisiting the original work of Dykstra and Parsons (Dykstra and Parsons 1950). Second, this analytical solution will be employed to construct upscaled dynamic unsteady-state relative permeability curves. Third, results obtained with this new unsteady-state method will be compared with a widely used steady-state solution. Fourth, the accuracy of the experimental JBN relative permeability determination method (Johnson et al. 1959) will be evaluated on the same field data-based system.

Mathematical Model
The mathematical model developed in this study considers multiphase flow in N non-communicating layers interconnected by a well at the inlet (Fig. 1). Although applicable to a wide range of scenarios, it is illustrated here for CO 2 injection into a storage site. Supercritical CO 2 is injected at a rate q inj (t) while the outlet boundary is held at constant pressure, p out .
The stratified model is composed of N layers, each characterized by its own absolute permeability (k), porosity ( ), relative permeability (k r ), and residual saturation (s r ). It is assumed that there is no cross flow between layers. This assumption is valid in the presence of shale interlayers (Xu et al. 2020) or if the vertical permeability is negligible compared to the horizontal one. The reader is referred to the analysis performed by Debbabi et al. (2017) for a more discussion of this assumption. Also, it is assumed that both the gravitational and the capillary forces are ignored. The pressure equation in its dimensionless form of an arbitrary layer i ∈ {1, 2, ..., N} , derived from the conservation law for volume and appropriate boundary condition reads as follows The saturation equation (Chen 2007 Saturation distributions (dotted lines) before boundary breakthrough in a linear (not radial) multilayer flow model. The Buckley-Leverett profiles mark the advancing injection fronts in the layers followed by the characteristic rarefaction wave. Each layer is characterized by its characteristic frontal advance (l) and petrophysical parameters, see notation table and saturation functions. Like in the Dykstra-Parsons' model, the layers are arranged in descending order from fastest to slowest. Injection rate inflow on the left (q inj ) , partitions into the layers according to the average total mobility of each one. The outlet (right) boundary has a fixed pressure (p out ).
where t denotes time, p fluid pressure, with subscript p in for the inlet boundary. The derived variables are the total mobility ( T ) , the fractional flow (f ) , and the storage capacity (c s = i Δs i ) , while Δs = 1 − s ini CO 2 − s r w quantifies the saturation range bounded by initial CO 2 saturation (s ini CO 2 ) and 1 − s r w . The product k e = kk e r w is the endpoint effective water permeability, and k e r w the water endpoint relative permeability. The superscript ^ indicates variables in dimensionless form (see Appendix A for their definition). Importantly, saturation-dependent parameters are defined as functions of s CO 2 instead of brine saturation as CO 2 is the injected phase which makes it more convenient for the application of Welge's tangent method for frontal advance (Welge 1951).
Solving this elliptic PDE subjected to the Dirichlet boundary condition in the system (1), the pressure drop (difference between inlet and outlet pressure) is the same across any layer and can be quantified as where h is the thickness, and ̂a vg T is the dimensionless average total mobility that is calculated as the harmonic mean of the layers

Flow stages
The key element of our mathematical analysis is how to correlate the time-dependent layer flow rate to a suitable time-dependent variable such that the corresponding average total mobility of the layer is a function of it. If done correctly, the system's effective total mobility can be evaluated at any instant during the flooding process. For this purpose, we distinguish three sequential flow stages in each layer. In each stage, the flow rate follows a characteristic relation with a specific time-dependent variable. The flow stages are distinguished as follows: 1. shock propagation: the passage of the saturation front through the layer (l i ≤ 1); 2. rarefaction wave: the outlet saturation gradually increases up to ŝ i = 1; 3. time span from when the layer reaches full saturation to the end of the injection process attained when the whole system reaches its maximum allowable saturation.
The stages are reached independently on the variation of the boundary conditions. The overlap between the stages, in different layers, can be calculated to obtain a series of solutions for the overall saturation change of the system. The characteristic variable ( ) of the stratum i in each stage is ̂l i ,ŝ * i and Q i , respectively, where ̂l is the flood front position, ŝ * is the saturation at the production boundary, and Q is the cumulative injection volume. Once the performance of the whole system has been established, its upscaled dynamic relative permeability can be calculated.

Stage 1: before breakthrough
Before CO 2 breaks through the outlet boundary, all layers are in the first stage. The flow rate in any layer is derived from frontal advance theory: The relation between the shock front in the i th layer and the fastest (first) layer is dŝ 2 , and is Kronecker delta (see Appendix B for the full derivation). In this study, reflects the mobility ratio (M i =̂e CO 2 i ) Although the solution (6) is given as a function of flood front advance in the fastest layer, it can be easily modified to correlate any two different layers in terms of their flood front locations. We prefer to correlate with the fastest layer as it simplifies the overall algorithm that is described later.

Stage 2: following breakthrough
When the flood front in a layer reaches the far side, the outlet saturation starts to increase above its shock saturation. Given the aforementioned ordering of the layers, this process evolves in a sequential manner starting with the fastest layer and continuing until all layers reach breakthrough. After breakthrough, the interest shifts to the outlet end saturation or the cumulative injection volume of this layer.
The following analysis considers possible combinations of layer stages after the breakthrough but before the third stage. Eq (5) is modified for those layers that reached breakthrough such that their flow rate reflects the second stage (1 ≤ j ≤ n < i ≤ N) described as follows is the dimensionless saturation at the outlet boundary. The performance relation between the fastest layer and these layers indexed with j is is the dimensionless outlet saturation of the fastest layer at the time of breakthrough in layer j and is a placeholder variable for saturation. Eq (9) reveals that the outlet saturation of a second-stage layer is an implicit function of the corresponding outlet saturation of the other second-stage layers. Turning our attention to the relation between the i th layer that still is in stage one and the fastest layer, it can be shown that (Appendix B) where ̂l i | | |l 1 =1 is the flood front position in layer i at the time of the breakthrough in the fastest layer. ̂l i | | |l 1 =1 is estimated using Eq (6) by assigning ̂l 1 = 1 . Equation (10) estimates the locus of the flood front positions in first stage layers as a function of the dimensionless saturation at the outlet of the first layer. Again, a special case arises when M i = 1 and Eq (10) reduces to Continued injection leads some layers into the third stage. While m layers are in the third stage, (n − m) layers are in the second stage, and the remaining layers are still evolving through the first stage (1 ≤ a ≤ m < j ≤ n < i ≤ N) . The flow rate of any layer that is in the third stage is where Q is the dimensionless injected volume. In this scenario, the relation between layer 1 and an i layer that is still in the first stage is given by is the shock front advance in layer i and Q 1 | | |ŝ 1 =1 is the cumulative injected volume in the fastest layer at the beginning of stage three in the fastest layer.
is estimated from the Buckley-Leverett solution. Equation (13) describes the relation between the fastest layer, which is now in the third stage, cumulative injected fluid volume and the shock-front advance in a layer still at the first stage. On the other hand, the relation between the fastest layer and a layer j at the second stage can be described as follows where ŝ * |̂s 1 =1 is the outlet saturation at the start of the third stage in the first layer. Equation (14) shows how the outlet saturation of the layers at the second stage varies with the cumulative injected volume of a layer at the third stage.
At any instant either before or after breakthrough, the dimensionless time (t) can be estimated from the conservation of volume as follow The right-hand side of Eq (15) consists of three sums: The first one on the left pertains to the layers at the third stage, the middle one concerns the layers at the second stage, and the right one accounts for the layers in the first stage.
In summary, the time-dependent saturation in the different layers is estimated following these steps: 1. arrange layers in descending order from the fastest to the slowest one using Eq (6) iteratively; 2. track the position of the shock front in the fastest layer ( ̂l 1 ) until breakthrough; 3. estimate the loci of the shock fronts in the remaining layers using Eq (6) or Eq (7) based on the mobility ratio within the layer of interest; 4. specify different values of ŝ * 1 between ŝ f 1 and one; 5. estimate the corresponding positions of the shock front at first stage layers using Eq (10) or Eq (11) from the mobility ratio, and estimate the outlet saturation of second stage layers using Eq (9); 6. use the same approach to correlate between Q 1 and the characteristic parameters of the other layers based on the evolution stage that they are at using Eqs (13) and (14); 7. Evaluate the corresponding t at each step using Eq (15).

Upscaled dynamic relative permeability
At a first glance, the abstraction of the upscaled relative permeability from the described relationships appears like a straightforward task. However, this is not so due to the complex mobility variations in the different layers, manifest in the changing layer pressures, and pressure gradients. It follows that purely analytical approaches (e.g., Potashev 2014) are suitable only for estimating the upscaled relative permeability of a slice of layers and not for their ensemble. A workaround is finding the upscaled relative permeability by history matching of results that are readily generated analytically. As the upscaled relative permeability is generated considering the passage of the shock fronts through the model, it is named "dynamic relative permeability". Before proceeding to history matching, additional analysis is warranted, if necessary, to reduce the number of parameters that need to be matched. Here we perform, two preprocessing steps to treat flow instabilities at different scales; (1) the sub-layer scale, and (2) the multi-layer system scale.

Micro (sub-layer) instability treatment
Our analytical model relies on unique distinct petrophysical parameters for each layer. However, in nature, it is more common to have property variations at the sub-layer scale. Bear (1972) suggested that below the characteristic length of representative elementary volume (REV), petrophysical parameters fluctuate between two extremes; an upper bound and a lower one. These limits converge as the REV scale is approached. Small-scale heterogeneity also promotes unstable displacement including channeling. This phenomenon has a strong impact on the upscaled dynamic relative permeability. Using simulation, Doorwar and Ambastha (2020), recently showed that the viscosity ratio between displaced and injected fluid ( w CO 2 ) is strongly correlated with the modification that is applied on the finescale (input) relative permeability parameters to incorporate the instability. In this study, we follow their recommendations to modify the input relative permeability of each layer before using it in our analysis. Doorwar and Ambastha (2020) proposed that to account for the effects of channeling, the displaced phase relative permeability (k r w ) needs to be modified while the displacing phase k r CO 2 can be retained. They also observed that the residual saturation of the displaced phase (s r w ) needs to be increased following a power-law trend with ( w CO 2 ) to yield the same mobility ratio as for the actual (theoretical) s r w . For unstable displacement, at any w CO 2 > 1 this empirical correlation yields the same mobility ratio as stable displacement at w CO 2 = 1 . Figure 2 presents an application of the Doorwar and Ambastha (2020) approach to estimate adjusted s r w . In this example, CO 2 w at s r w is set to 10 6 which corresponds to a f = 99.9999 % at unit viscosity ratio. This value of fractional flow is selected such that the adjusted s r w at unit mobility ratio is within acceptable proximity to the theoretical s r w . Figure 2b confirms that the achievable s r w follows the power-law trend as expected. The empirical correlation follows where a is a fitting parameter. This fitting process forces the achievable s r w to return the theoretical s r w at that unit viscosity ratio. The other parameter affected by unstable displacement is the curve shape parameter (e.g., Brooks-Corey's parameter in the Brooks-Corey model). Its modification is accomplished by matching with the fine-scale simulation results obtained for the heterogeneous rock. In this study, however, the curve shape parameter is not modified as the goal is to present a general framework for generating dynamic relative permeability curves for non-communicating layered media. Separate studies are needed to justify and identify the functional form of this modification based on pore-scale displacement patterns at different flow rates.

Macro (multi-layer) instability treatment
The contrast between layer properties (Fig. 1) and the degree of heterogeneity in individual layers directly affects the propagation velocity of saturation fronts. For large reservoir volumes, thousands of pore volumes need to be injected to sweep the whole system. This is not attainable in practice. It follows that the achievable s r w of the whole system will be higher than the analytical s r w that is reached when all the mobile fluid has been produced. At the system scale, the achievable s r w can be evaluated using Jones-Roszelle's graphical method (Jones and Roszelle 1978). Jones and Roszelle (1978) showed that linear extrapolation of the displaced fluid (brine) saturation (< s w >) to infinite injected pore volume ( 1 t → 0) gives a realistic estimate of the practical residual saturation. For unstable displacement, this approximation deviates from the analytical s r w and produces an achievable value that is more suitable for the upscaled dynamic k r (Luo et al. 2016).
Eventually, s r w is subjected to two modifications made to address micro and macro instabilities. For both, the term "achievable residual saturation" is adopted. The reader can easily distinguish which type of instability this qualifier refers to from the context in which it is being used.

History matching
The purpose of the upscaled dynamic k r is to mimic the overall behaviour of the system. Usually, the history matching is performed on production data (pressure drop and average CO 2 saturation), and, in case of availability, the saturation profile (Berg et al. 2021a(Berg et al. , 2021b. In our analysis, we avoid using the saturation profile. To evaluate the accuracy of the matches made herein, the least square error is defined as the normalized squared error sum of pressure and average saturation data as follows (c.f., Berg et al. 2021a) where Δp is the pressure drop and < s CO 2 > is the average CO 2 saturation. This normalization balances the contribution of the pressure drop squared error to 2 with the contribution of the saturation average squared error to 2 .

5. Results
The analysis methodology is applied to field data collected from the CRC-3 well at the Otway International Research Facility, Victoria, Australia (e.g., Mishra et al. 2020). Detailed air mini-permeameter measurements collected on split core nearly every 10 cm, along the injection interval of the Paaratte formations parasequence 1 (several coarsening-downward sedimentation cycles buried in the Otway basin underneath shale formations, Boon et al. 2022) show orders of magnitudes variations in the permeability (Fig. 3). The geo-statistical permeability distribution across the site indicates that the mode of the vertical permeability of mudstone, considered a non-reservoir volume, is around 1 mD (not shown). This value is considered a cut-off assuming that the vertical permeability of the mudstone is the smallest among all possible directions. Applying Variation of permeability estimated from mini-perm measurements (solid blue) with depth along injection interval in CRC3. The permeability cut-off value (1 mD) is shown by the red dotted line while red solid lines show the flowing intervals and their average permeability this cut-off value yields 11 flow units, with different thicknesses and average permeabilities, separated by streaks of mudstone (Fig. 3). Table 1 summarizes the relevant essential parameters needed to estimate their commingled frontal advance. The estimated Dykstra-Parsons coefficient, considering only the flowing layers, is 0.5 which reflects a moderate degree of heterogeneity. The fine-scale relative permeability is estimated using the Brooks-Corey model (Brooks and Corey 1965) where is the Brooks-Corey's parameter and k r CO 2 (s w = s w r ) is the CO 2 endpoint k r ( Table 1). All the flowing layers are assumed to have = 0.28, s r w = 0.159, s ini CO 2 = 0 and s co 2 r = 0.05 . Two scenarios are considered: (1) stable drainage with w CO 2 = 1 , and (2) an unstable drainage with w CO 2 = 30 . Table 1 also shows the estimated fitting parameter a appearing in Eq (16) which is used to calculate the achievable residual brine saturation accounting for potential sub-layer instabilities.
Considering the CO 2 -brine viscosity ratio and after shifting s r w to the new achievable values in the presence of potential micro-instabilities, the layers undertake the order in Table 1. Figure 4 shows the successive overlap of the flow stages in different layers with the first layer flow stages. Figure 4a shows the performance when all layers are in the first stage. After that, the fastest layer enters the second stage following its breakthrough while all other layers are still in Stage 1 (Fig. 4b). In the end, all layers reach breakthrough, and the focus shifts to the outlet saturations (Fig. 4c). It is clear from Fig. 4a that an increase in the viscosity ratio would magnify shock front propagation velocity differences across the layers. Moreover, Fig. 4b indicates that dl i dŝ * 1 directly after the breakthrough of the first layer is much steeper when the viscosity ratio is 1 compared to 30. These observations again confirm that an increase in the viscosity ratio increases the propagation velocity differences between layers. Figure 5 presents the variation in the average brine saturation with the inverse of dimensionless time. Extrapolation of the early-time linear best fits determines macroscopic achievable residual brine saturations as 0.31 and 0.57 for the viscosity ratios 1 and 30, respectively. The plot shows the important impact of the fluid viscosity contrast on the sweep even if micro-instabilities are ignored.
To introduce more flexibility to the history matching procedure, the upscaled dynamic k r is generated using a more generic form of the ]. This model is modified by allowing the empirical exponent parameter to vary between for fitting k r w and k r CO 2 , i.e., w may be different from CO 2 . Optimization was performed using the 'lsqnonlin' function embedded in MATLAB (MATLAB and Statistics Toolbox Release, 2021b, The MathWorks, Inc., Natick, Massachusetts, United States) and applying the 'trust-region-reflective' algorithm. The upper limit of both w and CO 2 was set to 10. The  Table 2. The comparison between the upscaled dynamic k r with the input k r shown in Fig. 6 reveals that the upscaled dynamic k r is shifted to higher s w . This shift increases as w CO 2 increases. Compared to the analytical solution, the upscaled k r predicts both the pressure drop and the breakthrough time with reasonable accuracy (Fig. 7). However, for both viscosity ratios, the fractional flow estimated at the outlet is higher than that expected from the analytical solution. Consequently, the upscaled dynamic k r diminishes the volume of CO 2 stored after breakthrough (Fig. 7c). Moreover, its smooth shape is unsuitable for predicting the discontinuous stair-step nature resulting from sequential breakthroughs (Fig. 7b). These results suggest that a more flexible functional form is needed to represent the ensemble relative permeability. Discussion about this point follows in the analysis of the laboratory experiment (objective 4). The results corroborate the assessment already made by Durlofsky (1997) on the basis of numeric experiments. He agreed that "the traditional functional form of relative permeability is too limited to capture behaviour" and "an extended description is needed taking into account the observed velocity -saturation covariance".  . 6 Upscaled dynamic k r for the two viscosity contrasts considered as compared with input k r of flow zones (Table 1) shown in grey in the background Fig. 7 a Dimensionless pressure differential between inlet and outlet, b outlet CO 2 fractional flow, and c average CO 2 saturation, generated using pseudo k r (dashed lines) and unsteady-state analytical solution (solid lines) at the two viscosity ratios considered

Comparison with steady-state upscaling
Steady-state upscaling is widely applied in upscaling studies of heterogeneous porous media (e.g., Dale et al. 1997;Virnovsky et al. 2004), with VL parameters generated assuming that fractional flow is constant along streamlines. It follows that this approach is highly sensitive to inlet boundary conditions. Usually, a constant-uniform fractional flow is assigned to the inlet (e.g., Ekrann and Aasen 2000). As the application target for our analysis is the unstable displacement of brine by supercritical CO 2 , it is instructive to evaluate the adequacy of steady-state upscaling to this problem in comparison with the unsteady state method proposed herein. For flow in the horizontal direction, this comparison is legitimate because the non-communicating layers enforce parallel flow. Figure 8 shows the steady-state upscaled k r of the layered sequence from Fig. 3. Unlike the dynamic k r , the steady-state k r is almost insensitive to the fluid viscosity contrast and controlled primarily by input k r . Furthermore, to avoid non-uniqueness, the steady-state k r must be generated between fixed endpoint saturations and then extrapolated from the critical gas saturation to s w = 1 . Its usage always overestimates breakthrough time (Fig. 9), and micro-displacement efficiency as compared to the values estimated by the proposed unsteady-state analysis. Figure 10 shows that the difference in < s CO 2 > between the two approaches increases with the reduction of w CO 2 . Also, < s CO 2 > estimated by unsteady-state analysis correlates with w CO 2 following the correlation function Fig. 8 Upscaled steady-state relative permeability at different fluid viscosity contrasts as compared with the input k r for the flow zones (Table 1) shown in grey Fig. 9 a Dimensionless pressure differential between inlet and outlet, b outlet CO 2 fractional flow, and c average CO 2 saturation, generated using the upscaled steady-state k r (dashed lines) and the unsteady-state analytical solutions (solid lines) for the two viscosity ratios considered where a and b are fitting parameters, < s CO 2 > max and < s CO 2 > min are the maximum and the minimum achievable average saturation at breakthrough, respectively. These saturation extremes can be determined analytically.

Laboratory (special core) analysis
Dynamic core flooding, often referred to as JBN method (Johnson et al. 1959), is often used to determine relative permeability experimentally. This widely used approach automatically homogenizes heterogeneous cores. CT imaging of laminated cores during such experiments indicates preferential flow in certain layers (e.g., Zhang et al. 2014). As the current analysis does not consider cross-flow between layers, it is possible only to compare the JBN approach as compared with the Potashev approach (Potashev 2014). The latter is a convention for analysing unsteady-state flow experiments conducted on such types of cores. Both approaches assume the dominance of the viscous force over both capillary and gravity forces. This is justifiable for unsteady-state experiments conducted at high flow rates suppressing the impact of capillary force while the influence of gravity force is nonsensible over the short and thin cores.

Potashev approach
The approach proposed by Potashev determines dynamic relative permeability of slices through a stack of non-communicating layers. Usually, the slice location is at the outlet of the sample (Johnson et al. 1959). Such results are useful also for the determination of fractional flows on the field scale (Hearn 1971). The Potashev algorithm for the estimation of dynamic relative permeability comprises the following steps: 1. Estimate average saturation as < s >= 3. Determine average total mobility as < T >= Fig. 10 Variation of the CO 2 average saturation at the time of the breakthrough with the fluids viscosity ratio as determined by the unsteady-state analysis (blue dots) versus the upscaled steadystate values (black squares). The solid red curve represents the best-fit curve based on Eq (19) with fitting parameters a and b are 0.23 and -0.19, respectively 4. Calculate average phases mobilities as < CO 2 > = < f >< T > and From the estimated average mobilities, the dynamic relative permeability of the target slice is deduced by knowing fluid viscosity.

Comparison between JBN and Potashev
Neglecting the micro-instabilities, and assuming the outlet saturation is the same for both approaches and equivalent to the one estimated from the unsteady-state analysis presented (Fig. 4). By contrast with the steady-state k r (Fig. 8a), both JBN and Potashev techniques k r curves are not necessarily bounded by input k r (c.f., Hamon and Roy 2000). Figure 12 reveals that the differences between these unsteady-state approaches become more obvious as w CO 2 approaches unity. Figure 11 also shows that the unsteady k r of the stratified medium is characterized by multiple saturation jumps reflecting jumps in the outlet fractional flow corresponding to breakthrough events (Fig. 7b). As already highlighted by El-Khatib (2001), the unsteady k r can also show non-monotonic variation with saturation.

Discussion
In the current analysis, the mobility ratio is defined in terms of endpoint mobilities. This is not the only definition in the literature. Some authors prefer to use the mobility ratio at the flood front, especially for unstable displacement that arises from perturbations of the pressure field (cf. Berg and Ott 2012). Here, we stick to the endpoint ratio definition because it Fig. 11 Outlet k r estimated using Potashev (solid lines) and JBN (dashed lines) for the two viscosity ratios considered; input layer k r is shown in the grey Fig. 12 Impact of mobility ratio on the slope of the linear relation between reciprocal of equivalent total mobility and flood front location before breakthrough is more telling with respect to the change in average total mobility, determining the sign of the right-hand side of Eq (6). Figure 12 shows that the mobility ratio controls the slope of the linear relation between the dimensionless flood front location and 1 avg T i Eq (22). If the slope is positive (M < 1) , the left-hand side in Eq (6) is also positive and vice versa. Our definition of i ascertains that both sides have the same sign and the calculated ̂l i is always positive.
Equation (6) can be easily reduced back to Eq (A15) derived by El-Khatib (2001) by choosing i = 1 and assuming that all layers have the same relative permeability, porosity, and initial fluid distributions equivalent to the residual saturation. Figure 10 also retrieves the same average saturation at the time of breakthrough as predicted by Dykstra and Parsons' method (Dykstra and Parsons 1950) for w CO 2 → 0 , when there are only two flow zones separated by a shock front without any rarefaction wave. This permits a quick verification of the proposed solution.
In this study, a simple approach to include approximations of flow instabilities into VL is proposed, treating micro instabilities with the technique presented by Doorwar and Ambastha (2020). This approach is based on results obtained from flooding simulations conducted on a model with a bimodal permeability distribution. However, the displacement can also be unstable if there is no or little heterogeneity as in the case of viscous fingering (e.g., Berg and Ott 2012). Their study also confirmed that the shock front mobility ratio is a reasonable criterion to recognize when such type of instabilities should occur. Hence, it is recommended to study micro instability on each layer separately, whether homogeneous or heterogeneous, and modify input k r as necessary. Moreover, the current approach still fails to address the effect of coarsening instabilities which grow over time and correlates with the domain extent (Berg and Ott 2012).
Application of our unsteady-state method to field data from Otway shows that the permeability contrast between layers plays a significant role. The presented example demonstrates that macro instabilities occur even if the flow at the layer scale is stable. Figure 6 shows how the residual brine saturation can be modified to capture this behaviour even when there is no viscosity contrast between the fluids. In this case, the maximum permeability contrast is about 130. This is less than that used by Doorwar and Ambastha (2020). This consideration highlights the important effect of baffles separating the layers and preventing cross-flow as shown in Figure 13.
We propose to estimate the dynamic k r by history matching assuming that the functional form of k r (s w ) is known. Extended fitted Brooks-Corey functions do an excellent job in capturing pressure differentials except at breakthrough. Yet, the match quality is lower for the average saturation. That is mainly due to the piecewise nature of the analytical solution which cannot be predicted using a smooth relative permeability function. Even when more are ordered from fastest to slowest from top to bottom of the plot. For the sake of visualization, all the layers were assigned the same thickness fitting parameters are included as in the LET model (Lomeland et al. 2005), one always obtains a continuous fractional flow. By contrast, as seen in the field, fractional flow often changes in a stairstep fashion at gas or water breakthrough to a producer well. This suggests that rigorously upscaled k r (s w ) curves do not necessarily possess continuous differentiability. Also, the fitted parameters are pseudos only and do not reflect the same physics attained from homogenous media. Therefore, the upscaled k r (s w ) parameters can be out of the expected range, i.e. when the viscosity contrast equals 30, both w and CO 2 are higher than 7.3 which is the value recorded for mono-disperse glass beads (Brooks and Corey 1965).
This study only addressed the VL flow neglecting the impact of other forces; capillarity and buoyancy, assuming that viscous forces dominate them everywhere in the system. However, accounting for the influence of these forces is crucial to predicting the multiphase flow behaviour in non-idealized porous media. Thus, where the flow velocity in each layer is a product of viscous and capillary forces (Boon and Benson 2021), therefore, the achievable sweep will depend on injection velocity. This limitation will be addressed in future studies that extend the presented framework to more general situations. In general, the validity of the dominance of viscous force assumption is more relevant to relatively thin and high permeability layers where the influence of the gravitational force is non-sensible, and the flow velocity is large such that it suppresses the capillarity.

Conclusions
The current study combines techniques from the literature with a new analytic solution for 2-phase unsteady immiscible flow in non-communicating layered porous media. This formulation serves as a semi-analytical framework that facilitates the quick estimation of the dynamic ensemble k r of the system. This solution is constructed by observing that.
• multiphase flow through non-communicating layered media occurs in three stages, each characterized by a unique time-dependent group of variables; • flow instabilities at the sublayer or the domain scale, both shift the residual saturation of the displaced phase to higher values reflecting the maximum anticipated microsweep efficiency; • the complexities associated with multiphase flow in non-communicating layers reduce the accuracy of classical models for ensemble k r especially with regard to predicting system behaviour. Hence, new k r (s w ) formulations are needed to accurately capture ensemble behaviour; • the results confirm that the widely used Dykstra-Parsons approach applied to predict the performance of non-communicating stratified reservoirs is only adequate at a very low viscosity ratio between displacing and displaced fluids promoting separation by a sharp shock front.
For the Stage 2 layer, again utilizing the frontal advance theory, the reciprocal of the average total mobility is formulated as a function of the outlet saturation as follows Substituting by Eqs (24) and (26) into Eq (25) then applying the method of separation of variables and integration, the solution (10) can be easily figured out. Following the same technique, it is handy to deduce the solution (9).
For a layer that has passed the first and the second stage, the saturation along the whole layer is no longer changing. Therefore, its average total mobility is constant and is equivalent to the mobility ratio. Substituting by Eqs (12) and (5) or (8) into Eq (20), and again applying the separation of variables, solutions (13) and (14)