Upscaling multiphase viscous-to-capillary transitions in heterogeneous porous media

Abstract Upscaling the effect of heterogeneities in porous media is crucial for macroscopic flow predictions, with numerous applications in energy and environmental settings. In this study, we derive simple semi-analytical expressions for the upscaling of multiphase flow in a porous medium with a range of vertical heterogeneities. We use this upscaling to give insight into how the flow transitions between a viscous flow regime, in which macroscopic pressure gradients dominate over heterogeneity-driven capillary forces, and a capillary flow regime, in which these capillary forces dominate and set the saturation distribution of the flow. In particular, by studying the dynamics of flow in an aquifer, we demonstrate that different regions lie within the viscous and capillary flow regimes whilst other regions lie in between these regimes. By modifying the classic Buckley–Leverett problem for fluid displacement we demonstrate where and when the flow transitions between these regimes and how this affects flooding speeds. Then, we discuss the implications of these results in the case of carbon dioxide sequestration, making comparisons with field data.


Introduction
The flow of immiscible fluids in heterogeneous porous media has widespread applications in energy and the environment. Nearly all subsurface rocks have a significantly heterogeneous structure, often in the form of sedimentary layers, and it is well known that such heterogeneities play an important role in the resultant flow properties (Corey & Rathjens 1956;Jackson et al. 2018;Nijjer, Hewitt & For example, in two-phase flow, the non-wetting phase tends to be preferentially drawn to regions of larger pore space. The effect of the heterogeneities also depends on how they are distributed. The most common type of heterogeneity is sedimentary layering parallel to the aquifer and flow direction, but the complexities of the processes responsible for the deposition of sediments and their subsequent geological history frequently impose much more complex permeability structures. Hence, we focus on the simpler horizontally layered case. Additionally, in pressure driven flows, heterogeneities result in the unstable displacement of phases (so long as capillary forces are large enough to overcome the driving pressure), and fingering (Dawe, Wheat & Bidner 1992;Dawe, Caruana & Grattoni 2011). Hence, an analogy can be drawn between the capillary-driven mixing of immiscible fluids, and the classic diffusion/dispersion-driven mixing of miscible fluids (Tchelepi et al. 1993;Nijjer et al. 2019). However, for this study we focus on the case of immiscible fluid flow in a layered porous medium.
The role of heterogeneities is often characterised by the non-dimensional capillary number, which is given as the ratio between typical horizontal pressure gradients p/L (over length scale L), and typical vertical gradients in the pore entry pressure p e /H (over length scale H), giving (1.1) At small N c , the background flow is sufficiently weak that the flow of fluid phases is largely dominated by the heterogeneity-driven capillary forces, whereas at large N c , the background flow dominates, such that capillary forces due to heterogeneities can be largely ignored. Hence, the limit N c → 0 is known as the capillary limit and N c → ∞ is known as the viscous limit. To model the flow in any case which is far away from the viscous limit, one needs detailed knowledge of the structure of the heterogeneities to describe the flow, which presents a significant challenge.
Recently, there has been renewed emphasis on attempting to upscale the effect of heterogeneities in porous media Boon, Bijeljic & Krevor 2017;Jackson et al. 2018;Jackson & Krevor 2020). One of the key difficulties lies in the sheer number of measurements, either experimental or numerical, needed to characterise the effect of rock layers across a broad range of flow conditions. For example, in the case of immiscible flow of wetting and non-wetting phases, the effect of the heterogeneities not only depends on the capillary number, as described above, but also on the fractional flow of each phase (Woods 2015). Furthermore, since each type of rock heterogeneity is different, it is difficult to transpose results without performing experiments and simulations for each specific case.
One successful approach involves using X-ray computerised tomography (CT) scans of flow in layered rocks, in conjunction with detailed numerical simulations. The recent study by Jackson et al. (2018) presents a systematic approach to estimate the global effect of rock layers on the flow. A set of CT scan experiments is first performed in the viscous limit at high capillary number to determine the intrinsic properties of the flow, such as the relative permeabilities and capillary pressure (which are both typically functions of the saturation). Then, a similar set of experiments is performed at low capillary number to characterise the heterogeneity of the rock by means of fitting a set of capillary pressure scaling factors (one for every scanned voxel) to match numerical simulations to the CT scans. Having performed this two-stage analysis, Jackson et al. then use the fitted numerical model to describe the flow at intermediate capillary numbers, thereby enabling a systematic upscaling of the heterogeneities. In this way, relationships for the equivalent properties of the flow are derived, such as equivalent relative permeability, which are particularly useful when employed in conjunction with flow simulators to make predictions in the field. However, without being able to perform CT scans of flow in the rock samples, such analysis is impossible. Furthermore, there exists no general upscaled theory for the flow regime in between the viscous and capillary limits.
The objectives of the current study are to develop a simple theoretical tool that can be used to upscale the effect of heterogeneities in arbitrary flow conditions, where the heterogeneity can be given as a model input. The ultimate goal is to be able to study a vast range of scenarios to provide ensemble forecasts for the migration of immiscible fluids in porous media. Hence, this tool needs to be computationally inexpensive, and needs to be able to predict where and when heterogeneities affect the flow in the aquifer via the transition between viscous and capillary limiting regimes. Such a tool can be used not only to pinpoint optimal sites and predict trapping efficiencies for CO 2 sequestration, for example, but also for inverse modelling of flow properties given field measurements.
In the present study, we restrict our attention to a layered porous medium, with heterogeneity varying in the vertical direction and flow driven in the horizontal direction only. Furthermore, we focus on drainage flows, where a non-wetting phase drives out a wetting phase, although the analysis can easily be extended to imbibition flows. Using a combination of asymptotic analysis and numerical simulations of steady-state flow conditions, similar to Ekrann & Aasen (2000), we derive semi-analytical relationships for the equivalent relative permeabilities that are valid across all capillary numbers and saturations. We then use the upscaled properties to describe the dynamic flooding of an aquifer with small-scale heterogeneities. The latter is an extension of the classic model of Buckley & Leverett (1942), where a one-dimensional system is used to model the displacement of immiscible fluids in a long, thin porous medium.
Whilst there are numerous papers on the Buckley-Leverett problem and its variants (McWhorter & Sunada 1990;Schmid & Geiger 2012;Deng & King 2015;Zheng & Neufeld 2019), there are none which address the transition between the viscous and capillary limits in the case of a heterogeneous medium. In the present study, we use our simplified semi-analytical expressions to address the dynamics of this transition, showing that regions of the aquifer near the injection point (or at early times) lie within the viscous limit, whereas regions far away from the injection point (or at late times) lie within the capillary limit. Finally, we use this approach to quantify the effect of heterogeneities on the injection of CO 2 , making comparisons with field data from the Salt Creek case study (Bickle et al. 2017).
Section 2 describes the heterogeneous system we consider, and derives relationships for the upscaled flow properties in the viscous and capillary limits. In the case of intermediate capillary numbers, numerical simulations are used to characterise the viscous-capillary transition and semi-analytical composite expressions for the upscaled properties are proposed. Then, § 3 uses the upscaled flow properties to study flooding dynamics via the Buckley-Leverett problem, extended to heterogeneous media. In § 4 we compare our upscaling predictions with the experimental measurements of other authors, and finally we close by summarising the results.

Upscaling heterogeneities
The general approach taken here is as follows: we start by summarising the governing equations and boundary conditions for two-phase flow in a layered porous medium; next, we define upscaled quantities, such as the equivalent relative permeabilities; then we derive expressions for these upscaled quantities in each of the two limiting viscous and capillary cases, using some simple examples for illustration; finally, we use numerical p Figure 1. Schematic diagram of a long, thin two-dimensional aquifer with steady, pressure-driven flow of wetting and non-wetting phases. Vertical heterogeneity is given by variation in the pore entry pressure p e (z) and permeability k(z), which here is illustrated in the case of a two-layered system.
simulations to calculate the upscaled quantities for intermediate capillary numbers, showing how to incorporate all regimes using simple semi-analytical parameterisations.

Immiscible two-phase flow in porous media
We consider the flow of a non-wetting phase driving out a wetting phase (e.g. carbon dioxide driving out water) in a two-dimensional aquifer of length L, height H, and whose intrinsic properties (e.g. porosity φ, permeability k, pore entry pressure p e ) vary in the vertical direction z (see figure 1). We model the flow behaviour at the continuum scale (but below the scale of the heterogeneities) using conservation of mass and the multiphase extension to Darcy's law under gravity (Bear 2013). Hence, the governing equations for the flow are where subscripts n and w indicate non-wetting and wetting phases, and we require the fluids to fill the pore spaces, such that S n + S w = 1. The parameters μ i and ρ i are the viscosities and densities of either phase, g is the gravitational acceleration constant, k ri (S i ) are the relative permeabilities and p i are the pressures of each phase, which differ by an amount where p c is the capillary pressure associated with the micro-scale capillary forces between phases. Although k ri and p c depend on many factors in general, they are often assumed to be functions of the saturation alone (Golding et al. 2011). A simple, commonly used empirical relationship for the capillary pressure is that proposed by Brooks & Corey (1964), where p e (z) is the pore entry pressure, λ 1 represents the pore size distribution and is the rescaled saturation, such that s ∈ [0, 1]. The irreducible wetting phase saturation S wi represents the amount of wetting phase that cannot be removed, and is therefore always trapped in the pores by capillary forces. The pore entry pressure p e describes the minimum pressure required to allow any non-wetting phase into the pore spaces. For p n − p w = p e , only the largest pore spaces are filled with non-wetting phase, and for p n − p w > p e , smaller and smaller pore sizes are invaded. Clearly, the pore entry pressure depends on the porosity and geometry of the pores, as does the permeability, and we assume these vary in the vertical direction. Therefore, in this study, heterogeneities are defined solely by φ(z), p e (z) and k(z). It is often assumed that p e (z) and k(z) depend on the porosity under some power law that reflects the geometry of the pore spaces (Leverett 1941). Hence, we have p e ∝ φ −a , k ∝ φ b , for parameters a, b. We therefore take the pore entry pressure and permeability to be related according to where p e 0 and k 0 are typical dimensional scalings, and B = a/b > 0 is a positive constant, since larger pore spaces should correspond to lower pore entry pressure. It has long been argued that such power law relationships do not apply generally (Cloud 1941), but specific power laws are often used for particular rock types (e.g. see Nelson 1994). For example, using b = 2 and the scaling proposed by Leverett (1941), where p e ∼ (φ/k) 1/2 , gives a value of B = 1/4. There are a vast number of different empirical relationships for the relative permeabilities k ri which have been proposed by various authors (Krevor et al. 2012), and the appropriate choice depends on the specific rock type and fluid phases. The relative permeabilities are monotonic functions of their respective phase saturations, and lie between 0 and 1. In the limiting case where the flow becomes single phase, the relative permeability of that phase should be 1 (and 0 for the other phase). But as we have already discussed, there may be an irreducible wetting phase saturation, and hence we have k rn (s = 1) = k rn0 , for some 0 k rn0 1. In this paper, we propose a general framework which is not limited by a specific choice of empirical relationship. However, we make comparisons with several commonly used laws, including those proposed by Corey (1954) and Chierici (1984), which we give explicitly in appendix A.
Finally, to complete the model, we require a set of boundary conditions. There are many possible choices of boundary conditions for such flows, as discussed by Krause (2012). We note that after some simple rearranging, it is possible to convert (2.1)-(2.3) to equations for the pressure and saturation of one of the phases only. Therefore, without loss of generality, we formulate our model focussing on the non-wetting phase, and we consider a pressure driven flow, with boundary condition Note that (2.7) determines the flow rate of non-wetting phase at the inlet. Similarly, (2.8) and (2.9) determine the flow rate of the wetting phase (or equivalently the pressure drop of the wetting phase). Hence, it is often useful to replace (2.7)-(2.9) by flow conditions where the inflow parameters U n , U w are related to s i and p by the multiphase flow model. To summarise, the model consists of the governing equations (2.1)-(2.3), as well as boundary conditions (2.7)-(2.11), and some initial conditions for p n and s. The heterogeneity is characterised by φ(z), k(z), and p e (z), which are related by (2.6).

Upscaling
As discussed by numerous authors (Krause & Benson 2015;Rabinovich et al. 2016), heterogeneities have the capability of changing the overall flow properties of porous media. In particular, in the presence of heterogeneities the empirical relative permeability relationships discussed earlier tend to become wholly inaccurate as we deviate away from the classic homogeneous or viscous limiting case. Typically, parallel layering (as studied here) tends to segregate phases in such a way as to increase the overall flow of non-wetting phase, and decrease the flow of wetting phase (Krause & Benson 2015). For this reason, and as a method of reducing the requirement to resolve individual heterogeneities, it is useful to define so-called equivalent properties instead which give a description of the flow that upscales the effects of these heterogeneities.
For the purpose of upscaling, we restrict our attention to the steady-state case. Therefore, similarly to Jackson et al. (2018), we define the equivalent relative permeabilities as where the pressure changes p i refer to the difference between inlet and outlet for each respective phase, and the operator · refers to a type of spatial averaging, which we leave in general terms for now but discuss later in § § 2.4, 2.5 and 2.7. Similarly, we define the equivalent capillary pressure as p c eq = p c p e , (2.14) which is a dimensionless quantity. As discussed earlier, the effect of heterogeneities is often characterised by the so-called capillary number N c (1.1), which is given as the ratio between typical horizontal pressure gradients, and typical vertical gradients in the pore entry pressure. For the horizontal pressure change in (1.1), we choose the constant non-wetting pressure difference (2.7), although we could equally choose the wetting pressure, or some kind of combination. As we will discuss later, this choice is satisfactory The equivalent properties (2.13)-(2.14), which are the main focus of this paper, depend on the following different quantities: (i) The underlying heterogeneity of the rock, characterised by p e (z) and k(z) via (2.6).
(ii) The flow-driving pressure drop across the aquifer p.
(iii) The aspect ratio of the domain δ.
(iv) The inlet conditions of the saturation s i .
Note, the capillary number N c contains all of (i)-(iii), but has no notion of (iv). Furthermore, it does not describe the spatial variation of the heterogeneity, only the typical variation scale p e . In addition, the definition of N c depends on the choice of length scales H and L, which are not necessarily well defined in real applications. Therefore, even though N c is not sufficient on its own to characterise the complete flow picture, we use it primarily as a metric for describing the type of flow regime (horizontal pressure-driven flow versus vertical capillary-driven flow), a task for which it is well suited.

Non-dimensionalisation and asymptotic analysis
Before we address each of the viscous and capillary limits it is useful to convert to dimensionless variables. Let us attribute the following scalings to each variable where δ = H/L is the aspect ratio, which we assume to be small, and w i is the vertical velocity component of each phase. Written in terms of these new non-dimensional variables, the governing equations (2.1)-(2.3) (in the steady state) becomê ∇ ·û n = 0, (2.17)

20)
where we have introduced the non-dimensional variables M = μ w /μ n (mobility ratio), σ P = p e /p e 0 , ψ i = ρ i gH/ p, andÑ c = p/ p e = N c /δ is the reduced capillary number. For this study, we restrict our attention to thin aquifers ψ i 1, in which gravity can be neglected, similarly to the core flooding experiments of Jackson et al. (2018). The boundary conditions (2.7)-(2.11) becomê p n |ˆx =0 −p n |ˆx =1 = 1, (2.24) Likewise, the inflow of each phase is given bŷ where we have introduced the two non-dimensional flow parameters which represent the flow of non-wetting phase and the flow fraction, respectively. Finally, the power law describing the scaling between permeability and pore entry pressure, (2.6), becomes 1 + σ Ppe =k −B . (2.33) We choose the dimensional scaling k 0 as the vertical average of the permeability, such that k averages to unity but note that 1 + σ Ppe may not.

Capillary limit
To find solutions in the capillary limit, we consider an asymptotic expansion in the scaled capillary numberÑ c 1. We assume that the statistical properties of the heterogeneity are fixed, such that σ P remains O(1) (i.e. we consider a weak overarching pressure gradient that is independent of the rock properties). In addition, we restrict our attention to the case where the aspect ratio is much smaller than the flow perturbation, such that δ Ñ c 1. From the capillary pressure equation (2.23), it is clear that both wetting and non-wetting pressure should scale likep i ∼ 1/Ñ c . Therefore, the variables s,p n andp w are expanded inÑ c as s = s 0 +Ñ c s 1 + . . . , (2.34) p n =Ñ −1 cp n −1 +p n 0 + . . . , (2.35) Hence, (2.19)-(2.22) indicate that the pressures in both phases must be constant to leading order, such thatp n −1 −p w −1 = γ , for some value of γ . This is consistent with the definition of capillary limit given by other authors (Ekrann & Aasen 2000;Rabinovich et al. 2016). From (2.23) we therefore derive a leading-order expression for the saturation where we writeP e = 1 + σ Ppe for convenience. Given the form of (2.13) and (2.14), we would like to express (2.37) in terms of the averaged saturation. Since, to leading order, the capillary limit solution only depends onẑ, we select our averaging operator here as the vertical average · = 1 0 · dẑ, which we henceforth represent with an overline. In this way, (2.37) becomes (1 −s). (2.38) Note that the solution (2.38) also satisfies the outlet condition (2.26) and the impermeability condition (2.28). The inlet condition (2.25) is not satisfied, which will lead to a boundary layer over which the saturation transitions to the outlet state, as we discuss later.
To calculate the equivalent relative permeabilities (2.13), we first need the averaged Darcy velocities, which only appear at first order. These are obtained by vertically integrating (2.19), (2.21) and using (2.29), (2.30), to give (2.40) By integrating (2.39) and (2.40) along the channel length, we arrive at expressions for the total changes in pressure along the channel, which we then insert into (2.13) to finally arrive at the capillary limit for the equivalent relative permeabilities The expressions (2.41) and (2.42) are a generalisation of the arithmetic mean expressions derived by Rabinovich et al. (2016) in the case where the heterogeneity consists of a set of horizontal layers. The equivalent capillary pressure is found by inserting (2.38) into (2.14), giving It should be noted that the capillary limit solution (2.38) may lead to negative saturation values fors which is clearly unphysical. In such situations, the saturation profile is instead given by and consequently there are regions of space devoid of non-wetting phase, a phenomenon associated with very strong heterogeneities. In this case, it is less straightforward to relate the capillary pressure constant γ to the mean saturation analytically. However, a nonlinear relationship can be established numerically instead. Note that we could go to higher order in the asymptotic expansions to capture near capillary limit behaviour. However, for the purposes of understanding the leading-order impact of capillary heterogeneity on the flow, we find leading-order solutions sufficient.
2.5. Viscous limit In contrast to the capillary limit, the viscous limit relates to the regime where the flow-driving pressure gradient is much larger than the capillary forces, such that the heterogeneities in capillary pressure do not affect the flow. Therefore, to address this limit we consider a small capillary correction p e / p =Ñ −1 c 1. Note that the pore entry pressure is related to the scaled capillary number via the parameter σ P = CÑ −1 c , where C = p/p e 0 . For this analysis, we assume that the overarching pressure gradient is fixed, such that C remains O(1) (i.e. we consider a weak heterogeneity p e independently of the pressure gradient). Furthermore, we assume that the aspect ratio is much smaller than the heterogeneity perturbation, such that δ Ñ −1 c 1. Given the power law relationship (2.33), we also havek Similarly to the capillary limit, here we seek an asymptotic solution, except now this is given in terms of powers ofÑ −1 c , such that s = s 0 +Ñ −1 c s 1 + . . . , (2.47) In this way, (2.19) and (2.21) indicate that there are no leading-order vertical pressure gradients ∂p n 0 /∂ẑ = ∂p w 0 /∂ẑ = 0. Furthermore, (2.23) indicates that, to leading order, which implies that s 0 must also be independent ofẑ. This also ensures that the impermeability condition (2.28) is satisfied at leading order. The Darcy velocities are obtained by vertically integrating the system (2.17)-(2.23) and using (2.29), (2.30), to give (2.52) Due to (2.52), the zero gradient boundary condition (2.26) can only be satisfied if s 0 is constant. This is equivalent to the condition which enforces a relationship between the flow fraction f 0 and the saturation s 0 . Therefore, since the viscous limit solution is constant to leading order, the averaging operator in (2.13) and (2.14) is trivial. With this taken into account, the viscous limit expressions for the equivalent relative permeabilities are Furthermore, the equivalent capillary pressure is given by The viscous limit expressions (2.54)-(2.56) are identical to the original expressions for relative permeability and capillary pressure, which is expected in the limit of vanishing heterogeneity. Note that this analysis can be extended to higher-order terms to approximate the case of a large but finite capillary number. However, we find a leading-order analysis satisfactory for our purposes.
2.6. Types of heterogeneity Whilst the above analysis applies for any given vertical heterogeneity and empirical relative permeability relationships k rn , k rw , we shall now discuss how our predictions manifest in an example scenario. We choose a simple background heterogeneity which consists of a sinusoidal perturbation on a uniform permeability profilê for some amplitude A and wavenumber n ∈ N. Meanwhile, the pore entry pressure is given by (2.33), in terms of some power B. For the intrinsic relative permeabilities k ri , we use the classic empirical power law of Corey (1954), which is given by (A1) and (A2), with a quadratic power law. A full list of parameter values is found in appendix A. Although in reality more complex permeability profiles may be present than a sinusoidal variation, we use (2.57) because, as a canonical function, it illustrates the fundamental effects of amplitude and wavelength on the flow properties. Furthermore, any sufficiently smooth and continuous permeability profile can always be decomposed into a Fourier series of such modes. Nevertheless, we have also investigated other permeability profiles, such as a layered profile which we illustrate in figure 12 in appendix B.
In figure 2 we plot the viscous limit (which is independent of heterogeneity) and the capillary limit for different values of A and B (for a fixed value of n = 1). The plots confirm that heterogeneity has the overall effect of lowering the flow of the wetting phase, and raising the flow of non-wetting phase. This can be explained by (2.38), which indicates that s is larger in places where the pore entry pressure is smaller, and hence in regions of larger pore space. Hence, capillary pressure forces the non-wetting saturation to preferentially segregate to regions of larger space, where it is easier to flow. Increasing the amplitude A accentuates this effect, since this corresponds to stronger heterogeneity. It is also accentuated by increasing the power law B, since this increases the strength of the pore entry pressure heterogeneity.
Note in some cases it is possible to derive analytical formulae for the equivalent relative permeabilities in the capillary limit. For example, in the simple case where B = 1, Viscous and capillary limits of equivalent relative permeability (2.13) (note the non-wetting relative permeability is normalised by k rn 0 = 0.116) for a sinusoidal heterogeneity (2.57) and a power law relationship for the pore entry pressure (2.33). The capillary limit is shown for different values of the heterogeneity amplitude A (fixing B = 1/2) (a) and power law B (fixing A = 0.8) (b). Experimental data taken from Bennion & Bachu (2005) in the viscous limit. (c,d) Grey scale maps of the percentage difference between viscous and capillary limit predictions for a heterogeneity with two wavenumbers n 1 , n 2 (2.60).
the resulting expressions are (2.59) The expressions (2.58) and (2.59) are valid for amplitudes A < 1, although only for values ofs large enough so that (2.38) does not have s = 0 anywhere (or according to (2.44), for In situations where there are regions of zero saturation, an analytical formula is still possible, though the expressions are more complicated so we do not display them here. In contrast to A and B, varying the wavenumber of the perturbation n ∈ N does not have a significant effect on k rn cap , k rw cap . However, more interesting effects are observed when two different wavelengths are introduced, such that the permeabilitŷ where the factor F is chosen such that the difference between the maximum and minimum perturbation (and hence the capillary number) is kept the same. In figure 2(c,d) we display grey scale plots of the percentage difference in equivalent relative permeability between the viscous and capillary limits, for different values of n 1 and n 2 . Since the plots are symmetric about n 1 ↔ n 2 , we only display half of the phase space. Clearly, the maximum difference occurs when n 2 = n 1 (at constant values of 55 % and 31 %), but there are also streaks near n 2 = n 1 /3, n 2 = n 1 /2, n 2 = n 1 /4, and so on (in descending order of magnitude). Whilst these heterogeneities are idealised, this simple investigation serves as an illustration for the different types of permeability and pore entry pressure one might encounter in the field. In particular, we have indicated how upscaled quantities depend on model parameters in the two limiting viscous and capillary limits, which will be useful throughout the paper. Next, we move on to model situations which are not in either of these two limits, but instead lie somewhere in between.

Intermediate capillary number
In the case of intermediate capillary number, there are two possible approaches: either we can perform numerical simulations of steady Darcy flow (2.17)-(2.23) with boundary conditions (2.24)-(2.28) and then calculate the equivalent properties (2.13); or we can go to higher-order terms in the asymptotic expansion of each of the viscous limit or the capillary limit. We prefer to use the numerical approach here, similarly to Virnovsky et al. (2004), since it gives a complete description that is valid across all capillary numbers, and this is more convenient than patching together asymptotic solutions from different regimes. However, in contrast to Virnovsky et al., we use our numerical simulations together with our previous analytical results to form composite expressions for the equivalent properties which can be readily applied elsewhere.
Although the previous analysis related to the scaled capillary numberÑ c , here, we keep everything in terms of the original capillary number, N c , since this is more common in the literature, and therefore makes our results more accessible.
We have calculated numerical solutions for capillary numbers N c between 1 and 10 4 and a heterogeneity (2.57) with amplitude A = 0.6 and power law B = 1/2. In addition, we set the aspect ratio as δ = 0.1. The numerical solutions are calculated using a fourth-order central difference scheme in space (with 80 × 20 grid points in the (x, z) directions) and a pseudo-time-stepping method that converges iteratively. We use the method of continuation to advance quickly through several orders of magnitude of the capillary number.
In figure 3(a) we display colour plots of both the wetting and non-wetting saturations, overlaid with streamlines given by the Darcy velocitiesû i for three different values of the capillary number. For small capillary numbers, the flow segregates into two separate streams, where all the non-wetting phase moves to the more permeable regions, and vice versa. There is a small region of strong transverse flow of wetting phase near the inlet due to sharp saturation gradients. For larger capillary numbers, the saturation profile is more uniform throughout. The segregation of phases is less pronounced, and there is little transverse flow near the inlet.
There is a kind of horizontal boundary layer in saturation distribution that exists near the inlet, over which the saturation transitions from the constant inflow value s i to the capillary limit solution downstream. The boundary layer thickness, which we denote δ visc , grows with capillary number. By defining δ visc as the distance needed to reach the capillary limit solution (2.38) to 90 % accuracy, we can plot the variation with capillary number, as can be seen in figure 3(b). Hence, we find that the boundary layer thickness δ visc is approximately proportional to N 3/5 c . Note that, if we were to extend the aquifer sufficiently, all cases would eventually reach the capillary limit. This is evident by noticing that the only solution to (2.17) and (2.28)  which is independent ofx is the capillary limit solution (p c = constant). Therefore, in the transition between the viscous and capillary limits, the inlet condition s i is of critical importance. Indeed, if we were to choose the inlet profile as (2.38), then any capillary number would result in the capillary limit solution. To mitigate this, we have chosen s i as a constant value so that both viscous and capillary limits can be recovered in the limit of large and small capillary numbers, respectively. In addition to the capillary number, the boundary layer thickness must clearly depend on the aspect ratio δ, and we have plotted this dependence in figure 3(c), holding the capillary number fixed at N c = 8. In this case, we see that δ visc grows linearly with aspect ratio. This is expected due to a uniform stretching of the domain. Clearly, the choice of the domain dimensions for upscaling has a significant impact on the resulting upscaled quantities, presenting a challenge for creating a general theory of upscaling. Later in § 4.3 we discuss how varying the choice of domain size may affect predictions.
To calculate equivalent properties of the flow, it is necessary to choose an appropriate averaging operator · in (2.13) and (2.14). We are dissuaded from choosing a core average, since undesirable boundary layer effects from the inlet make it impossible to recover the capillary limit solution, (2.41) and (2.42), as we decrease N c . Instead, we find the most convenient choice is a vertical average at the aquifer outlet · = 1 0 · dẑ|ˆx =1 . Since we have chosen zero gradient conditions (2.9), this removes boundary effects from the averaging process as much as possible. In the case of the pressure drop in (2.13), we use an average of the non-dimensional pressure gradient p i = ∂p i /∂x. Using this averaging method allows the solution to converge to both capillary and viscous limit solutions consistently.
The equivalent relative permeabilities and capillary pressure are shown in figure 4(a,b). The points on each coloured line have the same capillary number and different values of the inlet saturation s i (or equivalently the flow fraction f 0 = U w /U n ). In this way, it is possible to observe how the equivalent relative permeabilities vary over both saturation and capillary number, as illustrated in figure 4(c,d). As indicated in the plots, the equivalent relative permeabilities are very well approximated by the transition function with parameter values N c t = 394, Δ = 5.5 and k ri ± = k ri visc ± k ri cap , where the viscous and capillary limits are given by (2.41), (2.42), (2.54) and (2.55). The composite expression (2.61) captures the numerical results with mean relative error of around ∼1 %. Although an even better fit can be attained by allowing N c t and Δ to vary with saturation and capillary number, we take them as constants here for the sake of simplicity. The transition capillary number N c t represents the capillary number that lies logarithmically as a midpoint between the viscous and capillary regimes. The parameter Δ represents one logarithmic folding scale. As we can see in figure 4(c,d), the viscous and capillary limits are little more than one folding scale away from the transition capillary number on either side. These two parameters N c t and δ fully characterise the flow regime for intermediate capillary numbers, and they are subtly related to the boundary layer thickness discussed earlier. Hence, they are not universal for every scenario, since we have shown that the boundary layer thickness depends on the choice of domain aspect ratio and inlet conditions s i . Therefore, great care must be taken when choosing the domain for upscaling, as we discuss later in § 4.3.
Apart from the sinusoidal permeability profile investigated here, we have also tried numerous other types of permeability profiles (e.g. step function, Gaussian, . . . ) and different power law values B, and in each case (2.61) gave good comparison with the numerics, indicating the robustness of our current approach. For example, in figure 12 in appendix B we display equivalent relative permeabilities for a step-layered profile, together with the corresponding fit (2.61), showing close agreement.
Note that we could have equally fit the data to the capillary number defined in terms of the wetting pressure change instead of the non-wetting pressure change (see (1.1)). However, we observe that the ratio of these pressure changes is (2.62) Hence, the two definitions are not independent, and would just result in a different form of (2.61). Therefore, without loss of generality, we keep the capillary number defined in terms of non-wetting pressure difference. Variation in the equivalent capillary pressure (2.14) is much less significant, since p c cap /p c visc = 1.06. This can be seen in figure 4(b), where the capillary and viscous limit curves lie almost on top of each other. Therefore, there is not a great need to model the transition behaviour, and it is sufficient to assume the viscous limit everywhere In the next part of the study, we use the equivalent properties derived here to study dynamic flooding in an aquifer.

The Buckley-Leverett problem for heterogeneous media
3.1. Problem summary Now that we have analytical expressions for the equivalent relative permeabilities in the viscous and capillary limits (2.41), (2.42), (2.54), (2.55), and a composite expression (2.61) for intermediate capillary numbers fitted against numerical data, we have a full description of the equivalent properties across all flow conditions. Next, following the classic study by Buckley & Leverett (1942) of the displacement of immiscible flows in a long, thin aquifer, we extend this to the case of heterogeneous media, using our upscaled equivalent properties.
In the classic Buckley-Leverett problem, a one-dimensional porous medium, initially filled with saturation s ∞ , is flooded with a saturation s i at the inlet x = 0 (see figure 5a). While the problem is time dependent, we make the key assumption that the equivalent properties derived in § 2 still apply even when the flow is unsteady, which follows the approach taken in industrial applications. Our analysis here can be interpreted as the macroscopic flow picture of an aquifer with an underlying heterogeneity, where the length scale of the heterogeneity is much smaller than the flow length scale (see figure 5c).
A complete discussion of the Buckley-Leverett problem can be found in any standard porous media textbook, such as Bear (2013) and Woods (2015) for example. Here, we simply summarise the problem and describe how it can be extended to heterogeneous media. In the original problem formulation (for homogeneous media), the governing dimensional equation for the saturation is where the advective and diffusive terms are given by which can be derived by combining (2.1) and (2.2), where V tot = u n + u w is the total Darcy flow (conserved). Note that we have rescaled time in (3.1) by a factor of φ(1 − S wi ) for convenience. To extend to heterogeneous media, we replace the relative permeabilities and capillary pressure in (3.2) and (3.3) by their equivalent counterparts derived earlier, and the saturation s is interpreted as an upscaled saturation. (Note that, in the case where relative permeability depends on the capillary number (2.61), the advective velocity (3.2) contains a partial derivative with respect to N c . However, due to the logarithmic dependence this contribution is very small (e.g. O(10 −9 )-O(10 −3 ) for typical parameter values) and so we ignore it here.) Hence, this extension to the Buckley-Leverett problem, although it is one-dimensional, contains information about the vertical variation of saturation in the rock and flow properties. Furthermore, the rock heterogeneities only manifest in these upscaled quantities and their typical scalings (φ 0 , p e 0 , k 0 ). In figure 5(d-f ) we plot the advective and diffusive components, given in non-dimensional termsV = VLμ w /k 0 p e 0 ,K = Kμ w /k 0 p e 0 , for both the capillary and viscous limits. We also plot the nonlinear Péclet number Pe =V/K. For the purposes of this comparison we define a non-dimensional flow rate

4)
and we use typical parameter values, giving U = 3167 and a viscosity ratio of M = 30. A full list of dimensional parameters is given in table 1 (taken from the Salt Creek case study, which we discuss later). Several observations can be made immediately. Firstly, for these typical parameter values the diffusive term is much smaller than the advective term (indicated by the Péclet number), such that the diffusive term can be neglected, except perhaps when saturation gradients are very large (e.g. for shock solutions Woods 2015), or when s is very close to 1. Secondly, the faster limit (between viscous and capillary) depends on the saturation value. Finally, the slight kink in the capillary limit advection velocity curve in figure 5(e) is due to non-smooth changes in the saturation distribution due to (2.45).
It is well known that the non-monotone behaviour of V can result in multi-valued saturation distributions, as illustrated in figure 5(b). This is often dealt with by introducing a shock at some intermediary saturation s s , where the saturation value is found by solving the equation in terms of the advective flux J = V ds and the initial saturation s ∞ . The shock equation (3.5) can be derived by a conservation of mass balance across the shock (Woods 2015). A typical shock solution is illustrated in figure 5(b), where the original multi-valued solution is overlaid as a dashed line. In reality, the steep saturation gradients present in such a shock solution would be softened by the diffusive term (3.3) over a growing length scale ∝ (t/Pe) 1/2 . For typical situations, this results in a diffusive boundary layer of approximately 1 %-5 % of the total aquifer length. The solution behaviour of the Buckley-Leverett problem is characterised by several saturation values: the inlet saturation s i , the initial far-field saturation s ∞ and, should a shock develop, the shock saturation s s . Since we restrict our attention to drainage flows (e.g. CO 2 driving out water), we confine our analysis to s i > s ∞ . To understand the different flow regimes, it is useful to introduce the stationary point saturation value s m , which corresponds to the saturation at which the maximum advection velocity is achieved (e.g. see figure 5e). A multivalued saturation profile never develops (i.e. no shocks) for parameter values s m s ∞ s i , as illustrated by a yellow region in the phase diagram in figure 6(d). Hence, in the absence of shocks, the flooding front moves at the far-field saturation speed, which is V = V(s ∞ ). Likewise, a shock will always develop for s ∞ s m s i , and the flooding front moves at the shock speed V = V(s s ).
We note that (3.5) may result in a shock saturation value that lies outside of the range [s ∞ , s i ]. Therefore, in such cases (3.5) is replaced by the condition s s = s i , such that the shock value is simply equal to the inlet value, as illustrated with dark blue colouring in figure 6(d).

Viscous and capillary limits
Now that we have summarised the Buckley-Leverett problem, the next step is to discuss the two limiting viscous and capillary cases. In figure 6(a,b) we display a colour plot of the front velocity values for each of these limits V visc , V cap (normalised by their maximum values) over all possible values of s i , s ∞ . In figure 6(c) we plot the ratio between these two  limits V visc /V cap . Wherever the far-field saturation is larger than the stationary point s ∞ > s m , viscous advection speeds dominate, whereas in regions with s ∞ near zero (leading to shocks), capillary advection speeds dominate. The maximum and minimum values of the speed ratio V visc /V cap are 1.44 and 0.13, indicating that neglecting heterogeneities may lead to substantial error in flooding predictions. For modelling carbon sequestration, where s ∞ is expected to be near zero (CO 2 is typically injected into brine-saturated aquifers), the implications are that in situations where the capillary number is small, heterogeneities cause an overall acceleration of the advancing front. This will play an important role in trapping mechanisms and storage efficiency.
To illustrate these findings, in figure 7 we display two solutions to the extended Buckley-Leverett problem with and without shocks. In the first case, figure 7(a,c,e), we flood an aquifer which is initially saturated with a substantial fraction of gas s ∞ = 0.5. In the second case 7(b,d, f ), the aquifer is initially saturated with the minimum possible gas amount s ∞ = 0 (see the markers in figure 6c), causing a shock to develop. In each case we plot the saturation s and velocity V at both the initial time, and at a single later time, indicating both capillary (red) and viscous (black) predictions. We display all results in non-dimensional form, where a suitable non-dimensional time scale is (3.6) The saturation profiles are obtained by solving the characteristic equation for each x value, such that No shocks (s ∞ = 0.5) where the saturation value s is conserved along characteristics. As initial conditions, we use a localised initial saturation distribution where x * /L = 10 −3 . In figure 7 this initial saturation profile is advected according to either the capillary or viscous limit speed, V cap (s) or V visc (s) (c,d).
In (e, f ) we also plot the position of the leading edge of the flood X(t), which increases linearly with time, with slope V = V(s ∞ ) or V(s s ). The speed ratio is V visc /V cap = 1.44 in the case without shocks, and V visc /V cap = 0.82 in the case with shocks. For applications such as CO 2 sequestration, this indicates that a model which neglects the effects of heterogeneities may predict flooding speeds with nearly 50 % inaccuracy. where we have used the definition in terms of the non-wetting pressure gradient. The local pressure gradients are given by

Transition between viscous and capillary limits
(3.10) We note that the capillary number used here (3.9) is defined differently to (1.1), which was used to perform steady-state upscaling earlier. However, (3.9) can be interpreted as the local capillary number for a macroscopic flow description, whereas (1.1) can be interpreted as the bulk capillary number for a small-scale study. Hence, the two definitions become equivalent by zooming in or out of the aquifer appropriately.
Since the pressure gradient (3.10), and consequently the capillary number, are both functions of s, they are conserved along characteristics. Hence, the capillary number at the flooding front x = X(t) is the same for all time (though different to the capillary number at the inlet x = 0, for example). Hence, the capillary number at a given saturation value is calculated for all time by solving the nonlinear implicit equation . (3.11) We summarise the steps for modelling the transition between the viscous and capillary limits as follows: first the capillary number must be calculated for some initial saturation data (e.g. (3.8)) using the implicit equation (3.11); then, if shocks are present, the shock saturation value s s must be calculated using (3.5), where the advection velocity V (3.2) and flux J use the composite expressions (2.61) for the equivalent relative permeabilities; if no shocks are present, the flooding front simply corresponds to saturation value s ∞ ; finally, the solution is calculated for all time via the characteristic equation (3.7), where V (3.2) contains the composite expressions (2.61). For example, using the same parameter values as in figure 7(b,d, f ), and setting N c t = 394, Δ = 5.5 in (2.61), as before, we calculate that a shock develops at saturation value s s = 0.38 (which lies in between the capillary and viscous limit shock values s s = 0.29 and s s = 0.47) at a capillary number of N c = 303. Then, the shock (which remains at this capillary number for all time) is advected at a constant velocity which is approximately 10 % faster than the viscous limit and 10 % slower than the capillary limit.
Interestingly, if one were to consider an axisymmetric flooding instead of two-dimensional plane flooding, the flow speed and pressure gradients would decay radially due to conservation of mass. Hence, the capillary number would also decay radially, such that different regions of the aquifer switch between viscous and capillary limits over time.
An axisymmetric model is more realistic than the two-dimensional case for cases where injection occurs at a single point source, as is often the case in industry. In the context of our current modelling approach heterogeneities are below the continuum scale, and consequently the equivalent relative permeabilities derived in § 2 can equally be used to describe a two-dimensional (as above) or axisymmetric setting. Hence, in the next section we extend the above analysis to axisymmetric flow. which can be re-written as d dt

Axisymmetric flooding
(3.14) The pressure gradients are given by .
( 3.15) which are no longer constant along characteristics (since (3.15) contains r), and so the capillary number (now defined in terms of ∂p n /∂r) must be calculated at each radial value. Hence, given some initial data for s, such as (3.8), the solution is found by time integrating the coupled system d dt where Q = Q tot μ w /k 0 p e 0 . In figure 8 we display solutions to (3.16) and (3.17) using the parameters s i = 1 and s ∞ = 0.35 (i.e. no shocks). In figure 8(a) we display the capillary number N c (r, t) at four times, which decays like ∼1/r as r → ∞. In figure 8(b) we display the position of the flooding front R(t) for each case, also indicating the capillary and viscous limit predictions for comparison. Unlike the two-dimensional case, here, the front moves like the square root of time (instead of linearly). Also, unlike the two-dimensional case, the capillary number at the the flooding front changes over time. At early times, the entire flow is close to the viscous limit, whereas at late times, nearly all the flow is close to the capillary limit, except for a small region near the origin. At intermediate times the flow straddles between the two limits. This can be seen in figure 8(b), where the front evolution switches between viscous-like behaviour to capillary-like behaviour over time.
We also display surface plots of the saturation at different times in figure 8(c). The colouring in each plot is chosen as a binary value depending on whether the local capillary number is above or below the transition value N c t (see also figure 8a, where one folding scale is illustrated). The result is that the flow near the source is in the viscous limit, and is consequently unaffected by heterogeneities. However, as the flood spreads through the aquifer the heterogeneities play a strong role far away from the origin. The overall effect is a deceleration, driven largely at the leading edge of the injection. Note that if we were to choose a smaller value of the far-field saturation, such as s ∞ = 0, a shock would develop and the advection speed in the capillary limit would be faster than that of the viscous limit (as in figure 7b,d, f ).
Similarly to the two-dimensional case, by neglecting the effects of heterogeneities, flooding speeds can be misrepresented by as much as 50 %. Therefore, for applications in CO 2 sequestration, modelling the transition of the flow between the viscous and capillary limits is critical for accurately predicting how far the injection has spread, and this is important both from safety and efficiency perspectives.

Comparisons with experimental data
In this section we compare some of our results to different sources of experimental data from other authors. Firstly, we compare the results of our steady-state upscaling from § 2 to some X-ray CT scan experiments. Then, we compare our dynamic predictions from § 3 to field measurements from a CO 2 injection experiment in Salt Creek, USA.

Steady-state upscaling
We now quantitatively compare our results to data taken from core flooding experiments. The recent study of Jackson et al. (2018) calculates equivalent relative permeabilities using X-ray CT scans of Bentheimer sandstone with parallel layers (Peksa, Wolf & Zitha 2015). Their analysis provides a three-dimensional map of the pore entry pressure in a rock core, a two-dimensional slice of which is illustrated in figure 9(a). To upscale the observed heterogeneities, the intrinsic relative permeabilities k ri were first approximated by fitting the empirical relationship proposed by Chierici (1984), which is given explicitly by (A3) and (A4), to CT scans at very high capillary number. Then, a set of experiments at very low capillary number was used to iteratively fit a numerical model of the core to experimentally observed saturation data. A full list of the parameter values is given in appendix A.
Unlike their three-dimensional data, heterogeneities discussed here depend on the vertical dimension alone. Therefore, we take an average of the experimental data, p e (z) = p e exp (x,ŷ,ẑ) dx dŷ, which is illustrated in figure 9(b). Evidently, the experimental rock core has some longitudinal variation, so we do not expect our comparison to be perfect. Indeed, the mean relative error in approximating the pore entry pressure data in figure 9(a) as the simplified transverse/longitudinal average in figure 9(b) is ∼15 %. However, a good approximation should be attained, since the layering is predominantly parallel to the flow.
To compare with these experiments, we start with the two viscous and capillary limiting cases, since all other cases must lie between these. The capillary and viscous limits derived  figure 9(c). Spatial variation in the permeability k is not provided, so we fit our power law relationship (2.6) against their capillary limit data, giving B = 1/10, with a mean relative error of 23 %, which is most likely attributed to our approximation of heterogeneity by a simple vertical variation. For each of the pore entry pressure and permeability, we calculate the standard deviation divided by the mean, giving σ ( p e )/μ( p e ) = 0.16 (which is the same as quoted by Jackson et al.) and σ (k)/μ(k) = 0.74 (which is similar to field observations from Salt Creek, discussed later). The next step is to compare equivalent relative permeabilities for intermediate capillary numbers. To do so, we use our numerical simulations, as described earlier. Our calculated equivalent relative permeabilities are shown in figure 9(d), compared against the data of Jackson et al. (2018) in figure 9(c). Each coloured line on the plot has the same value of the total Darcy flow U tot = U n + U w and different values of the flow fraction f 0 = U w /U n . Consequently, the capillary number varies greatly over one value of U tot and so, following Jackson et al., we quote the value at f 0 = 0.5. To ensure that the quoted capillary numbers are the same, we use the same definition as Jackson et al. for the capillary number, where the pressure change in (1.1) is over the whole core. Overall, the comparison is good, with our data points varying between the viscous and capillary limits in a similar manner to Jackson et al. However, the slight differences in the curve shapes are most likely attributed to our one-dimensional approximation of the heterogeneity (which is only accurate to around ∼15 %, as described earlier).
For some core samples, the heterogeneities are not aligned with the flow direction at all, such as the Bunter sandstone analysed by Jackson et al. (2018), where perpendicular layering is present. In such situations, the details of the upscaling described here for flow along parallel layers may not be appropriate. However, the continual arrangement/rearrangement of the saturation as it moves downstream is similar to the flow rearrangement observed in the boundary layer described earlier. This can therefore be interpreted as an effective raising of the capillary number, since the flow is always close to the viscous limit by virtue of the fact that it does not have the spatial freedom to align with the capillary limit easily.  Krevor et al. (2012) for Paaratte sandstone in the viscous limit). (d) Upscaled predictions of the volume fraction of CO 2 at the observation well (4.1), compared with field measurements. The CO 2 volume fraction of the produced fluids (red solid curve) is calculated from the temperature (e), assuming adiabatic cooling, given the variation of density and coefficient of thermal expansion of CO 2 with pressure and temperature from Dubacq, Bickle & Evans (2013) and specific heats of CO 2 and water from Holland & Powell (2011). Temperature drops at days 15, 47-48 and 143-144 are related to reductions in production rates. High reported volumes of produced CO 2 between days 107-113 do not coincide with any changes in production rate or temperature fluctuations and are disregarded.

Dynamic flooding
To compare our extension to the Buckley-Leverett problem for heterogeneous media to field data, we use the Salt Creek CO 2 injection experiments from 2010, as detailed by Bickle et al. (2017). CO 2 was injected into a sandstone aquifer with vertical permeability structure as shown in figure 10(a), and aspect ratio δ ≈ 25 m/200 m. Injection was performed in several rows of wells, so that a two-dimensional model is probably more accurate than a radially symmetric one. Variations in the topography are neglected.
Relative permeability curves are not available for this sandstone, so to model this case study we use the curves of a similar sandstone called the Paaratte formation located in SE Australia, as detailed by Krevor et al. (2012). We display the empirical relationships (A5) and (A6) in appendix A. Likewise, pore entry pressure variation is not available, so we try using several different values of the power law B (2.6). We display the equivalent relative permeability curves for both the viscous limit, and the capillary limit (for several different values of B) in figure 10(c). Power laws 1/20 B 1/10 seem to give reasonable results. Moreover, for these B values, the value of the ratio between the pore entry pressure standard deviation and mean is σ ( p e )/μ( p e ) ∈ [0.1, 0.2], as compared to the Bentheimer sandstone of Jackson et al. (2018) which has σ ( p e )/μ( p e ) = 0.16. For such pore entry pressure distributions, we display the corresponding capillary limit saturation distributions (2.38) in figure 10(b). For comparison with field data, we use a mid-range value of B = 1/15. Following our extension to the Buckley-Leverett problem (ignoring diffusion), we use (3.1) to describe the temporal evolution of an injection of CO 2 , with s i = 1, s ∞ = 0.
We use (2.61) for the equivalent relative permeabilities with N c t = 394 and Δ = 5.5, as before. We choose a driving flow of V tot = 1.6 × 10 −6 m s −1 which results in a pressure drop across the aquifer between 4-8 MPa, which is consistent with field measurements. The full list of parameter values for this problem is given in appendix B.
Using all of the above information, we can compare our model predictions to field measurements. One useful metric for comparison is the volume fraction of CO 2 at the observation well, for which field data are available. The predicted volume fraction of CO 2 at any given saturation value and capillary number is Mk rn eq (s, N c ) + k rw eq (s, N c ) , (4.1) which we calculate at observation well 28WC2NW05 (200 m from the injection well) and plot in figure 10(d) (dotted blue curve). Due to the small far-field saturation value, a shock develops, creating a sharp advection front which moves at constant velocity through the aquifer, such that arrival at the observation well manifests as a discontinuous jump in CO 2 volume fraction. The diffusion term (3.3) which we neglected would smooth out the saturation profile near the shock in a diffusive boundary layer of growing width ∝ (t/Pe) 1/2 . However, since the Péclet number for this flow is so large, this manifests in a very small error margin, as illustrated with blue shading in figure 10(d).
In figure 10(d) we compare these predictions to field measurements of the volume fraction of CO 2 in the produced fluids. We consider that the volume fraction given at reservoir temperature and pressure (red solid curve), which is calculated from the temperature (figure 10e) of the produced fluids (assuming adiabatic cooling), is more reliable than the reported CO 2 production based on spot measurements (black curve).
Our modelling predicts breakthrough of CO 2 at volume fraction J ≈ 75 %, after 66 days, whereas the observations suggest significant CO 2 (J ≈ 10-20 %) arriving between 65 and 86 days after the start of injection. The breakthrough times for the capillary and viscous limits (which we plot in figure 13 in appendix B) are 50 and 83 days, indicating a significant effect of heterogeneities.
It should be noted that, whilst the field measurements only detected significant CO 2 breakthrough after ∼65 days, small quantities of noble gas tracers ( 3 He and 129 Xe) added to the CO 2 stream at the start of injection were detected only 10 days later. This suggests that regions of the aquifer, such as the high permeability zone at mid-depth, may advect CO 2 at much greater velocity than the bulk. This would also explain why the field data have a much lower, more spread out volume fraction than our predicted curve. Therefore, this motivates a slightly more resolved upscaled model that breaks the aquifer up into smaller regions. We discuss this and other questions regarding the choice of length scales in the next section.

A note on the choice of length scales
One of the key difficulties, and still an open question in the process of upscaling, is the choice of length scales. We demonstrated this earlier in figure 3 by showing that the boundary layer thickness depends on the aspect ratio of the upscaling domain, independently of the capillary number. Therefore, the viscous-capillary transition, characterised by the parameter N c t clearly depends on the domain over which upscaling is performed. However, as illustrated with regions (i)-(iii) in figure 11(a) for the Salt Creek permeability data, the aquifer can sometimes be naturally divided into subdomains. For the Salt Creek site, there is clearly a mid-depth region of very high permeability between regions of relatively low permeability. As the field data in figure 10(d) suggest, this may be After upscaling the heterogeneities within each of the three regions, the corresponding upscaled advection velocities V in each region are illustrated in (c,e). We also illustrate the standard upscaled velocities V cap , V visc , as well as the average velocity after upscaling the three regions independently (V cap ,V visc ).
responsible for a more distributed arrival of CO 2 at lower volume fraction than predicted by our upscaled model. Hence, it is not obvious whether it is more accurate to think of the aquifer as a single medium or three vertically stacked media, each to be upscaled separately.
In figure 11(b,d) we illustrate how the saturation of CO 2 would be distributed vertically in the aquifer in each of the capillary and viscous limits (for a mean value ofs = 0.2). The viscous limit has a uniform distribution, whereas the capillary limit is given by (2.38), leading to a focusing of CO 2 in the high-permeability region, and mean saturation values within each of the three subdomains ass = 0.082, 0.225 and 0.081. Now, if we upscale each of the three subdomains separately, we get three sets of equivalent flow properties k ri eq , and three different advection coefficientsV in the Buckley-Leverett problem. In figure 11(c,e) we illustrate how each of the three individual upscaled advection speeds would vary between subdomains, compared to the original viscous and capillary limits for the whole domain. The high-permeability region has a high-speed finger of CO 2 which precedes the low-permeability regions on either side, which is consistent with field observations (Bickle et al. 2017). Indeed this CO 2 finger may travel at almost double the speed of the bulk in the case of the capillary limit, and at almost quadruple the speed of the bulk in the viscous limit. In the viscous case, the mean advection speed of the three upscaled regionsV visc is equal to the upscaled speed of the whole region V visc , as expected. By contrast, this is not the case for the capillary limit, withV cap being approximately 65 % of the original upscaled advection speed V cap .
We compare this three-layered upscaling approach to the Salt Creek field data in figure 10(d) (blue dashed curve). Now, instead of a single bulk arrival of CO 2 , we see more of a staircase structure, with the mid-depth region of the aquifer delivering a CO 2 volume fraction of 12 % at 20 days after injection, followed by the other two regions at 128 and 148 days. This gives much better comparison with the field observations, indicating that a three-layered model is more appropriate than a single-layered model if one is interested in predicting the first arrival of CO 2 (e.g. the first tracers of CO 2 were detected at Salt Creek 10 days after injection) and the arrival distribution, but less useful if one is only interested in predicting the breakthrough of the bulk CO 2 quantity (65 days). More generally, there is an interesting question about how many upscaled layers are needed to accurately capture the CO 2 injection in a given aquifer. By breaking the aquifer up into smaller and smaller subdomains, we can achieve better and better comparison with field data, but at some point this defeats the point of upscaling, since our original objective was to avoid resolving all the heterogeneities.
The main implications from the comparison with Salt Creek are threefold: Firstly, we have shown that bulk CO 2 breakthrough times can be reasonably well predicted by our single-layered upscaling approach, although with an over-predicted volume fraction. Secondly, we illustrated that by breaking the aquifer into three subdomains, a much better comparison with field data is achieved, including realistic predictions of CO 2 volume fraction at the observation well. Finally, we have shown that there is clearly significant difference between capillary and viscous limit predictions, indicating that an accurate flow description requires careful modelling of the heterogeneities.
It should be noted that, instead of treating the heterogeneities by upscaling, one could alternatively use a three layer viscous limit model (such as figure 11d,e), which effectively assumes each layer is homogeneous with different mean permeability, and achieve good comparison with the data by fitting the other parameters of the problem. However, this would result in misleading parameter values that account for the heterogeneities indirectly. Therefore, a thorough upscaling approach is more advisable.

Concluding remarks
We have studied the effect of a vertical heterogeneity in a porous medium on the overall flow properties by way of upscaling. This is characterised by the two limiting cases of large capillary number (viscous limit), where heterogeneities play a weak role, small capillary number (capillary limit), where heterogeneities play a dominant role, and intermediate capillary number, for which a balance is sustained. In the former limiting cases we derived analytical expressions for the upscaled equivalent relative permeabilities using asymptotic analysis. For intermediate capillary numbers we used numerical simulations to suggest a composite (heuristic) form for the equivalent relative permeabilities that remains accurate across all flow regimes. These composite expressions enable fast, efficient simulations of dynamic flow configurations, which highlight where and when heterogeneities play an important role in the dynamics. The CT scan experiments of Jackson et al. (2018) were used for comparison with some of these upscaling results.
Using an analysis that stemmed from the classic Buckley-Leverett problem (Buckley & Leverett 1942), we applied the upscaled quantities to describe the flooding of an aquifer with heterogeneities. We illustrated how and when heterogeneities accelerate/decelerate the dynamic flow. By extending this analysis to the case of a radially symmetric injection, we illustrated how the capillary number at the flooding front changes over time. At early times, near the source, the front is in the viscous limit regime (where heterogeneities are unimportant), whereas later on, far away from the source, it is in the capillary limit regime (where heterogeneities dominate the flooding speed). The implications for CO 2 sequestration are that heterogeneities can alter advection of CO 2 by as much as 50 %, indicating the need for modelling such effects, as illustrated by our comparisons with field data from the injection experiments at Salt Creek, Wyoming. Finally, we illustrated how the choice of length scales for upscaling significantly affects predictions, underlining one of the key outstanding challenges in this field.
For future work, the effects of a dynamic flow on the equivalent properties could be investigated (i.e. instead of steady-state upscaling), using some canonical time-dependent case studies. This would be particularly useful for understanding when steady-state upscaling is an accurate approach, and when more detailed models are necessary.
Perhaps the clearest direction for future work would be to include the effects of gravity. One can quantify the importance of such effects by considering the ratio between buoyancy forces (due to gravity) and capillary forces (due to heterogeneities). Hence, gravitational effects are characterised by the so-called Bond number, which is defined as Bo = (ρ w − ρ n )gH/p e 0 . Inputting parameter values that are typical for CO 2 storage applications, we find that gravitational effects can be neglected (Bo < 1) for aquifers that are thinner than around H ∼ 1 m. For situations other than the thin aquifers we have studied here, gravity plays a significant role (Nordbotten & Celia 2011). This manifests in a saturation distribution that is determined not only by the heterogeneities (as we have shown here) but also by a stratification due to gravity (pushing the buoyant CO 2 upwards). This stratification would in turn affect the upscaled flow properties, such as the equivalent relative permeabilities, and consequently the upscaled flooding speeds.
For aquifers that are sufficiently wide, or which have sufficiently large Bond number, vertical variations in the flow become so significant that an averaged approach, such as the one taken here (in the Buckley-Leverett problem), is no longer applicable. Indeed, if the CO 2 is sufficiently buoyant that it rises and detaches from the lower aquifer boundary, a gravity current forms on the upper boundary and spreads laterally (Golding et al. 2011). The effects of heterogeneities on such flows have recently been investigated numerically by Jackson & Krevor (2020). As a future step, the simple upscaling work presented here could be extended to account for such cases. In other situations, the interface may remain confined above and below, and instead propagates through the aquifer as a buoyant plume. Recent studies have shown how vertical heterogeneities can alter such flows in the case of miscible fluids (Hinton & Woods 2018). It would be interesting to compare and contrast such results to the immiscible case, for which a simple upscaling approach as taken here would be fitting.
Another common challenge in hydrology applications is estimating rock heterogeneities, where it is often only possible to obtain very sparse measurements. It would be interesting to use our analysis here to explore the inverse problem of estimating rock heterogeneities from a small number of data points of the equivalent properties of the flow (and mean saturation). This would be easiest in the case of small capillary number, where one could use the function p e (z) and the power law B to fit the equivalent relative permeability curves to measurements. This approach is unlikely to be well posed, since multiple types of rock heterogeneity may give the same upscaled properties, but still one could develop an ensemble of likely heterogeneity profiles as an informative tool for geoscientists. Firstly, the model of Corey (1954) used by Golding et al. (2011) for Ellerslie sandstone is given by where the parameters are given by k rn 0 = 0.116, α = 2, β = 2, S wi = 0.651, λ = 1. The value of p e 0 is not given. Secondly, the model of Chierici (1984) used by Jackson et al. (2018) for Bentheimer sandstone is given by where the parameters are given by M = 0.65, L = 0.75, A = 3, B = 5, S wi = 0.081, λ = 2.3 and p e 0 = 3.51 kPa. Finally, the Brooks-Corey model (Dullien 2012) used by Krevor et al. (2012) for the Paaratte sandstone is given by where the parameters are given by k rn 0 = 0.95, α = 2, β = 8, S wi = 0.05, λ = 0.9 and p e 0 = 2.1 kPa.

Appendix B. Parameter values and extra plots
In this appendix we present the parameter values used for the Salt Creek case study in table 1 as well as two extra plots: figure 12 shows the equivalent relative permeabilities calculated for a step permeability profile, and figure 13 shows the viscous and capillary limit predictions for the time-dependant volume fraction of CO 2 for the Salt Creek case study.
For figure 12 we apply a similar methodology as outlined in § 2 except, instead of using a sinusoidal permeability profile, here we use a smoothed step function given (in using parameter values A = 0.5,ẑ 1 = 1/4,ẑ 2 = 3/4 and ζ = 40, which is illustrated in figure 12(c). All other parameters are taken as the same as in figure 2 for the sinusoidal permeability case. In figure 12(a,b) we display the equivalent relative permeabilities calculated numerically for different values of the capillary number and mean saturation, as well as best fit curves using the composite expression (2.61). For figure 13 we use the same approach as outlined in § 4.2 for modelling injection at the Salt Creek case study except, instead of using the composite expression (2.61) for the equivalent relative permeabilities, we use the viscous limit (figure 13a) and the capillary limit (figure 13b).