Exact Solutions for Gravity-Segregated Flows in Porous Media

: The review is devoted to exact analytical solutions for quasi-2D gravity segregated flows or gravity currents in subterranean porous formations. The problems under consideration are quasi-linear. The driving forces are two components of the buoyancy — one exerting the bulk of the light fluid and one due to the curvilinearity of the interface between the fluids. In the case of homogeneous formation or where the seal slope is negligible, the transport equation is parabolic and allows for a wide set of self-similar solutions. In a large-scale approximation of the buoyancy domination, the governing equation is hyperbolic; the method of characteristics allows for a detailed analytical description of gravity current propagation with final accumulation in the geological trap. Analytical models for leakage via the caprock seal are also discussed. The work was completed by formulating some unsolved problems in segregated flows in porous media.


Introduction
Segregated flows or gravity currents occur in numerous natural and industrial processes. The incomplete list includes secondary migration of gas/oil in geological basins, underground storage of CO2, hydrogen, and natural gas in aquifers or depleted oilfields, water-and gas-flooding of oilfields, vadose zone flows, and plant irrigation. The understanding and decision-making in those processes are based on mathematical modelling. Figure 1 shows the propagation of the heavy-fluid bank over the solid substrate with upward injection. The flow is three-dimensional (3D), reflecting the timely propagation of the bank in the substrate plane and the vertical profile along the slope (Figure 1). Depending on the interplay between viscous, capillary, and gravitational forces, several asymptotic flow regimes can be distinguished [1][2][3]. This includes viscous domination where the displacing phase moves in highly permeable layers, capillary-dominant filling of the layers in the order of increasing their permeability, gravity domination where the heavier/lighter fluid moves over the bottom/top of the reservoir, flow in commingled layers, and several mixed regimes. The above regimes allow averaging over the vertical, reducing 3D flow problems to 2D, and 2D to 1D.
In this paper, we consider 2D gravity-dominant flows. The lighter phases (gas, oil) flow above water, which are separated by the interface h = h (x,t) in unconfined and confined reservoirs (Figure 2). A heavy phase influx yields an unconfined flow at the early time and a confined flow later ( Figure 2); boundary conditions at x = 0 for the unconfined flow are transferred to the moving front x = xN2(t) for the confined flow. Figures 1 and 2 correspond to gravity-stabilised water injection, while Figure 3 shows gas migration below the seal. The gas migrates from the injection or expulsed zones towards the trap. Injection of the limited gas volume yields gas front propagation (drainage) at the advanced edge and thinning down the gas layer (imbibition) after completion of the gas pulse injection. To reflect the seal asperity, the cross-section of the seal stringer is assumed to have a power-law shape ( Figure 4). Fast gas propagation and a thicker gas stream are expected in the case of a sharp cross-section of the flow channel.  . Two-phase displacement under gravity domination: the current propagation is the unconfined reservoir at the early time; the current becomes confined at a late time (reproduced from [5]).
Exact solutions for porous-media flows allow for clear flux structures, fast multivariants, sensitivity-study calculations, and quasi-3D streamline modelling. Moreover, the exact solution for a direct problem may regularise the inverse ill-posed problem, providing a determination of the model functions from the laboratory coreflood or fieldscale tests, such as in two-phase or colloidal flow tests in porous media [6][7][8][9]. This explains the intensive current research on the exact integration of the equations of flow in porous media [3,[10][11][12][13].   The structure of the manuscript is as follows. Section 2 derives governing equations for segregated flows in porous media and describes types of the exact solutions. Section 3 presents the exact solutions for the general case of the parabolic transport equation. Section 4 reviews the exact solutions for the case of the large areal scale, where the governing equations are hyperbolic. Section 5 shows the solutions that account for leakage of the injected/expulsed fluid through the seal. Section 6 formulates some unsolved problems in the exact integration of the segregated flows in porous reservoirs. Section 7 summarises the critical review and the gaps in exact integration of gravity currents in subterranean porous formations.

Transport Equations for Water Displacement by Air/Gas/Oil in Sealed Channels
This section presents the basics for governing the equations of stratified gravity currents in porous reservoirs. It includes physical assumptions of the model (Section 2.1), derivation of the flow equation (Section 2.2), its dimensionalization (Section 2.3), the change of variables allowing simplification of the governing equation (Section 2.4), and some asymptotical forms of the governing equation (Section 2.5).

Assumptions of the Model
We discuss the two-phase migration of light fluid by buoyancy through the dipping sealed path. To be specific, we call this fluid "gas", while the model is fully applicable for air, CO2, H2, or oil too. Figure 3 shows the dipping top of the geological formation, primary expulsion of gas/oil below the seal, and propagation of gas inside the tortuous path due to buoyancy. This schematic is fully applicable for continuous or pulse injection of CO2 or other gases in the aquifer.
We assume incompressible gas and water. There is a steady-state flow of water in the carrier bed parallel to the seal. Gas expulsion and its buoyant propagation do not perturb permanent pressure distribution in the water flow zone below the gas-water contact (GWC). We neglect the flow in gas and water zones perpendicular to the seal, i.e., pressure is distributed in the flow cross-section hydrostatically. We account for gas dissolution in water but neglect water evaporation into the gas phase.

Derivation of Governing Equations
Following [14] and [10,11], here, we present the governing equation for segregated buoyant-driven flows. Axis x in Figure 3 is directed along the inclined path, and axis z is perpendicular to the seal; θ is the slope of the inclined layer. The thickness of the gas layer is h(x,t) and the cross-section of the path is A(h). Figure 4 shows different stream forms of the power-law cross-sections A(h) = h l for l < 1 and l > 1.
Darcy's law for water velocity is where W is the velocity of water, k is the permeability, p is the pressure, μ is the viscosity, ρ is the density, and g is the gravitational acceleration. From now on, the subscripts g and w correspond to gas and water, respectively. From the assumption of constant water velocity, the pressure gradient in the domain below GWC is constant, so the pressure is linear with x.
Pressure in the gas zone is determined from hydrostatic distribution in direction z: Here, ρ is the density difference between water and gas, and Pc is the capillary pressure that is equal to the pressure difference in gas and oil (non-wetting and wetting phases, respectively).
Expressions for the pressure gradient above GWC along with Darcy's law for gas allow expressing flow velocity in the gas zone via the WGC position h = h(x,t): sin cos where U is the gas velocity, and Krgcw is the relative permeability for gas at connate water. Introduce two reference velocities sin , for water and gas, respectively. The expression (3) for gas velocity becomes Gas velocity U is constant over the stream cross-section, so gas flux through the cross- Gas displaces water (drainage) during the increase in the stream thickness, and water displaces gas (imbibition) during the h(x,t) decrease. Figure 3 shows that drainage occurs at the advancing gas stream, where gas saturation below GWC is zero; gas saturation below GWC is equal to Sgr during imbibition after the pulse expulsion.
The continuity equation for gas is where residual gas during imbibition is accounted for by ,0 R is the equilibrium gas concentration in water, Q is the volumetric sink term that accounts for the dissolution of gas in water through gas-water menisci in GWC, and leakage through the seal.
Mass balance Equation (6) accounts for mobile gas in the stream, for gas dissolved in connate water, for residual gas entrapped in advancing water during imbibition, and for gas dissolution in water via menisci in GWC.
The expression for the sink term is where the saturation in seal s is determined by equality between capillary and overburden pressures, ωl is the volumetric leakage rate,  is the kinetic dissolution constant, B(h) is the bed width, and the leakage constant ωl corresponds to permanent gas flow through the breached seal. For no overburden h = 0, there is no leakage, i.e., the seal conductivity ωl = 0. The higher the overburden, the larger the conductivity, i.e., ωl(h) is a monotonic function. The seal conductivity becomes constant where the overburden exceeds maximum capillary pressure, and the saturation is stabilised at the level 1−sgr. The above justifies the assumption that the seal conductivity remains constant while the current thickness varies.
Dissolution through molecular diffusion is a relatively slow process. However, CO2 dissolution in formation fluid induces density changes in the fluid that can trigger the density-driven convective flow and, thus, enhance the dissolution rate. Indeed, carbonised water is heavier than formation water; the transition processes of gravity segregation with interface brine-gas mass transfer can result in significant gas flow through the GWC surface.
Accounting for gas incompressibility allows cancelling gas density; substitution of the migration flux expression (8) into mass balance (6) yields Equation (9) is an advective-diffusive equation with hysteresis and a sink-source term [15].

Dimensionless Variables and Governing Equation
Introduce the following dimensionless variables and parameters: where hL is a typical current thickness. Equation (10) presents the transformation of dimensional parameters to dimensionless. So, all parameters before Equation (10) are dimensional (and afterwards-dimensionless). Substituting dimensionless variables and parameters (10) into Equation (9) yields the dimensionless equation Now, we introduce a new unknown A = A(x,t), where C(A) is a non-linear function determined by parametric dependencies A(h) and C(h).
Expressing C(h) as parametric dependency C(A) makes A(x,t) unknown in Equation (14) ( ) The second term on the left-hand side of Equation (15) corresponds to buoyant flow below the dipping seal and disappears in the case of plume spreading below a horizontal caprock. In this case, and for C = A 2 , Equation (15) is attributed to Boussinesq and allows for numerous exact solutions ( [15][16][17]-see the second line of the third column in Table  1). Table 1. Similarity solutions for gravity-segregated flows in porous media.

General Form Early Time Late Time References
Governing equation Unconfined media

Transformations, Change of Variables
Now, we transfer to the reference system that moves with the unitary speed Substituting (16) into Equation (13) yields Consider the power law for the dimensionless stream cross-section From Equation (18), it follows that Equation (17) becomes Changing the time scale from t to  and introducing power k , 11 l k ll   == ++ (21) transforms Equation (20) to the form Self-similar solutions. Now, we look for the solution to Equation (22) Constant n in Equation (24) is determined from either boundary condition [15] 0: x h qt  == (25) or gas mass in the system [19] ( , ) The initial condition for the non-perturbed system is 0 :

Extensions to the Model
The 2D solution for the two-phase flow tends to the fully segregated system, given by the equations in the previous section, where the reservoir thickness tends to infinity. If the reservoir thickness tends to zero, the system tends to 1D Rapoport-Leas equations [1,2,14].
Introduction of the layer-cake heterogeneity k = k(z) yields the change from A to its non-linear functions in all terms of mass balance Equations (9) and (11), i.e., accumulation, advective, and dispersive fluxes become non-linear functions of gas stream thickness in thick reservoirs (saturation in thin reservoirs).
It takes some delay time for those functional dependencies to be established, so mass balance Equations (9) and (11) contain the current values, which are linked to the stabilised values by kinetic relaxation equations [15]. In the case of stratified flows, the relaxation delay is the time to establish accumulation and flux terms: Here, u serves as the plume (migrating layer) thickness h in thick reservoirs and as saturation s in thin reservoirs; x is the dimensionless linear coordinate, t is the dimensionless time, Gc is the mass of accumulated matter, Fcv is the advective flux, and Φ is the potential for diffusion flux. Subscript cv corresponds to the current values of the accumulated density, flux, and potential that reach the equilibrium values after the delay times τk, where k = G, F, and Φ, respectively. The dependencies G(u), F(u), and Φ(u) are valid under the stabilised flow conditions. The dimensionless groups , , ,, are expressed via the flow velocity U, the reference length L, and the delay times τk. The initial conditions are given for all current values Inlet boundary conditions are given for either total flux The flow-through outlet boundary condition is The impermeable boundary generates the no-flux boundary condition The outlet boundary conditions for the flow in the semi-infinite reservoir correspond to not reaching infinity by finite times: The same conditions at minus infinity are set for the flow in the infinite reservoir. The asymptotic cases governing system (28) allows for three simplified asymptotic regimes: large-scale, low-velocity, and high-velocity approximations, where initialboundary problems allow for either exact or asymptotic solutions [24].
The large-scale approximation corresponds to such a large L that the groups tend to zero (are negligible). The dissipative term in Equation (36) disappears, unknowns G, F, and Φ take their equilibrium values, and system (28) degenerates into a scalar hyperbolic equation Boundary conditions (31) degenerate into Figure 5 presents a graphical solution for the Riemann problem with constant initial and boundary conditions. Different forms of fractional flow function F(u) are shown in Figure 5a. Upper/lower dashed curves correspond to reservoir scale displacement of higher/lower viscosity fluid than the displaced one, respectively. The middle continuous curve corresponds to core-scale displacement. Upper/lower continuous curves reflect topdown/bottom-up flows. Point J corresponds to the boundary condition (injection). Point I corresponds to the initial condition. Non-self-similar solutions for stream propagation until the trap and its accumulation are obtained by methods of characteristics and interactions of non-linear hyperbolic waves [24].
The low-velocity approximation corresponds to such a small U that the kinetics groups are negligible, i.e., the unknowns G, F, and Φ take their equilibrium values System (28) degenerates into the non-linear parabolic equation The high-velocity approximation corresponds to such a large U that group  is negligible, i.e., the diffusion term disappears System (28) Three kinetics equations for the current values separate from the first Equation (44) and are solved independently.

Parabolic Nonlinear Partial Differential Equations
The gravity-driven propagation of viscous fluids has been the subject of extensive research since the 1950s [11,17,[25][26][27]. Usually, the propagation of a denser fluid underneath a lighter one is considered. The spread of the heavier fluid is obtained by a balance of viscous and buoyancy forces at a large Bond and a small Reynolds number, which results in a system of nonlinear parabolic PDEs. The solution of this system predicts the spread of the gravity current shape h as a function of x and t.
The parabolic PDE of buoyant currents "forget" the initial conditions from which they evolve, for some intermediate period t; this behaviour allows for an intermediate asymptotic [15,25], which results in a self-similar current shape. The self-similarity conversion reduces the PDE to an ODE. If the similarity forms are found by a scaling analysis, the solution is known as a first-kind self-similar solution [25].
In the case that the similarity variable cannot be found by a dimensional (scaling) analysis, the problem denotes a second-kind self-similarity [25]. The second-kind selfsimilar solution was classified by Gratton and Minotti 1990 [41]. They reduced the governed nonlinear PDEs to autonomous nonlinear ODEs.
This section contains a thorough review of the work on the first-kind non-linear parabolic PDEs for Newtonian buoyant currents. Huppert (1982) [19] derived governing equations that describe the motion of gravity current generated by a line source in wide horizontal channels. The negligible shear on the gravity current's upper surface was assumed. Using the balance between the pressure drop and the viscous forces along with vertical integration of local mass conservation, he derived the following nonlinear PDE to predict the propagation of the gravity current h(x,t) where ϑ is the dynamic viscosity of the injected fluid. He considered that the volume of the injected fluid increases with time as t  , so the global continuity equation becomes Here, xN is the current front, i.e., h(xN) = 0 (see Figure 2),  = w − c, and q is the injected flux volume, which allows a self-similar solution. Equation (46) also covers the injection of fixed volume (α = 0) and a fixed flux (α = 1). The similarity solution was expressed in terms of the following similarity variables (47) and the solution has the form of where is the value of at x = xN and function is found from For the general form of , Equations (49) and (50) were solved numerically; however, for the release of constant volume, i.e.,  = 0, they can be solved analytically Lister 1992 extended the previous work by [19] and derived the system of equations to describe the flow of the gravity current from a point or line source in an inclined wide channel. He assumed that the volume released from the point or line sources is proportional to t  , resulting in the following systems of dimensionless equations for 2D and 1D flows, respectively. By looking at possible asymptotic balances, and changing the dimensionless coordinates, he found that, at the early times, Equations (53) and (54) could be approximated by In this case, the gradient of the interfacial thickness is much bigger than that of the plane, as on a horizontal plane, the fluid spreads symmetrically from the source, as described previously by Huppert 1982 [19].
For late times, Lister 1991 approximated Equations (53) and (54) as The similarity solutions for the long-time behaviour of the system of PDEs were obtained and the results showed that the current was mainly a downslope, with some cross-slope motions for the case of release from a point source.
The approximated equations for the case of the injection from the line source, Equation (58), are hyperbolic. Section 4 reviews the analytical solutions of the hyperbolic system of the equations. Table 1 summarises the analytical solution for early-and latetime asymptotic models.
As mentioned earlier, numerous studies have used a sharp interface approximation to describe the movements of buoyant currents in porous media, in which the areas inside and outside of the current are assumed to be fully saturated and these fluids are disconnected by a sharp interface. The current is confined if it fills the entire depth of the reservoir, usually close to the release point; in this case, the movement of the existing fluid has a significant dynamic role. However, when the height of the current (h) is less than that of the porous media (h/H << 1), the current is unconfined, and the dynamic role of the existing fluid is negligible [10,22,42,43]. Figure 2 shows the constant rate injection of a fluid into a porous media containing an ambient fluid with lower density; the current is unconfined at the early time of the injection and is confined at a late time [5].
For the confined current in a horizontal reservoir, the flow is driven by changes in the local hydrostatic pressure and by the background pressure gradient; the velocity of the current can be obtained using Darcy's law [17,44], along with vertical integration of the local mass conservation that yields ̂ is the unit normal in the direction of the linear coordinate x, ub = kg/c, and M = c/w.
The problem can be solved subject to the global mass conservation equation where D is the 2D region filled by the injected fluid.
For the unconfined current, where h/H << 1, the second term in Equation (59) can be neglected and the third term becomes a function of h only-independent of H and M-as was presented by Hesse et al. 2007 [22] for the spread of a fixed volume release in a 2D problem.
If the currents move along an inclined layer at angle θ, the upslope component of gravity must be considered, so that Equation (59) for the unconfined inclined current becomes Here x is the unit vector in the along-slope direction. Vella and Huppert 2006 [4] showed that for early time (t < t * = (q 2 /ub 3 tg) 1/(3−) ), the spreading current is axisymmetric. At the late time, the current moves downslope, and the upslope extent grows as t.
Many theoretical studies have focused on unconfined currents [4,10,15,21,[45][46][47][48][49][50][51][52][53][54]. The self-similar solution in which the scaling exponent can be obtained by balances of terms in the system of equations was presented by [15,55] to define the propagation of a fixed mass of fluid. Huppert and Woods 1995 [10] extended the previous work to predict the spread due to continuous fluid injection. They derived a series of solutions to predict several gravity-driven currents in unconfined porous media, accounting for permeability changes in the direction normal to the boundary. They assumed the permeability was related to the porosity by = 0 and porosity changes linearly in the direction normal to the boundary = 0 + 1 . The governing equations describing the motion of the gravity current for the small values ϕ1h/ϕ0 become Here, δ = 1/2ϕ1/ϕ0 and S = ϕ0 n−1 ρgk0/µc. For the case of injection of fixed volume and fixed flux of denser fluid in a horizontal,  = 0, 1D porous layer, Equation (62) reduces to For the case of α = 1/2, the system admits the similarity solution, where the current length increases as t 1/2 . For the limit →0 and the injected flux equal to qt α , Equation (64) satisfies a similarity solution where the current length increases as t (α+1)/3 . Further, Huppert and Woods 1995 [10] studied buoyant current propagation along a sloping layer, i.e.,  > 0. In the limit of constant permeability, Equation (64) reduces to For the release of constant volume α = 0, the current moves steadily in the downslope direction, and moves diffusively due to the gravity normal to the layer, as on a horizontal boundary.
For the release of constant volume α = 0, considering the permeability variation Equation (65) It was shown that at both early and late times, the propagation of the fixed injection of fluid was defined by self-similar solutions, where t * = q 2 /(2Sn 3  3 tgθ 2 sin θ) is the critical time at which both solutions calculate the same shape of the current. At early time, t << t *, the main motion relative to the moving frame consists of the gravitational slumping of the current caused by the gravitational component normal to the boundary and, therefore, they recovered the solution described for the case of constant permeability Equation (71). For the late time, t >> t *, the asymptotic behaviour was described as Vella and Huppert 2006 [4] studied the release of relatively heavy fluid from a point source into an unconfined porous media (Figure 1). The governing equation is as Equation (61) and the system was closed by posing that the volume of the current increases with time, such as qt α , as in Equation (60). For the case of fixed volume release α = 0, introducing the moving coordinate z = x − t transforms Equation (61) to the form that admits an axisymmetric similarity solution For the case of constant volume release, at a very long time, Equation (61) was reduced to Which predicts the steady-state shape of the current. Equation (71) has a similarity solution in the form where function satisfies Which has a solution of For the case where the injected volume increases with time as t α , it was shown that for α < 3 at early times, the current moves asymmetrically, with radius r ∼ t (α+1)/4 , while at long times, it moves mainly downslope. For α > 3, this situation is reversed (see Table 1). Heterogeneity always exists in porous reservoirs, and it can considerably affect the spread of the buoyant current; as mentioned earlier, Huppert and Wood 1995 [10], obtained a system of equations for the buoyant current accounting for a power law permeability change in a vertical direction. A similar assumption was also seen in work by [32,56,57]. However, Zheng et al. 2014 [23] argued that heterogeneity exists in a horizontal direction as well as the vertical one. They studied the effects of horizontal heterogeneities on the motion of buoyant currents, by assuming a heterogenous reservoir with the power-law variation in porosity and permeability where Ap = k0ρg/(µcϕ0), the global mass conservation takes the form by introducing a similarity variable The shape of the current becomes The approximated form of f near to = 1 was obtained While previous work focused on the cases of gravity currents with constant density, Peglar et al. 2016 studied the propagation of gravity currents with variable densities in porous media. New equations for unconfined currents were derived, similarity solutions accounting for the fixed volume release and the constant injection rate of fluid with a variable density distribution were obtained. The results show that the density variation significantly affects the vertical motion of the gravity current. An analytical solution for the constant volume release where the density contours stratify horizontally was obtained.
So far, we reviewed the work for gravity currents in unconfined geometries, yet there are critical cases where the porous reservoir is confined by a second boundary (h ∼ H); sharp interface models have been obtained for buoyant currents in the confined media, where the flows of both fluids are significant and, therefore, their viscosities play an important role [22,43,[58][59][60]. In contrast, in the unconfined reservoirs, the flow of the existing phase is insignificant and, therefore, models of their spread are not functions of the fluid viscosities [4,10,21]. Huppert and Woods 1995 showed that the limitation in the total height of the flow in a confined reservoir results in an extra pressure drop related to driving the fluid along the length of the reservoir. Early studies of confined currents studied the exchange flows between layers occupied with phases of different densities but the same viscosities [10,17]. By using thin layer assumptions, these studies defined a similarity solution, revealing a linear interface between the fluid layers. Nordbotten [18] presented general 3D governing equations for the gravity current flow in vertically confined porous media, showing that the motion of the current is dependent on two distinguish contributions: evolution due to changes in hydrostatic pressure and pressure gradients provided by the injection. They considered an injection of fluid with viscosity μ and density ρ into vertically confined porous media that contained an ambient fluid with viscosity μa and density ρa and ρa > ρ. Assuming a purely hydrostatic force balance, they determined the pressure in the upper and lower layers as where  = ρa − ρ and P(x,0,t) is the unknown. By substituting (81) into Darcy's law, they determined the horizontal velocities. Vertical integration of the continuity equation across the porous media and substituting them in the Darcy velocities yields which forms a system of equations for two unknowns, h and P. To model the interface propagation, they imposed 0, where XL and XU are contact lines along the lower and upper boundaries.
Equation (82) (84) where q is the volumetric flux per unit width along the length of the medium.
Obtaining ∂p/∂x from (84) Therefore, Equation (82) becomes which is similar to Equation (59). As Equation (87) shows the velocities u and ua are separated into two contributions: u p which is driven by the background pressure drop; u g is driven by the local hydrostatic buoyancy gradient. Depending on the relative magnitude of u p and u g , Equation (87)  The former provides a gradient of the interface at the source in unconfined porous media.
The governing Equation (87) in the dimensionless form becomes Here, The boundary condition for the unconfined flow becomes Conditions (94) are imposed only after the current contacts the lower boundary. It was shown that, at the early time, the current was much thinner than the height of the reservoir and the system of equations simplifies to those systems in unconfined porous media.
For this case, the solution is similar to the previous work by Huppert and Woods 1995 [10]. They assumed the interface between two fluids was between two contact lines along the lower and upper of the confined medium. For the late time, they approximated Equation (59) to the nonlinear hyperbolic equations. A similar hyperbolic system of the equation was presented in [5]. They also studied the flow gravity current behaviour in a confined porous medium, when a fluid with a higher density was injected into a porous medium saturated with a lower density fluid. For the injection of fluid from a point source, they derived a nonlinear convection-diffusion, as was presented by Huppert 2006 [11] (Equation (59)). For the early time, Equation (59) was reduced to Equation (95).
For the constant injection of the flux, the global mass constraint is written as The self-similar solution was constructed with a similarity variable = / 2/3 , and the front spread as Zheng et al. 2015 [5] provided the asymptotic behaviour of Equation (97) near the current front, s→1 They showed that in the late time, Equation (59) could be approximated by a nonlinear hyperbolic equation.
Pegler et al. 2014 [18] solved Equations (91)-(94) numerically; they used partially implicit finite-difference schemes of the second-order initialised using the early-time asymptotic solution described in [10]. The numerical results showed that, at the early time, the current was in contact only with the upper boundary and had a frontal contact line XU that grew in proportion to t 2/3 . At a certain time, the current contacted the lower surface, forming the second contact line XL which spread sharply from x = 0 with a finite velocity. While the interface h for M = 1 to 10 remained nearly linear to the later time, M = 0.1 instead formed a concave shape. In all cases of M, the two contact lines XU and XL eventually transitioned toward linear growth in time. They compared the early time and late time asymptotic analytical solutions with the numerical results; it was shown that the analytical solutions agreed well with the early time numerical solution at t = 0.123 for M = 1 and 10. There was a weaker agreement when M = 0.1. This was in accordance with the asymptotic condition t << M 3 becoming sufficiently restrictive-that is, not satisfied when t = 0.125. While the asymptote of Equation (95) is independent of M, the time scale on which the flow transitions away from is strongly dependent on M.
As mentioned earlier, numerous studies of buoyant currents in porous media assumed that the two fluids are separated by a sharp interface and the area with and without the current are fully saturated. However, these assumptions are only valid for two-phase buoyant currents when capillary forces are insignificant and related to viscous and gravitational forces. Golding et al. 2011 claimed that in the case of the two-phase gravity currents, capillary forces can cause a saturation transition region that is not small; therefore, the sharp interface assumption is not correct. They considered the spread of fluid in a semi-infinite horizontal reservoir; the local saturation was defined by the vertical balance between gravitational and capillary forces. They modelled a vertically integrated equation for the buoyant current shape that accounted for non-uniform saturation variation. This model is equivalent to the one-phase models with the two-phase effects described by new functions defining the saturation at the impermeable layer and the vertically integrated height-dependent injected flux [10].
Hesse et al. 2007 [22] investigated the propagation of a fixed injection of fluid into a confined, 2D horizontal reservoir saturated with fluids of different viscosities and densities. They showed that, at the early time, the injected fluid covers the entire height of the reservoir, and the governing equation allows for a self-similar solution that is dependent on the fluid viscosities. The solution shows the current moves as t 1/2 . At the late time, the current moves along the boundary, and its height is much smaller than the thickness of the layer. The system of equations simplifies and permits a different similarity solution that is not a function of the fluid's viscosity. The late solution showed that gravity current speared as t 1/3 . Golding and Huppert 2009 [49] studied the effect of confining boundaries on the buoyant current in a horizontal reservoir. The media is assumed to have a uniform crosssection with height d satisfying the relationship d = ab(y/a) l where y is the cross-section coordinate and ab is the length scale of the channel width (see Figure 4). Therefore, the current cross-sectional area is given by r c is a dimensionless number of proportionalities between the channel width and height. Local mass conservation in a porous medium yields They derived a similarity solution for the case when the injected volume increases as t α For the horizontal medium, Equation (102) reduces to That allows for the similarity solution in the form where f is found from (1) = 0 (108) Here, F1 = (l + 1)(1 − 2α)/(3l + 2) and F2 = (l + 1)(αl + l + 1/l(3l + 2)).
For the case of α = 0, the system allows for analytical solution Equation (106) shows the current spreads at t c , where c is a constant number to be determined. Golding and Huppert 2009 [49] analysed the dependency of constant c on the geometry of the confined boundaries, l, the injected volume of fluid, α. They found that independent of the boundaries' shape, there is a critical value for α = 1/2, above that increase in l causes a rise in c, and below that increase in l causes a decline in c. For the case of inclined porous media by a large angle, tg >> ℎ Equation (102)  and Equation (103) becomes Equations (111) and (112) are similar to the equation in the horizontal layer (with S change to Scos θ as presented by Huppert and Woods 1995). Therefore, the motion of gravity along the slope is independent of the channel shape.
Hesse and Khodabakhsh 2006 [65] obtained a 1D model that accounted for capillary trapping and the reservoir slope. The derived system of equations was solved numerically.
Recently Huppert and Pegler 2022 [66] presented a series of various transition times between the early time radial regime and the gravity current. They studied the spread of a relatively dense fluid from the injection point at the bottom of a horizontal unconfined porous medium. It was observed that in the 3D model, there was an initially growing hemispherical shift of a steadily growing radius r N (t) ∝ t 1/3 into a buoyant current moving over the impermeable base of the reservoir. After evolution of intrinsic time scales, the hemispherical regime transitions towards the axisymmetric thin-layer regime with radius r N (t) ∝ t 1/2 . In previous work, such as [5,18,47,58,67], the authors studied the final thinlayer current, in which stresses became hydrostatic to the leading order under the Dupuit assumption ( [17]).
Further to capillary trapping, dissolution of the injected fluid into the existing fluid has been a topic of several studies. Modelling of the dissolution phenomena needs incorporation of both diffusions of dissolved fluid far from the sharp interface, and evolution of convective mixing within the ambient fluid. A complete model for dissolution, which accounts for a background flow and a fixed horizontal extent of the current, is given by [68][69][70][71][72][73][74][75].
Other notable studies include Zheng et al. 2021 [76], who studied the propagation of a buoyant current in a variable width Hele-Shaw cell as an analogy for flow in fractures and flat heterogeneous porous media. Using a phase-plane analysis, they obtained a second-kind self-similar solution to define the propagation of the buoyant current during both the spreading and leveling periods.

Hyperbolic Nonlinear Partial Differential Equations
As we showed in the previous section, the mathematical modelling of the spread of gravity currents results in a system of nonlinear parabolic PDEs. In some asymptotic cases, the governed parabolic PDE can be approximated with the hyperbolic system of equations. This section reviews the hyperbolic form of equations for the gravity currents.
Lister 1991 [20] obtained the late-time behaviour of Equations (53) and (54), the approximated form for the case of an injection of a fluid from the line source is hyperbolic (see Equation (58)). He derived the exact similarity solution of Equation (58), and the obtained results showed that the current length spreads as t (2α+1)/3 (Table 1).
Zheng et al. 2015 [5] also showed that, at the late time, Equation (59) can be approximated to a hyperbolic form of If the injected phase is less viscous than the existing phase, the solution of Equation contains shock waves, Figure 2 shows the positions of xN1 and xN2.
In work by Pegler et al. 2014 [18], it was shown that for the late time, Equation (59) is reduced to the nonlinear hyperbolic equations similar to Equation (113) presented by Zheng et al. 2015. The similarity solution using similarity coordinate x/t was obtained Three late-time asymptotic solutions were presented, subjected to whether the viscosity ratio of the injected fluid to the ambient one was greater, less, or equal to one. For the ratio of less than one, the current motion due to the hydrostatic pressure gradient eventually becomes small, with the system containing a similarity solution that is not a function of permeability. For the ratio equal to unity, both contact lines approach the same order asymptotic position t. Yet, the lateral expansion of the interface grows as t 1/2 due to buoyancy spreading. For a ratio greater than one, a linear evolution of the contact lines develops; however, in a coordinate moving with the contact lines, the interface emerges toward a steady-state. The approximated equations and self-similar solutions for different cases of M are presented in Table 1 (rows 1-3). Golding and Huppert 2009 [49] showed that for the case of inclined porous media by a large angle , tg >> ℎ Equation (102) becomes hyperbolic They showed that independent of flux parameter α or channel shape parameter n, the current length is proportional to St sin θ. Hesse et al. 2008 [61] and Juanes and MacMinn 2008 [59] in separate work, studied the effect of capillary trapping on the shape of the gravity currents. They obtained exact solutions to the hyperbolic limit of sharp interface models. Hesse et al. 2008 studied a model for an inclined reservoir. They reported a similarity solution for the movement of the current with a mobility ratio equal to unity. Juanes and MacMinn 2008 [59] studied a horizontal reservoir with ambient groundwater flow and used the exact solutions to upscale capillary trapping to the reservoir scale. Junes et al. 2010 [43] extended their earlier work ( [59]), presenting a sharp interface model that accounted for gravity override, capillary trapping, ambient fluid flow, and the shape of the current during the injection. They presented a nonlinear advection-diffusion equation during the injection period as For the post-injection period, the equation has the form of For the case of the small Ng, Equation (121) simplifies to a hyperbolic equation They obtained the complete analytical solution, using the method of characteristics. During injection, they assumed that the injection was symmetric with respect to the well, and they solved the equation with respect to the following initial condition The solution was described as a function of the characteristic speed = / , as the G function has a concave shape, the solution of the Riemann problem contains a rarefaction wave During the retreat stage-after the injection stop-the current moves to the right subject to the ambient water flow. This solution continues a rarefaction wave. The solution for the left side is a convergent fan with a characteristic speed quicker than the drainage. The characteristics move onto each other at = 0.
Chase stage: After the front moves through = 0-because of the concave shape of function G-the imbibition front moves as a shock wave. A shock that spreads at a velocity described by The Hugoniot-Rankine condition Sweep stage: Once the current disconnects from the base of the reservoir, the solution is a continuous interaction of a rarefaction wave with a faster shock. The equation predicting the spread of this stag is where M' is the mobility ratio A similar hyperbolic system of the equations is presented in work by Szulczewski et al. 2012 [76]. They used the method of characteristics to solve Equation (124). The associated pressure equation was solved numerically.

System of Equations Accounting for Leakage of the Injected Fluid through the Seal
Acton et al. 2001 [77] and Pritchard et al. 2001 [78] considered leakage of a gravity current through the impermeable reservoir boundary. In their mathematical models, they assumed that the pressure gradient between the current and the ambient drives the fluid flow. Therefore, the leakage velocity has the following form where ωl is the leakage velocity, kB is the permeability of the bonding layer, and B is the thickness of this layer. Accounting for the leakage, Equation (61) This idea was improved by Neufeld and Huppert (2008) [52]; they assumed that the fluid remained continuously connected inside the porous reservoir; therefore, the leakage velocity was expressed as Here, l(r,t) is the area where the current drained past the low permeability layer. They expressed the steady-state radius of the current as In this model, the capillary pressure was neglected. However, capillary forces are important in the area between the media of differing geometries. Therefore, for leakage to occur, the hydrostatic pressure must overcome a capillary entry pressure pc. This results in a critical capillary height derived by Wood and Farces 2009 [79]: below which there is no leakage; above that, the fluid from the current can leak through the impermeable layer. Leakage through a range of different reservoir geometries has been studied in confined and unconfined reservoirs [4,[52][53][54]58,80].

Some Unsolved Problems
The main applications of segregated flows lay in gas storage, hydrology, plant irrigation, and geological formation of petroleum reservoirs. The propagation tips in fully segregated flows can be significantly thinner than the capillary transition zone; however, the analytical models ignore capillary pressure between gas and water.
Chemical reactions, such as oil biodegradation or CO2 and hydrogen reactions with the mineral rock components as well as rock dissolution are not accounted for. In this case, mass balance must be complemented by the chemical rate equation, such as the active mass law with reaction kinetics [81,82]. The analogous 1D two-phase flow problems allow for exact solutions [83][84][85][86][87][88][89]. The segregated models can be derived under the assumption of instant reactions.
Natural top-down water fluxes are not accounted for too. The analytical models also do not include pressure dependencies for viscosity and density. The analytical models with varying seal inclination angles and varying permeabilities along the flow trajectories are unavailable too. Another unexplored aspect is attic volume due to seal asperity and tortuosity, which can highly affect plume propagation.
Natural top-down water fluxes are not accounted for too. The analytical models do not include pressure dependencies for viscosity and density.
The segregated gas migrates below the geologically formed inclined asperous seals; water flows over the asperous tops of impermeable substrates. It makes the problem of flow in curvilinear streams with power-law shapes (A(h) = βh l ) realistic (see [11,19]); moreover, the l-value determines the form of self-similarity. The flows in sharp-shape channels are typical for secondary migration (b < 1), while for CO2 injection, thick layers (b > 1) are more typical flat streams. Despite the analytical solutions for any arbitrary l obtained [21], typical l-values are not determined from the basin data on petroleum accumulations. However, there was no attempt to determine l from the basin or field data. This topic is important for modelling-based predictions due to the strong l-dependency of the buoyant profiles. Besides, it would be useful to implement an exact solution for the flow in a single stream with a percolation-effective-medium network-flow model.
The dynamics of natural gravity streams can be determined by network modelling. The exact solutions can be determined using the percolation (for stream density) and effective medium (for flux) theories [90][91][92][93].
Gravity-driven segregated flows below the seals in infinite-thickness geological formations are typical for gas storage and secondary oil-gas migration, while oil displacement by water mainly occurs in thin reservoirs with limited thickness. The first case should be described as an asymptotic limit of the second case; besides thickness, the dimensionless groups should include gravity, permeability, density difference, etc. The above would determine the area of validity of segregated models and lead to richer analytical models for gravity stream propagation. Definition of the limits of dimensionless groups where the segregated models are valid would allow for the development of hybrid numerical models [94].
The composition of expulsed oil and gas vary from pulse to pulse; it can even vary inside one pulse. Adsorption on the rock can also separate different components and create a chromatographic flux. Numerous analytical models for a two-phase multicomponent 1D transport in porous media have been derived [13,24,95]. However, analytical models of segregated flows with multicomponent effects are currently not available.
One way around would be the 1D Riemann problem solution for the two-phase ncomponent flow along a seal, approximation of all rarefaction ways by jumps, so the Riemann n × n solution becomes a sequence of jumps. Then, the multifront stratified 2D model is assumed. The unknowns are positions of n fronts; phase saturation and component concentrations are assumed to be the same as in the 1D solution; n equations are conservation laws of masses for each component. For a three-component system of water-polymer-oil, and under the assumption of incompressibility, which reduces the number of conservation laws by one, two positions h(x,t) of water and polymer fronts are determined by a 2 × 2 hyperbolic system of conservation laws for water and polymer [24].
The important case of multicomponent flow is the so-called deep bed filtration of colloidal-suspension nanoparticles [96][97][98] widely occurring during gas injection and, consequently, in the segregated gas-water flows during gas storage. One way around derivation of the segregated flow model would be averaging using exact 1D solutions [85,[99][100][101].
Equation (9) contains a hysteretic accumulation term due to gas that resides in advancing water during imbibition. For the horizontal seal, the problem of gas plume propagation is mathematically equivalent to the so-called Barenblatt's second-type selfsimilar solution [15,102]. In a large-scale approximation, besides the hysteretic accumulation term, the hyperbolic flow Equation (37) contains end-point relative permeability that is also different for drainage and imbibition. The technique of Riemann problem and wave interactions developed for hysteretic hyperbolic problems [100,103] can be applied for segregated flow problems.
The segregated flow models solve the problems of sweep and storage volumes. The injectivity problems due to chemical reactions and fine migration can be solved by matching the analytical models for injectivity [101,104,105] with the segregated flow models.

Summary
The dynamics of segregated gravity streams in subterranean porous formations are described by a single 1D advective-diffusive equation with the sink term. This equation is obtained by averaging quasi-2D flows under the assumption of vertical equilibrium. The advective flux corresponds to a buoyant flow of a lighter-phase stream below a seal or in the network of conductive paths. The diffusive flux exterminates curvilinearity (straightening) of the interface boundary. The sink term corresponds to leakage through the breached seal and chemical reactions.
For any power-law of timely mass change in the migrating plume, the 1D problem allows for a self-similar solution. Power zero corresponds to the pulse injection. Power one corresponds to the constant rate injection. A self-similar solution can be formulated in terms of the inlet boundary condition for injection (Barenblatt 1996) or mass variation in the system (Huppert). Initial and boundary conditions for self-similar solutions have a power form, including the uniform one with the power equal to zero.
In a large-scale approximation, the diffusive term is ignored, and the governing equation becomes hyperbolic. Here, the shocks fulfill the conditions of mass balance, stability, and admissibility. Continuous solutions are obtained by the method of characteristics. The above yields an analytical solution for any initial-boundary problem (see Table 1).
Travelling wave solutions of a general advective-diffusive equation near the shocks in the reduced hyperbolic equation provide asymptotic matching of inner and outer expansions, reflecting the small-scale diffusive effects in a large-scale approximation model. The above technique is applicable for segregated flows with hysteresis, where the model coefficients are different for drainage and imbibition.
Exact solutions and related analytical models are extremely efficient in in situ gas storage, including CO2 geo sequestration and hydrogen storage, environmental engineering, and hydrology in partly saturated reservoirs and geological formations, agriculture, and secondary migration in geological basins. The exact solutions allow for predictive estimates of the system's behaviour, multi-variant sensitivity studies, optimising the numerical schemas, and a deeper understanding of the natural and technological processes.

Conflicts of Interest:
The authors declare no conflict of interest.

A
Current cross-sectional area, L 2 a Self-similar form of current cross-sectional area ab The length scale of the channel width, L B Seal thickness, L b Channel shape power parameter Initial values