1. Introduction
Underground gas storage plays a crucial role in addressing anthropogenic climate change, both by sequestering greenhouse gas emissions prior to atmospheric release and by enabling the transition to cleaner energy sources. The first approach is well established through carbon capture and storage, where large quantities of CO
$_2$
are injected underground for long-term containment. More recently, underground hydrogen storage has gained significant attention for its potential to enable large-scale integration of renewable energy into the power grid (Lord, Kobos & Borns Reference Lord, Kobos and Borns2014; Tarkowski Reference Tarkowski2019; Heinemann et al. Reference Heinemann2021; Zivar, Kumar & Foroozesh Reference Zivar, Kumar and Foroozesh2021; Muhammed et al. Reference Muhammed, Haq, Al Shehri, Al-Ahmed, Rahman and Zaman2022). The concept is to use surplus renewable energy generated during the summer months to produce hydrogen via electrolysis, store it in underground formations and subsequently recover it during the winter when demand exceeds supply. Hydrogen is a particularly attractive candidate for this large-scale energy storage owing to its high energy density and the fact that its combustion or use in a fuel cell produces minimal harmful emissions (Zivar et al. Reference Zivar, Kumar and Foroozesh2021). Proposed underground storage sites for gas include depleted reservoirs, aquifers and salt caverns, all of which are typically saturated with brine. Unlike CO
$_2$
storage, where the gas is intended for permanent sequestration, hydrogen storage requires recoverability, making accurate predictions of gas plume evolution crucial. Due to its exceptionally low viscosity relative to brine, injected hydrogen may spread rapidly along the top of the brine within an aquifer, forming a thin gas layer that extends far from the injection site (Hagemann et al. Reference Hagemann, Rasoulzadeh, Panfilov, Ganzer and Reitenbach2015). This behaviour not only reduces storage efficiency but also increases the risk of hydrogen loss.
The injected layer of hydrogen, through its displacement of heavier brine, is a form of viscous gravity current (Zheng & Stone Reference Zheng and Stone2022). Much of the research on viscous gravity currents has focused on incompressible flows, where the interaction between buoyancy and viscosity governs the rate and extent of spreading. In the context of hydrogen storage, the low viscosity of the injected gas, coupled with buoyancy effects, plays a critical role in determining how far and how quickly the hydrogen spreads underground. Early studies primarily addressed gravity-driven spreading, where a dense fluid advances along a rigid boundary within an infinitely deep porous medium (Barenblatt Reference Barenblatt1952; Huppert & Woods Reference Huppert and Woods1995; Lyle et al. Reference Lyle, Huppert, Hallworth, Bickle and Chadwick2005). However, when a fluid is actively injected into a confined porous layer, the resulting pressure gradient, caused by the pressure drop between the injection site and the far field, significantly alters the flow dynamics.
Huppert & Woods (Reference Huppert and Woods1995) examined the exchange of two confined incompressible fluids with differing densities and pressures, but equal viscosities, between two aquifers. By deriving a similarity solution, they showed that buoyancy transports the denser fluid along the lower boundary and the lighter fluid along the upper boundary, while an imposed pressure gradient drives a net flow. Subsequent research expanded on this framework to incorporate fluids with different viscosities, highlighting how the viscosity contrast influences spreading behaviour (Nordbotten & Celia Reference Nordbotten and Celia2006; Pegler, Huppert & Neufeld Reference Pegler, Huppert and Neufeld2014; Zheng et al. Reference Zheng, Guo, Christov, Celia and Stone2015). Nordbotten & Celia (Reference Nordbotten and Celia2006) generalised Huppert & Woods’ (Reference Huppert and Woods1995) study to investigate the effect of viscosity contrasts in an axisymmetric geometry. Pegler et al. (Reference Pegler, Huppert and Neufeld2014) conducted a comprehensive study of flow behaviour in a two-dimensional channel. Their findings revealed that, at early times, spreading is primarily buoyancy driven and can be approximated by the classical porous-medium equation of Barenblatt (Reference Barenblatt1952). At later times, low-viscosity injected fluid forms a thin layer along the upper boundary, with buoyancy becoming less dominant compared with viscous effects. By deriving a large-time similarity solution, they demonstrated that the movement of the upper contact line is controlled by the strength of the source. In this case, since the flow is incompressible, the length of the channel determines the pressure needed to drive the motion but does not affect the spreading dynamics. These theoretical predictions were validated through laboratory experiments in which freshwater was injected into a saltwater-saturated porous medium, confirming that the similarity solution accurately describes the interface evolution.
When the injected gas is compressible, however, the dynamics change: the source can in principle increase the gas pressure without advancing the gas–liquid interface, and the subsequent propagation of the upper contact line may be influenced by the length of the channel. In long channels, the viscous pressure drop between source and outlet may lead to spatial variations in pressure, which could in turn lead to noticeable variations in gas density along the channel. These variations could affect the mass flow rate and, consequently, influence the rate of gas spreading. Despite the potential for compressibility to significantly alter the flow dynamics, its effects in this context remain under-explored. While some studies have addressed compressibility, most have focused on fluid–structure interactions, particularly the impact of pressure buildup on the rock matrix integrity (Mathias et al. Reference Mathias, Hardisty, Trudell and Zimmerman2009) and fluid migration in neighbouring reservoirs (Jenkins, Foschi & MacMinn Reference Jenkins, Foschi and MacMinn2019). A common approach to coupling the mechanics of the pore structure with pressure buildup during injection is to introduce rock compressibility, defined as
$c_r = ({1}/{\phi }) {\mathrm{d} \phi }/{\mathrm{d} p^*}$
, where
$\phi$
is the porosity and
$p^*$
is the pressure. This parameter quantifies how the pore space deforms in response to pressure changes and is valid under conditions of constant vertical stress and negligible lateral strain (Jenkins et al. Reference Jenkins, Foschi and MacMinn2019).
Mathias et al. (Reference Mathias, Hardisty, Trudell and Zimmerman2009) examined pressure buildup during gas injection into an axisymmetric channel of infinite extent. Using matched asymptotics, they derived similarity solutions for pressure evolution in the case where the compressibilities of the gas and resident liquid are comparable (relative to rock and water), and where the interface propagates much more slowly than the diffusive pressure front. This work was later extended by Mathias et al. (Reference Mathias, González Martínez De Miguel, Thatcher and Zimmerman2011) to finite-length channels, offering further insight into pressure evolution under more realistic conditions. Jenkins et al. (Reference Jenkins, Foschi and MacMinn2019) investigated pressure dissipation in a layered aquifer system, where multiple confined channels are stacked vertically. They analysed how water leakage from the injection aquifer to the surrounding aquifers influences the pressure field. In their model, CO
$_2$
displaces resident water in the injection channel, and weak vertical flow allows liquid to escape through the confining boundaries, relieving pressure within the system. By incorporating the compressibility of both gas and liquid through a linearised equation of state, they found that pressure reduction due to water escape slows the gas plume’s spreading rate. In contrast to these studies, Cuttle, Morrow & MacMinn (Reference Cuttle, Morrow and MacMinn2023) examined how compressibility interacts with surface tension and viscous forces in the liquid, to influence viscous fingering instabilities, when a compressible gas is injected into a Hele-Shaw cell filled with water. By neglecting pressure variations in the gas due to viscous effects, they used Boyle’s law to couple the uniform pressure in the gas, which drives displacement, to the viscous effects in the liquid and the capillary forces at the interface. Their results showed that gas compressibility delays the onset of viscous fingering and reduces its severity, with the severity quantified by the isoperimetric ratio (which compares the length of the interface to the area it occupies). Additionally, compressibility was found to increase the breakthrough time, i.e. the moment when the interface reaches the outlet.
In this paper we extend the framework of Pegler et al. (Reference Pegler, Huppert and Neufeld2014) for injection of a gas into a confined porous channel containing brine (figure 1) to investigate how gas compressibility interacts with buoyancy and viscous effects in both fluids to determine the spreading dynamics. Using a long-wave approximation, we derive two coupled nonlinear evolution equations governing the gas/liquid interface height and the gas pressure field. These are mass conservation equations, with the gas pressure serving as a proxy for gas density via the equation of state. The viscous pressure drop ahead of the upper contact line is captured by a boundary condition, while two kinematic conditions determine the motion of the contact lines. The model admits several distinguished asymptotic limits that reveal the roles of different mechanisms across parameter space; we derive a set of reduced models that capture dominant balances. In particular, we estimate key quantities of practical interest, including the pressure scale, pressure rise time and the breakthrough time. In § 2 we present the governing equations and describe the scaling argument used to derive the long-wave model. The system is characterised by three primary dimensionless parameters, and we simplify the model in various sublimits, identifying relevant scaling relationships. In § 3 we present numerical solutions of the full model, illustrating how the key features of compressible flow emerge and comparing them with predictions from asymptotic analysis. Finally, in § 4 we revisit the dominant balances in each regime, interpret the underlying physical mechanisms and discuss the broader implications for hydrogen storage and related gas-injection technologies.
Figure 1. Schematic showing the displacement of an ambient liquid in a porous medium (regions II
$_w$
and III) due to the injection of a compressible gas (occupying regions I and II
$_g$
) from a line source at
$\boldsymbol{x}^* = \boldsymbol{x}_0^*$
. The interface separating the fluids is sharp and located at
${y}^* = {F}^*({x}^*, {t}^*)$
. The pressure in the ambient brine is hydrostatic at the outlet at
$x^*=L$
.
2. Model formulation
2.1. Governing equations
We consider a two-dimensional model of fluid displacement driven by the injection of a compressible gas into a planar horizontally confined porous medium with height
$H$
and length
$L$
, as illustrated in figure 1. Initially, the medium is predominantly saturated with a liquid, with only a small amount of gas present. The length of the region initially occupied by gas is denoted by
$L_{g0}$
; the aspect ratio of this region is assumed small, i.e.
$H/L_{g0} \ll 1$
. (We do not seek here to model the initial formation of the gas bubble.) The permeability of the medium is taken to be uniform, with value
$k_0$
.
At time
$t^* = 0$
, gas is injected into the channel from a line source located at
$\boldsymbol{x}^* \,\equiv (x^*,y^*) = \boldsymbol{x}^*_0$
. Variations in the injection rate are governed by the function
$q \, \mathcal{Q}(\omega {t^*})$
, where
$1/\omega$
represents the time scale of injection,
$q$
denotes the source strength per unit width of the channel and
$\mathcal{Q}$
is a dimensionless function. The pressure at the outlet (
$x^* = L$
) is taken to be hydrostatic with a baseline value of
$p_{g0}$
. Using a sharp-interface model, disregarding capillary or miscibility effects, the fluids are segregated into two distinct phases and the domain is divided into four regions (figure 1): the gas occupies regions I and II
$_g$
, while the liquid occupies regions II
$_w$
and III. We use the subscripts
$g$
and
$w$
to represent quantities in the gas and liquid (water) phase, respectively. We revisit all modelling assumptions in § 4.
Figure 2. (a) Isothermal equations of state for hydrogen, methane, nitrogen and carbon dioxide at 333 K, sourced from Lemmon et al. (Reference Lemmon, Bell, Huber, McLinden, Linstrom and Mallard2023). (b) The dimensionless function
$\mathcal{P}$
(2.2) for hydrogen (solid line), computed using the reference pressure
$p_{g0} = 23$
MPa and reference density
$\rho _{g0} = 15$
$\mathrm{kg\,m^{-3}}$
. The dashed line represents the tangent at the reference state, with slope
$\mathcal{P}^{\prime }(1)$
.
We denote the height of the gas–liquid interface by
$y^* = F^*(x^*, t^*)$
and the points of contact with the upper and lower boundaries by
$X^*_u({t^*})$
and
$X^*_l({t^*})$
, respectively. For later convenience, we extend the definition of
$F^*$
to the whole domain such that
\begin{align} F^*(x^*,t^*)=\begin{cases} 0, & 0\leqslant x^*\leqslant X^*_l(t^*), \\ H, & X^*_u(t^*) \leqslant x^* \leqslant L. \end{cases} \end{align}
The liquid is assumed to be incompressible with constant density
${\rho }_w$
. The initial gas and liquid volumes are assumed to be in thermal equilibrium with the surrounding pore matrix. Gas compression is treated as isothermal, an assumption that is valid when the thermal inertia of the solid pore matrix is much greater than that of the gas (Kushnir et al. Reference Kushnir, Dayan and Ullmann2012a
). At a fixed temperature, we define
$\rho _{g0}$
to be the equilibrium gas density when at the outlet pressure
$p_{g0}$
. Then, using
$p_{g0}$
and
$\rho _{g0}$
as a reference pressure and a reference density, we write the equation of state for the gas (relating density
$\rho ^*_g$
to pressure
$p^*_g$
) and the equilibrium speed of sound in the gas (
$c$
) as
\begin{align} p^*_g = p_{g0} \mathcal{P}\left (\frac {\rho ^*_g}{\rho _{g0}}\right ), \qquad c^2 \equiv \left . \frac {\partial {p^*_g}}{\partial \rho ^*_g} \right |_{\rho _{g0}} = \frac {p_{g0}}{\rho _{g0}} \mathcal{P}^{\prime }(1), \end{align}
respectively, for a suitable dimensionless function
$\mathcal{P}$
, for which
$\mathcal{P}(1)=1$
and
$\mathcal{P}'(1)\gt 0$
. For hydrogen,
$\mathcal{P}'(1)\geqslant 1$
and
$\mathcal P''(1)\geqslant 0$
at temperatures and pressures of interest (figure 2
$b$
), so that
$c^2$
rises with increasing pressure, although this may not be the case for gases such as CO
$_2$
(figure 2
$a$
).
The fluid velocities
$\boldsymbol{u}_i^* = (u^*_i, v^*_i)$
in both phases are then given by
\begin{align} \boldsymbol{u}^*_i = -\frac {k_0}{\mu _i} \: ( \boldsymbol{\nabla }^* p^*_i - \rho ^*_i \boldsymbol{g}) \quad (i = w,g), \end{align}
where
$\boldsymbol{g}$
is the acceleration due to gravity,
$\rho _w^*\equiv \rho _w$
and
$\mu _g$
and
$\mu _w$
are gas and liquid viscosities, assumed constant. The pore-scale Reynolds number,
$Re = q d / (H\mu _g)$
, is assumed sufficiently small to neglect inertia, where
$d$
is a typical pore diameter. Taking
$H = 10\, \mathrm{m}$
,
$\mu _g\approx 10^{-6}\,\mathrm{kg\,m^{-1}\,s^{-1}}$
and
$d\approx 10\,\unicode{x03BC} \mathrm{m}$
(Hashemi, Blunt & Hajibeygi Reference Hashemi, Blunt and Hajibeygi2021), the condition
$Re\leqslant 1$
suggests that inertia can be neglected for injection rates up to
$ q \approx 1\, \mathrm{kg\,m^{-1}\,s^{-1}}$
.
The continuity equations for the gas and liquid phases are
\begin{align} \rho ^*_{g, t^*} + {\nabla }^* \boldsymbol{\cdot }(\rho ^*_g\boldsymbol{u}^*_g) = q \,\mathcal{Q}(\omega t^*) \, \delta ^*(\boldsymbol{x}^* - \boldsymbol{x}_0^*), \quad{\nabla }^* \boldsymbol{\cdot }\boldsymbol{u}_w^* = 0, \end{align}
respectively, where
$\delta ^*(\boldsymbol{x}^*)$
is the Dirac delta function. The source at
$\boldsymbol{x}^* = \boldsymbol{x}_0^*$
lies within a distance of order
$H$
of an impermeable boundary at
$x^*=0$
. Equations (2.2), (2.3) and (2.4a
) form the governing equations for the gas and are applied in regions I and II
$_g$
; (2.3) and (2.4b
) apply in regions II
$_w$
and III. The kinematic condition
\begin{align} F^*_{t^*} + u^*_i F^*_{x^*} = v^*_i\quad \text{at} \ y^* = F^* \quad (i = w,g; \,0\leqslant x^*\leqslant L), \end{align}
relates the velocity of the interface to the fluid velocities for
$X^*_l\lt x^*\lt X^*_u$
. Subscripts
$x^*$
and
$t^*$
denote partial derivatives. No penetration at the walls implies that
$u^*_g = 0$
at
$x^* = 0$
,
$v^*_i = 0$
at
$y^* = H$
and
$v^*_i = 0$
at
$y^* = 0$
. The relative strength of buoyancy to capillary forces can be characterised by a Bond number. We adopt the definition of the Bond number from Golding et al. (Reference Golding, Neufeld, Hesse and Huppert2011) and Zheng & Neufeld (Reference Zheng and Neufeld2019), which measures the relative magnitude of buoyancy to the capillary entry pressure
$p_e$
in the largest pores. Using
$p_e = \gamma / k_0^{1/2}$
(Zheng & Neufeld Reference Zheng and Neufeld2019), the Bond number can be expressed as
$Bo = \Delta \rho g H k_0^{1/2}/\gamma$
, where
$\Delta \rho = \rho _w - \rho _{g0} \approx 10^3\,\mathrm{kg\,m^{-3}}$
is the density difference between the fluids,
$g$
is gravitational acceleration,
$H \approx 10$
–
$50\,\mathrm{m}$
is a characteristic aquifer height,
$k_0 \approx 10^{-14}$
–
$10^{-12}\,\mathrm{m^2}$
is the permeability and
$\gamma \approx 7 \times 10^{-2}\,\mathrm{N\,m^{-1}}$
is the hydrogen–liquid interfacial tension (Zivar et al. Reference Zivar, Kumar and Foroozesh2021). An order-of-magnitude estimate gives
$Bo \sim 0.1\,{-}\,10$
, indicating that capillary effects may be important in shallow, low-permeability reservoirs. A detailed treatment of the capillary fringe is beyond the scope of this study, although we discuss its potential influence in § 4. Accordingly, we assume a sharp interface and take the pressure to be continuous across it:
\begin{align} p^*_g = p^*_w \quad \text{at} \ y^* = F^* \ \text{for}\ X^*_l\lt x^*\lt X^*_u. \end{align}
The liquid pressure at the outlet is
\begin{align} p^*_w = p_{g0} - \rho _w g y^* \quad \text{at} \ x^* = L. \end{align}
The initial bubble length specifies the initial gas volume (per unit width) via
$V^*(0)=\int _0^{X^*_u} (H-F^*(x^*,0))\,\mathrm{d}x^*=HL_{g0}$
. The interface is assumed to have a linear initial profile over a length scale
$\Delta _{g0}$
, so that
\begin{align} F^*(x^*,0)=\frac {H}{{\Delta }_{g0}}(x^*-X_l^*(0))\,\mathrm{for}\, X_l^*(0)\equiv L_{g0}-\frac {1}{2}{\Delta }_{g0}\lt x^*\lt X_u^*(0)\equiv L_{g0}+\frac {1}{2}{\Delta }_{g0}. \end{align}
We integrate (2.4a
) over
$F^*\lt y^*\lt H$
for
$0\leqslant x^*\lt X^*_u$
and (2.4b
) across
$0\lt y^*\lt F^*$
for
$X^*_u\lt x\lt L$
, which, after applying (2.5) and the no-penetration conditions, respectively, yield the transport equations
\begin{align} \left ( \int _{F^*}^H \rho ^*_g \: \mathrm{d} y^* \right )_{t^*} + \left ( \int _{F^*}^H \rho ^*_g u^*_g \:\mathrm{d} y^* \right )_{x^*} &= q \mathcal{Q}(\omega t^*) \delta ^*(x^* - x^*_0) \quad (0 \lt x^* \lt X^*_u), \end{align}
\begin{align} {F^*}_{t^*} + \left (\int _0^{F^*} u^*_w \: \mathrm{d} y^* \right )_{x^*}& = 0 \quad (X^*_l \lt x^* \lt L). \end{align}
The total masses (per unit width) of gas and liquid in the channel at a given instant are then expressed using (2.1) as
\begin{align} M^*_g(t^*) & = \int _{0} ^{X^*_u(t^*)} \int _{F^*(x^*,t^*)}^H \rho ^*_g(x^*,y^*,t^*) \: \mathrm{d} y^* \mathrm{d} x^* , \\[-9pt] \nonumber \end{align}
\begin{align} M^*_w(t^*) &= \int _{X^*_l(t^*)}^{L} \int _0^{F^*(x^*,t^*)} \rho ^*_w \: \mathrm{d} y^* \mathrm{d} x^* , \end{align}
respectively. Integration of (2.9) along the channel, and using the no-penetration condition on the wall at
$x^* = 0$
, gives the rates of gas accumulation in the channel and liquid outflow as
\begin{align} \frac {\mathrm{d} M^*_g}{\mathrm{d} {t^*}} = q \,\mathcal{Q}(\omega t^*) ,\quad \frac {\mathrm{d} M^*_w}{\mathrm{d} {t^*}} = - \rho ^*_w \int _0^H u^*_w(L,y^*,t^*) \mathrm{d} y^*. \end{align}
Equation (2.11a
) shows that the change in mass of gas within the system is due to the influx from the source; (2.11b
) shows that the change in mass of the liquid is equal to the negative of the flux at the outlet. Gas compressibility allows
$\rho ^*_g$
to vary in space and time; however, in the incompressible limit, when
$\rho ^*_g$
is effectively constant, the fixed domain volume ensures that
\begin{align} \frac {{d}}{\mathrm{d}t^*} \left ( \frac {M^*_g}{\rho ^*_g}+\frac {M^*_w}{\rho _w}\right )=0, \end{align}
relating the source strength to the outflow directly via (2.11). In general, however, (2.12) will not hold while gas is compressed transiently in regions I and II
$_g$
.
Input parameters are summarised in table 1. We seek outputs such as the breakthrough time
$t^*_b$
(at which
$X^*_u(t^*_b)=L$
) and the total mass of gas delivered at this time
$M_g^*(t_b^*)$
. (The model applies only for
$0\lt t^*\leqslant t_b^*$
.) We now present an asymptotic reduction of (2.9) to a coupled system of evolution equations (see (2.29a
,
b
)) for the interface position
$F^*$
and the gas pressure
$p^*_g$
, beginning with a scaling argument.
Table 1. Parameters used in the governing equations, their approximate values and corresponding sources.
2.2. Scaling argument
The gas mass flux per unit width from the source
$q$
generates a volume flux per unit length of order
$q/\rho _{g0}$
and a horizontal velocity of order
$q/(\rho _{g0} H)$
. The corresponding transit time along the channel is of order
$\rho _{g0} H L/q$
. Assuming for the time being that
$\mu _g$
and
$\mu _w$
are comparable, the viscous pressure drop along the channel is of order
$\mu _g q L/(k_0 \rho _{g0} H)$
. Vertical hydrostatic pressure variations in each phase are of order
$\rho _{g0} g H$
and
$\rho _w gH$
, with the difference
$\Delta \rho g H$
providing the buoyancy force that seeks to flatten the gas–liquid interface. In the long-wave limit, vertical velocities are expected to be a factor
$H/L_{g0}$
smaller than the horizontal velocity, allowing the dominant contributions to the gas and liquid pressure fields to be hydrostatic plus a field that varies with
${x}^*$
and
${t}^*$
to leading order, independent of
${y}^*$
. Assuming that
$\rho _{g0} \ll \rho _w\approx \Delta \rho$
, the hydrostatic field in the liquid dominates that in the gas.
Suppose that the buoyant and compressible pressure variations are comparable to
$p_{g0}$
, but that the viscous pressure scale is larger, so that
\begin{align} \frac {\mu _g}{\mu _w} {\,=\,} O(1),\quad \frac {q \mu _g}{k_0 p_{g0} \rho _{g0}}\frac {L}{H}\gg 1, \quad \frac {\rho _{g0} g H}{p_{g0}}\ll \frac {\Delta \rho g H}{p_{g0}} {\,= \,} O(1), \quad \frac {\rho _{g0} c^2}{p_{g0}} \gtrsim 1. \end{align}
We assume that the equation of state (2.2) can be approximated by the truncated Taylor series
$\mathcal{P}(1+\vartheta )\approx 1+\vartheta \mathcal{P}'(1)+( {1}/{2})\vartheta ^2\mathcal{P}''(1)$
, with
$\mathcal{P}'(1)$
and
$\mathcal{P}''(1)$
both being of order unity. (Figure 2b
shows that, for hydrogen, we can assume that
$\mathcal{P}'(1)\approx 1$
and
$\mathcal{P}''(1)\ll 1$
.) Thus, relative density variations in (2.2),
$\rho ^*_g/\rho _{g0}=1+\vartheta$
, generate gas pressure variations of order
$p_{g0}$
. For
$\vartheta \ll 1$
, for example, (2.2) can be approximated by the linear relationship
\begin{align} {p^*_g\approx p_{g0}+c^2(\rho ^*_g-\rho _{g0})}. \end{align}
We then use (2.14) to write the transport equation for gas density (2.9) as an evolution equation for gas pressure. Anticipating no vertical gradients in
$p^*_g$
to leading order (the hydrostatic component being subdominant), so that
$p_g^*=p_g^*(x^*,t^*)$
, the transported quantity
$\int \rho _g^* \,\mathrm{d}y^*$
in (2.9a
) is proportional to
$(H-F^*)(\rho _{g0} + (p^*_g-p_{g0})/c^2)$
. The time derivative (in (2.9a
)) then includes
$-\rho _{g0}F^*_{t^*}$
and
$Hp^*_{g,t^*}/c^2$
.
Eliminating
$p_{g0}$
from ratios in the second, third and fourth terms of (2.13) identifies the dimensionless parameters
\begin{align} \mathcal{M}=\frac {\mu _g}{\mu _w}\lesssim 1, \quad \mathcal{L} \equiv \frac {q \mu _g L}{k_0 \rho _{g0} \Delta \rho g H^2}\gg 1,\quad \zeta \equiv \frac {\Delta \rho g H}{\rho _{g0} c^2} \lesssim 1, \quad \beta =\frac {p_{g0}}{\rho _{g0} c^2}\equiv \frac {1}{\mathcal{P}'(1)}\lesssim 1. \end{align}
We take the viscosity ratio
$\mathcal{M}$
to be of order unity for now, specialising later to the limit
$\mathcal{M}\ll 1$
. Buoyancy-driven flattening of the interface takes place over a length scale
$L_b$
for which the pressure gradient
$\Delta \rho g H/L_b$
generates a horizontal velocity
$k_0 \Delta \rho g H/(\mu _g L_b)$
that we assume balances the source-driven velocity
$q/(\rho _{g0} H)$
. It follows that
\begin{align} L_b=L/\mathcal{L} \ll L. \end{align}
Thus, we can interpret
$\mathcal{L}$
as a measure of the channel length relative to the length scale over which buoyancy forces flatten the interface. The corresponding time scale for buoyancy-driven flow is
$T_b=L_b \rho _{g0} H/q$
, which is smaller than the transit time by a factor
$1/\mathcal{L} \ll 1$
. Taking the pressure scale
$\Delta \rho g H$
, the two time derivatives in the gas transport equation (
$-\rho _{g0}F^*_{t^*}$
and
$Hp^*_{g,t^*}/c^2$
) differ in magnitude by
$\zeta$
, a parameter that measures hydrostatic pressure changes relative to those associated with gas compressibility. Thus, for
$\zeta \lesssim O(1)$
with
$\mathcal{M}=O(1)$
, we expect compressible effects to operate on shorter time scales than source- and buoyancy-driven spreading, which in turn may take a long time to reach the end of the channel. The parameter
$\beta$
in (2.15) measures the global convexity of the equation of state at the baseline state: as illustrated in figure 2(b) for hydrogen,
$\beta$
is slightly below unity because the equation of state is almost linear but with a small positive curvature.
To model these multiple physical effects, we therefore scale
$x^*$
by
$L_b$
,
$y^*$
and
$F^*$
by
$H$
,
$t^*$
by
$T_b$
,
$p^*_g$
by
$\Delta \rho g H$
and
$\rho ^*_g$
by
$\rho _{g0}$
. The dimensionless domain length is then
$\mathcal{L}$
. Assuming that
$\mathcal{L}\gg 1$
allows us to track the evolution of the interface over long distances and long times. At the early stages of spreading, when compressible effects are likely most dominant, the dimensionless viscously generated gas pressure is of order
$\mathcal{L}$
, which yields dimensionless density variations (
$\vartheta$
) of order
$\zeta \mathcal{L}$
. To capture compressible effects, we seek to accommodate the distinguished limit
\begin{align} \zeta \rightarrow 0, \quad \mathcal{L}\rightarrow \infty \quad \mathrm{with}\quad \theta \equiv \, \zeta \mathcal{L} \equiv \frac {q \mu _g L}{k_0 \rho _{g0}^2 H c^2} {\,=\,} O(1), \end{align}
where
$\theta$
measures the relative magnitude of viscously generated pressure compared with that due to compressibility; this is analogous to the compressibility number defined in Cuttle et al. (Reference Cuttle, Morrow and MacMinn2023). Here the long domain length (
$\mathcal{L}\gg 1$
) generates large pressure deviations that cause appreciable density variations, even in a weakly compressible gas (
$\zeta \ll 1$
).
Shortly, we extend the model to incorporate the limit
$\mathcal{M}\ll 1$
, representing the motion of a gas into a liquid. In this case, a thin film of gas can spread rapidly over the liquid, reducing the time taken for the gas to travel to the channel outlet. We show how compressible effects can remain dominant throughout this spreading process. For now, we proceed assuming that
$\mathcal{M}=O(1)$
.
2.3. Model equations
2.3.1. Non-dimensionalisation
Adopting the proposed scalings (
$x^*=L_bx$
,
$y^*=Hy$
,
$t^*=T_bt$
,
$\rho ^*_g=\rho _{g0} \rho _{g}$
,
$p^*_g=(\Delta \rho g H) p_g$
,
$p^*_w=\Delta \rho g H(p_w-y)$
,
$V^*=HL_bV$
,
$M_g^*=\rho _{g0}HL_b M_g$
),
$(u_g^*,v_g^*)=(q/H\rho _{g0})(u_g,\epsilon v_g)$
,
$(u_w^*,v_w^*)=(q/H\rho _{g0})(u_w,\epsilon v_w)$
with
$\epsilon \equiv H/L_b\ll 1$
, the Darcy equations (2.3) become
\begin{align} u_g & =-p_{g,x},\qquad\qquad\qquad\qquad u_w=-\mathcal{M} p_{w,x}, \end{align}
\begin{align} \epsilon ^2 v_g & =-p_{g,y}+(\rho _{g0}/\Delta \rho ),\qquad \epsilon ^2 v_w=\mathcal{M}\left [-p_{w,y}+(\rho _{g0}/\Delta \rho )\right ]. \end{align}
Thus, assuming that
$\rho _{g0}\ll \rho _w$
, ensuring that
$\rho _{g0}\ll \Delta \rho$
,
$p_g$
and
$p_w$
are independent of
$y$
at leading order. We also assume, for simplicity, that
$\beta \mathcal{P}''(1) \ll 1$
: with a linearised equation of state (2.14), the gas density
$\rho _g$
is represented in terms of gas pressure by
$1-\beta +\zeta p_g$
. With the flow configured as in figure 1, we then recover from (2.6), (2.7), (2.9):
\begin{align} \left [(1-F)(1-\beta +\zeta p_g )\right ]_t& = \left [ (1-F)(1-\beta +\zeta p_g)p_{g,x}\right ]_x \quad(0\lt x\lt X_u),\\[-10pt] \nonumber \end{align}
\begin{align} F_t& = \mathcal{M} \left [F p_{w,x} \right ]_x\quad (X_l\lt x\lt X_u), \\[-10pt] \nonumber \end{align}
\begin{align} 0 &= \mathcal{M} \left [p_{w,x}\right ]_x\quad (X_u\lt x\lt \mathcal{L}), \\[-10pt] \nonumber \end{align}
\begin{align} - (1-\beta +\zeta p_g)p_{g,x}&=\mathcal{Q}\quad (x=0), \\[-10pt] \nonumber \end{align}
\begin{align} p_w &=p_g+F\quad (X_l\lt x\leqslant X_u), \hphantom{fffffffffffff} \\[-10pt] \nonumber \end{align}
\begin{align} (p_{g,x}+F_x)\big \vert _{X_u-}&=p_{w,x}\big \vert _{X_u+}, \\[-10pt] \nonumber \end{align}
\begin{align} p_w & =\beta /\zeta \quad (x=\mathcal{L}, y=0). \\[9pt] \nonumber \end{align}
The problem is governed by two coupled evolution equations (2.19a
,
b
), unlike the incompressible problem in which a single evolution equation describes the dynamics. These equations describe respectively conservation of mass in the gas (of thickness
$1-F$
) and of the liquid beneath it (of thickness
$F$
). Equation (2.19c
) is a statement of mass conservation where the liquid fully fills the channel; (2.19d
) balances the mass flux of gas near the inlet (assuming that
$F(0,t)=0$
) with the imposed source; (2.19e
) ensures continuity of pressure across the gas–liquid interface; (2.19f
) ensures continuity of liquid flux across the upper contact line; and (2.19g
) imposes the hydrostatic pressure constraint at the channel outlet. Integrating (2.19c
) to find
$p_w$
, and applying (2.19e
,
f
) at
$x=X_u$
and (2.19g
), gives
\begin{align} p_g+ (p_{g,x}+F_x) (\mathcal{L} - X_u)=({\beta }/{\zeta }) -1 \quad (x=X_u-). \end{align}
Referring to figure 1, two boundary conditions are required in region I (for gas alone), two boundary conditions are required in region III (for liquid alone) and four are required for region II (for gas and liquid). Thus, we have a single inlet condition (2.19d
), a single outlet condition (2.19g
) plus continuity and kinematic conditions at internal free boundaries. Specifically, combining (2.19a
,
b
) with the constraints
$F(X_l,t)=0$
,
$F(X_u,t)=1$
requires that
\begin{align} X_{l,t}=-\mathcal{M} (p_{g,x}+F_x)\big \vert _{X_l+}, X_{u,t}=- p_{g,x}\big \vert _{X_u-}. \end{align}
We define the interface length as
$\Delta (t) \equiv X_u(t) - X_l(t)$
, with initial value
$\Delta (0)=\Delta _0\,\equiv {\Delta }_{g0}/L_b$
. The initial condition (2.8) becomes
\begin{align} F(x,0) \equiv \begin{cases}0, & 0\leqslant x\lt X_l(0)\equiv L_0-\dfrac {1}{2}\Delta _0, \\(x-X_l(0))/\Delta _0, & X_l(0)\leqslant x\lt X_u(0), \\ 1, & L_0+\dfrac {1}{2}\Delta _0\equiv X_u(0)\leqslant x\leqslant \mathcal{L},\end{cases} \end{align}
where
$L_0\equiv L_{g0}/L_b$
. An initial condition is also required for
$p_g$
; we assume for now that it is large enough in magnitude for the gas volume to increase. The volume of gas in the channel (per unit width) is given by the integral
\begin{align} V(t) = \int _0^{X_u(t)} (1 - F) \, \mathrm{d} x; \end{align}
thus, (2.22) defines the initial gas volume
$V(0)=( {1}/{2}) [X_u(0) + X_l(0)]=L_0$
. The mass of gas (2.11a
) evolves according to
\begin{align} \frac {\mathrm{d}M_g}{\mathrm{d}t}=\mathcal{Q}(\varOmega t) \quad \mathrm{where}\ \varOmega =\omega T_b. \end{align}
The full dimensionless problem is governed by seven dimensionless parameters (
$\beta$
,
$\varOmega$
,
$L_0$
,
$\Delta _0$
,
$\mathcal{L}$
,
$\mathcal{M}$
,
$\zeta$
), which are summarised in table 2. We show shortly how the convexity parameter
$\beta$
can be eliminated and focus attention primarily on steady injection with
$\mathcal{Q}=1$
, although we retain
$\mathcal{Q}$
in the following in order to comment later on time-varying injection rates parametrised by
$\varOmega$
. In many instances we expect details of the initial conditions
$L_0$
,
$\Delta _0$
to be unimportant at large times. This leaves the domain length
$\mathcal{L}$
, viscosity ratio
$\mathcal{M}$
and compressibility
$\zeta$
as parameters of primary interest.
Table 2. A summary of the seven dimensionless parameters in the governing equations (2.19).
2.3.2. Incompressible limit
The incompressible limit of (2.19) with
$\mathcal{L}=O(1)$
is recovered by taking
$\zeta \rightarrow 0$
,
$\beta \rightarrow 0$
with
$\beta /\zeta =p_{g0}/(\Delta \rho g H) \,{=\,} O(1)$
. Using (2.19e
) to eliminate
$p_w$
, the two evolution equations (2.19a
,
b
) become
\begin{align} -F_t &= \left [ (1-F)p_{g,x}\right ]_x \quad (0\lt x\lt X_u), \end{align}
\begin{align} F_t& = \mathcal{M} \left [F (p_{g,x}+F_x)\right ]_x \quad (X_l\lt x\lt X_u), \end{align}
which are subject to the boundary conditions (2.19d ), (2.20), (2.21), from which the flux constraint
\begin{align} (1-F)p_{g,x}+\mathcal{M}F(p_{g,x}+F_x)=-{\mathcal{Q}} \quad (0\lt x\lt X_u) \end{align}
emerges. Elimination of
$p_{g,x}$
between (2.25a
), (2.26) yields a single evolution equation for
$F$
,
\begin{align} F_t=\mathcal{M}\left [\frac {F(1-F)F_x-F\mathcal{Q}}{1-F+\mathcal{M}F}\right ]_x \quad (X_l\lt x\lt X_u), \end{align}
a system addressed by Pegler et al. (Reference Pegler, Huppert and Neufeld2014) and others (with
$\mathcal{Q}=1$
). The term proportional to
$F_x$
captures the role of buoyancy; the term proportional to
$\mathcal{Q}$
represents advection driven by injection. The source strength
$\mathcal{Q}$
appears only in a boundary condition in the compressible problem (2.19), but incompressibility enables its presence to be felt throughout the domain in (2.27). Assuming the (diffusive) buoyancy term is subdominant to (advective) viscous effects, so that
$F_t+\mathcal{M}\mathcal{Q} (1-F+\mathcal{M}F)^{-2}F_x=0$
, a solution using characteristics shows how the source drives both contact lines directly:
\begin{align} X_l(t)= X_l(0) + \mathcal{M} \int _0^t \mathcal{Q}\,\mathrm{d}t, \quad X_u(t)= X_u(0) + \frac {1}{\mathcal{M}}\int _0^t\mathcal{Q}\,\mathrm{d}t. \end{align}
For
$\mathcal{Q}=1$
, this becomes, in the large-time limit (Pegler et al. Reference Pegler, Huppert and Neufeld2014),
\begin{align} X_l(t)\approx {\mathcal{M} t}, \quad X_u(t)\approx {t}/{\mathcal{M}}. \end{align}
In this formulation,
$t=\mathcal{L}$
corresponds to the transit time of source-driven flow. The breakthrough time of
$X_u$
at the outlet predicted by (2.28b
),
$t_b=\mathcal{M}\mathcal{L}$
, lies below the transit time for
$\mathcal{M}\lt 1$
, because a thin film of gas spreads rapidly over the top of the liquid.
2.3.3. Compact formulation of the model
A more compact formulation of (2.19), that is independent of
$\beta$
, is obtained by writing the gas density as
$\zeta P=1-\beta +\zeta p_g$
(or
$p^*_g=p_{g0}-\rho _{g0} c^2+(\Delta \rho g H) P$
, so that the equation of state (2.14) reduces to
$\zeta P =\rho _{g}$
), and (2.19) simplifies to
\begin{align} \left [(1-F)P\right ]_t &= \left [ (1-F) P P_{x}\right ]_x \quad (0\lt x\lt X_u), \\[-9pt] \nonumber \end{align}
\begin{align} F_t& = \mathcal{M} \left [F(P_{x}+F_x)\right ]_x \quad (X_l\lt x\lt X_u), \\[-9pt] \nonumber \end{align}
\begin{align} - \zeta P P_x & =\mathcal{Q} \quad (x=0), \\[-9pt] \nonumber \end{align}
\begin{align} P+(P_{x}+F_x) (\mathcal{L} - X_u)& =(1/\zeta )-1 \quad (x=X_u), \\[-9pt] \nonumber \end{align}
\begin{align} X_{l,t}& =-\mathcal{M} (P_{x}+F_x) \quad (x={X_l+}), \\[-9pt] \nonumber \end{align}
\begin{align} X_{u,t}& =- P_{x} \quad (x={X_u-}). \\[9pt] \nonumber \end{align}
Taking
$\mathcal{Q}(0) = 1$
, the initial pressure field is chosen to be
\begin{align} P(x,0) = P_c - \frac {1}{ \zeta P_c} x \quad (0 \lt x \lt X_u), \end{align}
where the constant
$P_c$
satisfies the quadratic equation
\begin{align} P_c^2 + P_c \left ( \frac {\mathcal{L} - X_u(0)}{\Delta _0} + 1 - \frac {1}{\zeta } \right ) - \frac {\mathcal{L}}{\zeta } =0, \end{align}
where
$X_u(0)=L_0+({1}/{2})\Delta _0$
. Equation (2.31) results from substituting (2.30) into (2.29d
) so that the initial pressure field is consistent with the boundary conditions. The initial interface slope
$F_x(x,0) = 1/\Delta _0$
in (2.22) must be consistent with the inequality
$P_x \equiv -1/\zeta P_c \lt (1/\zeta - 1)/(\mathcal{L} - X_u(0)) - 1/\Delta _0$
to ensure compatibility with the boundary condition (2.29d
). Otherwise, the initial condition (2.30) would lead to negative densities (
$P\lt 0$
), rendering (2.29a
) ill-posed.
Having eliminated
$p_w$
and
$\beta$
from the model, we can interpret (2.29a
) as mass conservation for gas, (2.29b
) as mass conservation for liquid, (2.29c
) as the mass flux driving the flow, (2.29d
) as a pressure condition at the upper contact line (capturing the viscous resistance of the liquid-filled region of the aquifer) and (2.29e
,
f
) as kinematic conditions for the contact-line locations where
$F=0$
and
$F=1$
, respectively.
If, at any point during the flow, the slope of the interface at the lower contact line exceeds the magnitude of the pressure gradient at that location, buoyancy will cause the lower contact line to move backwards through (2.29e
). As a result, the lower contact line may reach the boundary at
$x = 0$
before the upper contact line reaches the outlet. To accommodate this, once
$X_l = 0$
we replace the boundary condition (2.29c
) with
\begin{align} \zeta (1-F) P P_x = -\mathcal{Q}, \qquad F_x = -P_x \quad (x = 0). \end{align}
This modification sets the liquid flux at the origin to zero and reduces the effective area through which the gas flux enters the channel by a factor of
$1 - F$
.
It can be verified from (2.29) that the mass of gas
\begin{align} M_g(t)=\zeta \int _0^{X_u}(1-F)P\,\mathrm{d}x \end{align}
satisfies (2.24). In particular, when
$\mathcal{Q}=1$
, this implies that
$M_g(t_b)=M_g(0)+t_b$
, i.e. the breakthrough time
$t_b$
reveals the total mass of gas delivered by the source.
In (2.29) the incompressible limit (2.25), (2.26) with
$\mathcal{L}=O(1)$
and small
$\zeta$
is recovered using the expansion
$P=(1/\zeta )+\hat {P}+\ldots$
for
$\zeta \ll 1$
, with the large mean pressure arising from the
$1/\zeta$
term in (2.29d
).
Figure 3. A schematic map of
$(\zeta , \mathcal{M})$
-parameter space. The full problem (2.29) is derived for
$\mathcal{M}\sim \mathcal{L} \sim \zeta \sim 1$
. The map specialises to the case
$\mathcal{L}\gg 1$
, focusing on high gas mobility and weak compressibility (
$\mathcal{M}\ll 1$
,
$\zeta \ll 1$
); features of the model in the shaded regions in the map become increasingly distinct as
$\mathcal{L}\rightarrow \infty$
. The distinguished limit
$\mathcal{M}\sim \mathcal{L} \zeta \sim 1$
(problem
$\varPi _a$
) is given by (2.37). An inner/outer structure given by (2.38), (2.39) emerges from this system along
$\mathcal{L}^{-1}\ll \mathcal{L} \zeta \sim \mathcal{M} \ll 1$
(problem
$\varPi _{12}$
). The outer problem simplifies to (2.40) for
$\max (\mathcal{L} \zeta ,\zeta ^{1/2})\ll \mathcal{M} \ll 1$
(shaded pink, region
$\varPi _1$
) and (2.44) for
$\mathcal{L}^{-1}\ll \mathcal{M}\ll \min (1,\mathcal{L} \zeta )$
(shaded blue, region
$\varPi _2$
). Buoyancy effects influence the inner region for
$\mathcal{L} \mathcal{M}^2\lesssim 1$
, allowing
$X_l$
to recede. Problem
$\varPi _b$
arises at
$\mathcal{L}\mathcal{M}^2\sim 1$
and
$\mathcal{L}\zeta \sim \mathcal{M}$
. Ultra-low-viscosity effects emerge via (C6) (problem
$\varPi _c$
) at
$\mathcal{M}\sim \mathcal{L}^{-1}$
and
$\zeta \sim \mathcal{L}^{-2}$
, from which emerge sublimits shown in green (region
$\varPi _3$
;
$\mathcal{LM}\ll 1$
,
$\mathcal{L}^{-2}\ll \zeta \ll 1$
) and orange (region
$\varPi _4$
;
$\mathcal{M}\ll \zeta ^{1/2}$
,
$\mathcal{L}^{2}\zeta \ll 1$
). Scales for the pressure at the origin
$P$
, breakthrough time
$t_b$
and pressure rise time
$t_r$
are shown in blue, green and magenta, respectively, in the coloured regions;
$t_r\ll t_b$
above the line
$\mathcal{L} \mathcal{M} \sim 1$
for
$\mathcal{L}^2 \zeta \gtrsim 1$
. Relative locations within the parameter space of results shown in figures 4–9 are indicated with symbols.
2.3.4. Summary of model formulation
To summarise, we have taken a long-wave limit to derive evolution equations (2.19)–(2.21) for the interface shape and gas pressure, generalising an existing model to incorporate gas compressibility. This system has been reformulated as (2.29), the core model of interest, which shows how the transport involves competing nonlinear diffusion (with fluxes proportional to
$(1-F)PP_x$
and
$FF_x$
) and advection (
$FP_x$
). Numerical solutions are given in § 3, obtained using methods described in Appendix A. The system (2.29) is parametrised by
$\mathcal{M}$
,
$\zeta$
,
$\mathcal{L}$
,
$L_0$
and
$\Delta _0$
(plus the parameter
$\varOmega \equiv \omega T_b$
relating to the source variation
$\mathcal{Q}$
). We assume henceforth that the flow domain is long (
$\mathcal{L}\gg 1$
). For fixed initial condition parameters
$\Delta _0$
,
$L_0$
, we will henceforth focus on behaviour across the
$(\zeta ,\mathcal{M})$
-parameter space, in particular, the quadrant
$\zeta \lesssim 1$
and
$\mathcal{M}\lesssim 1$
of relevance to weakly compressible gases.
The map in figure 3 will be helpful in navigating further reductions of this model. The boundaries of the map will be justified as we proceed; at present the reader is asked just to consider its overall organisation. As is typical of multi-parameter problems, distinguished limits arise in which different combinations of physical effects determine the dynamics. Figure 3 identifies codimension-0 regions (coloured areas) of parameter space, labelled
$\varPi _1$
–
$\varPi _4$
; these are separated by codimension-1 boundaries (lines we shall label
$\varPi _{ij}$
, bounding
$\varPi _i$
and
$\varPi _j$
), which intersect codimension-2 points (labelled
$\varPi _a$
–
$\varPi _c$
). Reduced forms of the model can be derived in each case, with the number of competing physical effects, and the number of relevant independent parameters, increasing with codimension. We will start in § 2.4 at the codimension-2 point
$\varPi _a$
(the limit
$\mathcal{M}\sim 1$
,
$\zeta \sim 1/\mathcal{L}\ll 1$
, parametrised by
$\mathcal{M}$
and
$\theta \equiv \zeta \mathcal{L}$
, as in (2.17)). Problem
$\varPi _a$
encompasses as a special case the codimension-1 problem
$\varPi _{12}$
, which we address in § 2.5 (the specialised limit
$\mathcal{M}\sim \theta \ll 1$
, parametrised by
$\mathcal{M}/\theta$
). Contained in problem
$\varPi _{12}$
are sublimits
$\varPi _1$
(
$\mathcal{M}/\theta \gg 1$
) and
$\varPi _2$
(
$\mathcal{M}/\theta \ll 1$
), represented respectively by the pink and blue regions of parameter space in figure 3. We will also highlight the codimension-2 problems
$\varPi _b$
(
$\mathcal{M}\sim \mathcal{L}^{-1/2}$
,
$\zeta \sim \mathcal{L}^{-3/2}$
) to assess the role of buoyancy (Appendix B) and
$\varPi _c$
(
$\mathcal{M}\sim \mathcal{L}^{-1}$
,
$\zeta \sim \mathcal{L}^{-2}$
) to assess flows at ultra-low-viscosity ratios (Appendix C). Within problem
$\varPi _c$
we identify additional codimension-1 problems
$\varPi _{23}$
,
$\varPi _{34}$
and
$\varPi _{14}$
that bound two remaining codimension-0 regions of parameter space (
$\varPi _3$
and
$\varPi _4$
, orange and green, respectively, in figure 3). Many of the subproblems will also involve consideration of distinct asymptotic regions in time (distinguishing early from late spreading) and space (near-source and far field). We will see that buoyancy effects (represented by the term
$F_x$
in (2.29b
)) are often subdominant to pressure gradients, and we will identify conditions under which spreading is driven directly by the volume flux
$\mathcal{Q}$
in (2.29c
) or indirectly by the evolving pressure field near the source. We now address the different physical balances within each limit; for later reference, these are summarised in table 3.
Table 3. A summary of asymptotic limits. Codimension-2 problems
$\varPi _a$
,
$\varPi _b$
,
$\varPi _c$
are formulated in distinguished asymptotic limits specified by two parameter balances (“
$\sim$
” denotes comparable magnitude), given using dimensionless parameters in figure 3 and here using dimensional parameters. These general problems capture multiple competing physical effects; only a subset of descriptors are given in the final column. Codimension-1 problems
$\varPi _{12}$
,
$\varPi _{23}$
,
$\varPi _{34}$
,
$\varPi _{14}$
are specified by a single parameter balance; descriptors are given in columns 2 and 3;
$P_b$
is shorthand for the near-source pressure of the gas bubble that connects to a thin gas film. Codimension-0 problems
$\varPi _1$
,
$\varPi _2$
,
$\varPi _3$
,
$\varPi _4$
occupy regions of
$(\zeta ,\mathcal{M})$
-parameter space (figure 3) and represent the most specific balances of physical effects; primary distinguishing features are summarised in columns 2–4.
2.4. Problem
$\varPi _a$
: weak compressibility
We first consider the distinguished limit (2.17) with
$\mathcal{M}=O(1)$
, large
$\mathcal{L}$
, small
$\zeta$
and
$\theta \equiv \zeta \mathcal{L}=O(1)$
. We distinguish the initial phase of spreading, when
$X_u=O(1)$
, from the later phase when
$X_u$
is comparable to
$\mathcal{L}\gg 1$
.
2.4.1. Early times
At early times, we pose the expansion
$P=P_0(t)/\zeta +P_1(x,t)+\ldots$
,
$F=F_0{(x,t)}+O(\zeta )$
, in which case (2.29) yields, at leading order,
\begin{align} P_0^{-1} \left [(1-F_0)P_0\right ]_t& = \left [ (1-F_0) P_{1x}\right ]_x \quad (0\lt x\lt X_u), \\[-9pt] \nonumber \end{align}
\begin{align} F_{0t}& = \mathcal{M} \left [F_0(P_{1x}+F_{0x})\right ]_x \quad (X_l\lt x\lt X_u), \\[-9pt] \nonumber \end{align}
\begin{align} - P_0 P_{1x}&=\mathcal{Q} \quad (x=0), \\[-9pt] \nonumber \end{align}
\begin{align} P_0+(P_{1x}+F_{0x}) \theta &=1 \quad (x=X_u), \\[-9pt] \nonumber \end{align}
\begin{align} X_{l,t}&=-\mathcal{M} (P_{1x}+F_{0x})\quad (x={X_l+}), \\[-9pt] \nonumber\end{align}
\begin{align} X_{u,t}&=- P_{1x}\quad (x={X_u-}). \\[9pt] \nonumber \end{align}
Functions in expansions will be assumed to remain
$O(1)$
in the relevant limit; here,
$P_0$
,
$P_1$
and
$F_0$
are assumed to be
$O(1)$
as
$\mathcal{L}\rightarrow \infty$
with
$\theta =O(1)$
and
$\mathcal{M}=O(1)$
. Initial conditions for
$P_0$
and
$P_1$
are obtained by expanding the initial condition for the full pressure field (2.30), yielding
\begin{align} P_0(0) = \frac {1}{2} - \frac {1}{{2 \Delta _0} }\left [ \theta - \left ( (1+4\theta ) \Delta _0^2 + \theta ^2 - 2 \theta \Delta _0\right )^{1/2} \right ]. \end{align}
The formulation (2.34) couples compressible effects, captured by the unsteady but spatially uniform pressure component
$P_0$
, with buoyancy and spreading effects captured by gradients of the smaller pressure component
$P_{1}$
.
To interpret (2.34), it is helpful to integrate (2.34a
) across regions I and II
$_g$
of figure 1, using boundary conditions (2.34c
,
f
) to give the mass balance for the gas phase
$[P_0 V(t)]_t=\mathcal{Q}$
; the gas volume
$V(t)$
is defined in (2.23). Integrating (2.34b
) across region II
$_w$
and using (2.34d
,
e
) gives a mass flux balance for the incompressible liquid phase
$V_t=\mathcal{M}(P_0-1)/\theta$
. Together, these yield
\begin{align} \frac {P_{0,t}}{P_0} V(t) = \frac {\mathcal{Q}}{P_0} - \frac {\mathcal{M} ( P_0 - 1)}{\theta }. \end{align}
This is an analogue of the differentiated form of Boyle’s law, which for a closed system with a fixed mass of gas would be
$(P_0 V)_t = 0$
. In contrast, (2.35) applies to an open system in which gas enters the channel; the first term on the right-hand side accounts for the volume flux of gas entering the channel, while the second term captures changes in gas volume due to liquid leaving the channel. Consequently, the gas compresses when the inflow exceeds the outflow. Equation (2.35) shows how
$P_0$
can equilibrate to a constant value
$P_\infty$
, say, as spreading initiates (while
$X_u$
remains small compared with
$\mathcal{L}$
), satisfying
\begin{align} P_\infty -\frac {\mathcal{Q} \mathcal{L}\zeta }{P_\infty \mathcal{M} } =1 \end{align}
(provided
$\mathcal{Q}$
is sufficiently slowly varying). Setting
$P_0=P_\infty$
in (2.34) recovers the incompressible flux constraint (2.26) with
$p_g$
replaced by
$P_1$
and
$\mathcal{Q}$
replaced by
$\mathcal{Q}/P_\infty$
. As anticipated in § 2.2, (2.35) shows how compressibility may lead to a rise in gas pressure prior to appreciable expansion of the gas domain.
Although derived for
$\mathcal{M}\sim 1$
and
$\zeta \mathcal{L} \sim 1$
, (2.35) extends to the limit
$\mathcal{M}\sim \zeta \mathcal{L} \ll 1$
, allowing
$P_0$
to equilibrate over a time scale
$t_r {\,=\,} O(1)$
while
$X_u$
remains small compared with
$\mathcal{L}$
. We note for later reference that when
$\mathcal{Q}=1$
for
$\mathcal{M}\ll \mathcal{L} \zeta \ll 1$
(region
$\varPi _2$
), (2.35) predicts that
$P$
rises to
$(\mathcal{L} /\mathcal{M} \zeta )^{1/2}$
over a time scale
$t_r\sim (\mathcal{L} \zeta /\mathcal{M})^{1/2}$
; in contrast, for
$ \mathcal{L} \zeta \ll \mathcal{M}\ll 1$
(region
$\varPi _1$
), (2.35) predicts that
$P$
rises to
$1/\zeta$
over a time scale
$t_r\sim \mathcal{L} \zeta /\mathcal{M}$
. Below, we compare
$t_r$
to the breakthrough time
$t_b$
at which the upper contact line first reaches the outlet.
2.4.2. Late times
Staying in the distinguished limit (2.17) of large
$\mathcal{L}$
and small
$\zeta$
with
$\mathcal{M}\sim 1$
, we now turn to spreading at later times when
$X_u$
becomes comparable to
$\mathcal{L}$
. Setting
$x=\mathcal{L} \bar {x}$
,
$t=\mathcal{L} \bar {t}$
,
$P=\mathcal{L} \bar {P}$
,
$F(x,t)=\bar {F}(\bar {x},\bar {t})$
, (2.29) becomes, at leading order with
$\theta \sim 1$
and
$\mathcal{M}\sim 1$
,
\begin{align} \left [(1-\bar {F)}\bar {P}\right ]_{\bar {t}} &= \left [ (1-\bar {F}) \bar {P} \bar {P}_{\bar {x}}\right ]_{\bar {x}} \quad (0\lt \bar {x}\lt \bar {X}_u), \\[-9pt] \nonumber \end{align}
\begin{align} \bar {F}_{\bar {t}}& = \mathcal{M} \left [ \bar {F}\bar {P}_{\bar {x}}\right ]_{\bar {x}}\quad (\bar {X}_l\lt \bar {x}\lt \bar {X}_u), \\[-9pt] \nonumber \end{align}
\begin{align} - \theta \bar {P} \bar {P}_{\bar {x}}&=\mathcal{Q}\quad (\bar {x}=0), \\[-9pt] \nonumber \end{align}
\begin{align} \theta \left [ \bar {P}+\bar {P}_{\bar {x}} (1-\bar {X}_u)\right ] &=1\quad (\bar {x}=\bar {X}_u), \\[-9pt] \nonumber \end{align}
\begin{align} \bar {X}_{l,\bar {t}}&=-\mathcal{M} \bar {P}_{\bar {x}}\quad (\bar {x}={\bar {X}_l+}), \\[-9pt] \nonumber\end{align}
\begin{align} \bar {X}_{u,\bar {t}}&=- \bar {P}_{\bar {x}}\quad (\bar {x}={\bar {X}_u-}). \\[9pt] \nonumber \end{align}
Over long spatial scales, buoyancy terms that would appear in (2.37b
,
d
,
e
), associated with tilting of the interface, are subdominant (of order
$1/\mathcal{L})$
, but may influence dynamics close to the lower contact line where the interface is likely to be steepest. Otherwise, (2.37) retains the viscous and compressible effects of the full model (2.29) throughout the region occupied by gas. Associated mass and volume constraints are
\begin{align} \left (\int _0^{\bar {X}_u} (1-\bar {F})\bar {P}\,\mathrm{d}\bar {x}\right )_{\bar {t}}=\frac {\mathcal{Q}}{\theta }, \quad \left (\int _0^{\bar {X}_u} (1-\bar {F})\,\mathrm{d}\bar {x}\right )_{\bar {t}}=\frac {\mathcal{M}}{1-\bar {X}_u}\left [\bar {P}\big \vert _{\bar {X}_u}-\frac {1}{\theta }\right ], \end{align}
implying that
$\bar {P}(\bar {X}_u,\bar {t})\rightarrow 1/\theta$
as
$\bar {X}_u\rightarrow 1$
. We next seek a reduction of this system when
$\mathcal{M}$
and
$\theta$
are both small, in which (2.37) develops a spatial inner/outer structure.
2.5. Problem
$\varPi _{12}$
: high gas mobility
We now refine the
$\varPi _a$
problem in the limit in which the viscosity ratio becomes small (
$\mathcal{L}^{-1}\ll \mathcal{M}\sim \zeta \mathcal{L}\ll 1$
), focusing on the late time structure (2.37) in which
$X_u$
moves over distances comparable to the channel length. Low
$\mathcal{M}$
can be expected to promote rapid spreading of a thin layer of gas above the liquid (Pegler et al. Reference Pegler, Huppert and Neufeld2014); likewise, (2.37e
) shows how motion of the lower contact line is suppressed. In the distinguished limit
$\theta \sim \mathcal{M} \ll 1$
, we therefore consider motion over short times
$\bar {t}=O(\mathcal{M})$
, during which a thin film extends from a shorter region in which the interface moves more slowly. Accordingly, we split the domain into an inner region over distances
$\bar {x}=O(\mathcal{M}^2)$
for which
$\bar {F}=O(1)$
(discussed in § 2.5.1), and an outer region over distances
$\bar {x}=O(1)$
for which
$1-\bar {F}=O(\mathcal{M})$
(discussed in § 2.5.2). The contact-line locations
$\bar {X}_l$
and
$\bar {X}_u$
lie in the inner and outer regions, respectively. Compressibility turns out to have a sustained effect in the outer thin film region.
2.5.1. Inner region
In the inner region, we set
$\bar {t}=\mathcal{M}\tilde {t}$
,
$\bar {x}=\mathcal{M}^2\tilde {x}$
and
$\bar {F}=\tilde {F}(\tilde {x},\tilde {t})+O(\mathcal{M})$
,
$\bar {P}=\mathcal{M}^{-1}\tilde {P}_0(\tilde {t})+\mathcal{M}^2\tilde {P}_1(\tilde {x},\tilde {t})+\ldots$
. Then (2.37) becomes
\begin{align} 0 &= \left [ (1-\tilde {F}) \tilde {P}_0 \tilde {P}_{1\tilde {x}}\right ]_{\tilde {x}}, \\[-9pt] \nonumber \end{align}
\begin{align} \tilde {F}_{\tilde {t}}& = \left [ \tilde {F}\tilde {P}_{1\tilde {x}}\right ]_{\tilde {x}}, \\[-9pt] \nonumber \end{align}
\begin{align} - \tilde {P}_0 \tilde {P}_{1\tilde {x}}&=\mathcal{Q}(\mathcal{M}/\theta ) \quad (\tilde {x}=0), \\[-9pt] \nonumber \end{align}
\begin{align} \tilde {X}_{l,\tilde {t}}&=-\tilde {P}_{1\tilde {x}} \quad (\tilde {x}={\tilde {X}_l+}). \\[9pt] \nonumber \end{align}
This is an unsteady incompressible viscous flow problem involving a so-far undetermined leading-order gas pressure
$\tilde {P}_0$
. The gas and liquid share the same pressure gradient
$\tilde {P}_{1\tilde {x}}$
in (2.38a
,
b
);
$\tilde {P}_0$
in (2.38a
,
c
), a proxy for gas density, converts volume fluxes in the gas into mass fluxes. It follows from (2.38a
,
c
) that
$(1-\tilde {F}) \tilde {P}_0 \tilde {P}_{1\tilde {x}}=-\mathcal{Q}(\mathcal{M}/\theta )$
throughout the inner region (the gas mass flux is uniform in this incompressible limit) and that
\begin{align} \tilde {F}_{\tilde {t}}+\frac {\mathcal{Q}}{\tilde {P}_0} \frac {\mathcal{M}}{\theta }\frac {\tilde {F}_{\tilde {x}}}{(1-\tilde {F})^2}=0. \end{align}
The interface height
$\tilde {F}$
can be integrated on characteristics. There will be some transient adjustment of
$\tilde {F}$
at early times via (2.34); for simplicity, we assume this is modest and impose the initial condition (2.8), giving
\begin{align} \tilde {F}= 1-\left [\frac {\mathcal{M}}{\theta }\frac {1}{\tilde {x}-\xi } \int _0^{\tilde {t}} \frac {\mathcal{Q}}{\tilde {P}_0}\,\mathrm{d}\tilde {t}'\right ]^{1/2}, \quad \tilde {F}=\frac {\xi -\tilde {X}_l(0)}{\tilde {X}_u(0)-\tilde {X}_l(0)} \end{align}
for
$\tilde {X}_l(0)\lt \xi \lt \tilde {X}_u(0)$
, which can be expressed as
\begin{align} (1-\tilde {F})^2 \left [\tilde {x}-\{\tilde {X}_l(0)+(\tilde {X}_u(0)-\tilde {X}_l(0))\tilde {F}\}\right ]= \frac {\mathcal{M}}{\theta }\int _0^{\tilde {t}} \frac {\mathcal{Q}}{\tilde {P}_0}\,\mathrm{d}\tilde {t} \end{align}
for
$0\lt \tilde {F}\lt 1$
. At the outer end of the inner region,
$\tilde {F}\rightarrow 1$
with
\begin{align} \tilde {F}\approx 1-\left [\frac {\mathcal{M}}{\theta }\frac {1}{\tilde {x}} \int _0^{\tilde {t}} \frac {\mathcal{Q}}{\tilde {P}_0}\,\mathrm{d}\tilde {t}'\right ]^{1/2}. \end{align}
This profile describes spreading of incompressible gas, from a source with a time-dependent volume flux proportional to
$\mathcal{Q}/\tilde {P}_0$
, into a thin gas film at the top of the aquifer. The inner solution depends on the unknown
$\tilde {P}_0(\tilde {t})$
, which must be determined by matching to the outer (thin-film) problem, but is independent of the precise initial condition.
2.5.2. Outer region
In the outer region, we set
$\bar {P}=\mathcal{M}^{-1}\hat {P}(\bar {x},\tilde {t})$
and
$\bar {F}=1-\mathcal{M}\hat {F}(\bar {x},\tilde {t})$
. Then (2.37) becomes, to leading order in
$\mathcal{M} \ll 1$
,
\begin{align} \left [\hat {F}\hat {P}\right ]_{\tilde {t}}& = \left [\hat {F} \hat {P} \hat {P}_{\bar {x}}\right ]_{\bar {x}} \quad (0\lt \bar {x}\lt \tilde {X}_u), \\[-9pt] \nonumber \end{align}
\begin{align} -\hat {F}_{\tilde {t}} &= \hat {P}_{\bar {x}\bar {x}} \quad (0\lt \bar {x}\lt \tilde {X}_u), \\[-9pt] \nonumber \end{align}
\begin{align} \hat {P}+\hat {P}_{\bar {x}} (1-\tilde {X}_u)& =\mathcal{M}/\theta \quad(\bar {x}=\tilde {X}_u), \\[-9pt] \nonumber \end{align}
\begin{align} \tilde {X}_{u,\tilde {t}}&=- \hat {P}_{\bar {x}}\quad (\bar {x}={\tilde {X}_u-}). \\[9pt] \nonumber \end{align}
The transport equation (2.39
$a$
) demonstrates that compressibility is significant along the thin gas film. Matching conditions at
$\bar {x}\rightarrow 0$
are provided by the outer limit (2.38h
) of the inner problem. Continuity of pressure, gas mass flux and of interface shape requires that
\begin{align} \hat {P}\rightarrow \tilde {P}_0,\qquad \hat {F}\hat {P}\hat {P}_{\bar {x}}\rightarrow -\frac {\mathcal{M}}{\theta }\mathcal{Q}, \qquad\hat {F}\approx \left [\frac {\mathcal{M}}{\theta }\frac {1}{\bar {x}} \int _0^{\tilde {t}} \frac {\mathcal{Q}}{\tilde {P}_0}\,\mathrm{d}\tilde {t}'\right ]^{1/2} \quad (\bar {x}\rightarrow 0). \end{align}
The associated mass and volume constraints are
\begin{align} \left (\int _0^{\tilde {X}_u} \hat {F}\hat {P}\,\mathrm{d}\bar {x}\right )_{\tilde {t}}=\frac {\mathcal{Q} \mathcal{M}}{\theta }, \qquad \left (\int _0^{\tilde {X}_u} \hat {F}\,\mathrm{d}\bar {x}\right )_{\tilde {t}}=\frac {1}{1-\tilde {X}_u}\left [\hat {P}\big \vert _{\tilde {X}_u}-\frac {\mathcal{M}}{\theta }\right ]. \end{align}
To match with the early time evolution (2.35), (2.36), we impose
$\hat {P}(0,0)=(\mathcal{M}/\theta )P_\infty$
in instances when
$P_0$
has had time to equilibrate.
2.6. Problems
$\varPi _1$
and
$\varPi _2$
The codimension-2 problem (2.37) (problem
$\varPi _a$
) is parametrised by
$\mathcal{M}$
and
$\theta$
; this problem retains compressibility but buoyancy becomes subdominant. The reduced codimension-1 problem
$\varPi _{12}$
emerging from
$\varPi _a$
, represented by (2.38), (2.39), is parametrised by the single parameter
$\mathcal{M}/\theta$
, which is formally of order unity in problem
$\varPi _{12}$
; here the low-viscosity ratio promotes formation of a thin film in which compressible effects remain significant, connected to an inner region where the flow is incompressible. We can recover two further simplifications, when
$\theta \ll \mathcal{M}\ll 1$
and
$\mathcal{M}\ll \theta \ll 1$
; these are the codimension-0 problems in regions
$\varPi _1$
and
$\varPi _2$
in figure 3. The problems differ in the outer region but share the same inner region structure, as described in § 2.5.1. We now address each outer region in turn.
(i) Region
$\varPi _1$
: for
$\theta \ll \mathcal{M}\ll 1$
, we expect compressibility effects to be suppressed. We set
$\hat {P}=(\mathcal{M}/\theta )+\hat {P}_1+\ldots$
, consistent with (2.36). The outer problem (2.39) becomes, to leading order in
$\theta /\mathcal{M}\ll 1$
,
\begin{align} \hat {F}_{\tilde {t}}& = \left [\hat {F} \hat {P}_{1\bar {x}}\right ]_{\bar {x}} \quad (0\lt \bar {x}\lt \tilde {X}_u), \\[-9pt] \nonumber \end{align}
\begin{align} -\hat{F}_{\tilde{t}} &= \hat {P}_{1\bar {x}\bar {x}}\quad (0\lt \bar {x}\lt \tilde {X}_u), \\[-9pt] \nonumber \end{align}
\begin{align} \hat {P}_1+\hat {P}_{1\bar {x}} (1-\tilde {X}_u)& =0\quad (\bar {x}=\tilde {X}_u), \\[-9pt] \nonumber \end{align}
\begin{align} \tilde {X}_{u,\tilde {t}}&=- \hat {P}_{1\bar {x}}\quad (\bar {x}={\tilde {X}_u-}), \\[-9pt] \nonumber \end{align}
\begin{align} \hat {F}\hat {P}_{1\bar {x}}&=-\mathcal{Q}, \qquad \hat {F}\approx \left [{\tilde {\mathcal{V}}}/{\bar {x}} \right ]^{1/2}\quad (\bar {x}\rightarrow 0). \\[9pt] \nonumber \end{align}
Here we have introduced
$\tilde {\mathcal{V}}(\tilde {t})\equiv \int _0^{\tilde {t}} \mathcal{Q}\,\mathrm{d}\tilde {t}'$
; (2.40e
) matches to (2.38h
). The mass and volume constraints are
\begin{align} \left (\int _0^{\tilde {X}_u} \hat {F}\,\mathrm{d}\bar {x}\right )_{\tilde {t}}=\mathcal{Q}, \qquad \left (\int _0^{\tilde {X}_u} \hat {F}\,\mathrm{d}\bar {x}\right )_{\tilde {t}}=\frac {1}{1-\tilde {X}_u}\hat {P}_1\big \vert _{\tilde {X}_u}, \end{align}
implying that
$\hat {P}_1\vert _{\tilde {X}_u}=\mathcal{Q}(1-\tilde {X}_u)$
. It follows from (2.40a
,
b
,
e
) that
$\hat {P}_{1\bar {x}}=-\mathcal{Q}/(1+\hat {F})$
and
\begin{align} \hat {F}_{\tilde {t}}+\frac {\mathcal{Q} \hat {F}_{\tilde {x}}}{(1+\hat {F})^2}=0. \end{align}
Solving using characteristics, assuming that at early times the interface is confined to a region
$\bar {x}=o(1)$
, gives
\begin{align} \hat {F}=\left [{\tilde {\mathcal{V}}}/{\bar {x}}\right ]^{1/2}-1 \quad (0\lt \bar {x}\lt \tilde {X}_u=\tilde {\mathcal{V}} ). \end{align}
In (2.40)–(2.42) we have recovered the incompressible spreading problem described by Pegler et al. (Reference Pegler, Huppert and Neufeld2014) and others, extended to account for time-dependent forcing, consistent with (2.28a ). The volume flux at the source has an instantaneous effect on the leading-order location of the upper contact line. The corresponding pressure field is
\begin{align} \bar {P}(\bar {x},\tilde {t})=\frac {1}{\theta }+\frac {\mathcal{Q}}{\mathcal{M}}\left [1-\frac {\tilde {\mathcal{V}}}{3}\left (1+2\left (\frac {\bar {x}}{\tilde {\mathcal{V}}}\right )^{3/2}\right ) \right ] \quad (0\lt \bar {x}\lt \tilde {\mathcal{V}} \leqslant 1). \end{align}
Here, with
$\theta \ll \mathcal{M}\ll 1$
, the large but uniform leading-order pressure
$1/\theta$
(obtained via rapid compression of the gas) effectively decouples from the pressure field driving spreading. However, (2.43) suggests that once
$\mathcal{M}$
and
$\theta$
are of comparable magnitude, viscous and compressible effects will together influence spreading.
(ii) Region
$\varPi _2$
: the other simpler case arising from (2.39) arises for
$\mathcal{M}\ll \theta$
, a limit in which we expect compressibility to be promoted. Here we set
$\hat {P}=(\mathcal{M}/\theta )^{1/2}\check {P}(\bar {x},\check {t})$
,
$\tilde {t}=(\theta /\mathcal{M})^{1/2}\check {t}$
,
$\hat {F}=\check {F}(\bar {x},\check {t})$
,
$\bar {X}_u=\check {X}_u(\check {t})$
, leading to
\begin{align} \left [\check {F}\check {P}\right ]_{\check {t}}& = \left [\check {F} \check {P} \check {P}_{\bar {x}}\right ]_{\bar {x}} \quad (0\lt \bar {x}\lt \check {X}_u), \\[-9pt] \nonumber \end{align}
\begin{align} -\check {F}_{\check {t}} &= \check {P}_{\bar {x}\bar {x}}\quad (0\lt \bar {x}\lt \check {X}_u), \\[-9pt] \nonumber \end{align}
\begin{align} \check {P}+\check {P}_{\bar {x}} (1-\check {X}_u)& =0\quad (\bar {x}=\check {X}_u), \\[-9pt] \nonumber \end{align}
\begin{align} \check {X}_{u,\check {t}}&=- \check {P}_{\bar {x}}\quad (\bar {x}={\check {X}_u-}), \\[-9pt] \nonumber\end{align}
\begin{align} \check {F}\check {P}\check {P}_{\bar {x}}\rightarrow -\mathcal{Q},& \qquad \check {F}\approx \left [\frac {1}{\bar {x}} \int _0^{\check {t}} \frac {\mathcal{Q}}{\check {P}(0,\check {t}^{\prime })}\,\mathrm{d}\check {t}'\right ]^{1/2}\quad (\bar {x}\rightarrow 0). \\[9pt] \nonumber \end{align}
The gas transport equation (2.44a
) balances viscous and compressible effects; it communicates with an incompressible inner region via (2.44e
), which matches to (2.38h
). The pressure field is smaller in magnitude and the spreading rate is slower than in problem
$\varPi _{1}$
. The system (2.44) satisfies the mass and volume constraints
\begin{align} \left (\int _0^{\check {X}_u}\check {F}\check {P}\,\mathrm{d}\check {x} \right )_{\check {t}}=\mathcal{Q},\qquad \check {V}_{\check {t}}\equiv \left (\int _0^{\check {X}_u}\check {F}\,\mathrm{d}\check {x} \right )_{\check {t}}=\frac {1}{1-\check {X}_u} \check {P}\big \vert _{\check {X}_u}. \end{align}
From (2.44c
,
d
,
f
) the volume constraint can be written as
$\check {V}_{\check {t}}= \check {X}_{u,\check {t}}$
. Integrating up to the breakthrough time
$\check {t}_b$
gives
$\check {V}(\check {t}_b) = \check {V}_0 + 1$
, where
$\check {V}_0=L_0/(\mathcal{LM})$
;
$\check {V}_0\ll 1$
in region
$\varPi _2$
. Thus, provided
$\mathcal{L}^{-1}\ll \mathcal{M}\ll 1$
, the dimensional delivered gas volume (per unit width)
\begin{align} V^*(t_b^*) \approx H L \mathcal{M} \end{align}
matches the volume delivered in the incompressible case
$\varPi _1$
, although the delivered mass
$\int _0^{t_b} \mathcal{Q}\,\mathrm{d}t$
will differ.
For
$\mathcal{Q}=1$
, the problem (2.44) is parameter free, yielding a solution that depends only on details of the initial condition. Numerical solutions of (2.44) are compared with solutions of the full problem (2.29) below.
2.7. Overview of asymptotic regimes
A summary of the asymptotic regions described so far is provided in figure 3. In the pink shaded region (
$\varPi _1$
), compressibility regulates rapid initial adjustment of the pressure via (2.35) but thereafter mass-flux-driven spreading is controlled by incompressible effects via (2.40). In the blue region (
$\varPi _2$
), viscous and compressible effects combine to determine the pressure scale and the spreading rate via (2.44). There is again a rapid adjustment of the pressure before spreading initiates. Further analysis at lower viscosity ratios reveals new balances. First, as explained in Appendix B, buoyancy effects allow the lower contact line to recede for
$\mathcal{L}\mathcal{M}^2\lesssim 1$
, without having an appreciable effect on spreading rates. Second, as explained in Appendix C, at even lower values of
$\mathcal{M}$
, the film of gas escaping from the bubble becomes so thin that it carries very low flux relative to that provided by the source. The bubble pressure rises at a rate determined by the source strength and compressive effects determine the breakthrough time. This is addressed in Appendix C by formulating the reduced problem at codimension-2 point
$\varPi _c$
, from which the codimension-1 problem
$\varPi _{12}$
is recovered, while revealing new codimension-0 regions
$\varPi _3$
and
$\varPi _4$
(shaded orange and green, respectively, in figure 3). In
$\varPi _3$
, compressible effects persist in the thin film and dynamic changes to the pressure of the main gas bubble influence spreading dynamics. In
$\varPi _4$
the thin-film flow is incompressible but compressive effects in the main gas bubble ensure that spreading is pressure-driven rather than directly driven by the imposed mass flux. The pressure rises to levels comparable to the incompressible problem in region
$\varPi _4$
(see (C12)), but to a level depending on
$\zeta$
when the gas is more compressible in region
$\varPi _3$
(see (C9)). These asymptotic predictions will now be evaluated against simulations of the full model (2.29).
3. Results
Figure 4. A simulation of (2.29) for
$\zeta = 10^{-4}$
,
$\mathcal{M} = 0.1$
and
$\mathcal{L} = 100$
, within the pink region of figure 3. Panels (a) and (b) show the pressure field
$P(x,t)$
and the interface height
$F(x,t)$
, respectively, at
$t = \{0, 0.1, 1, 4, 8, 10.4\}$
, following the light-to-dark colour gradient. Black dots in (a) show samples of the pressure at the upper contact line
$P(X_u,t)$
. Panel (c) shows the pressure evolution at the origin and at the upper contact line; the early time approximation given by (2.35) (dashed orange) captures the transient rise in pressure due to gas compressibility. Panel (d) shows contact-line locations;
$X_u$
approaches the late time solution (2.42) (dashed black line).
For a dimensionless domain length of
$\mathcal{L} = 100$
, we performed extensive simulations of (2.29) across a broad range of
$( \zeta , \mathcal{M})$
values. Throughout, we take
$L_0= \Delta _0 = 2$
in (2.22). In § 3.1 we present results for a steady injection rate (
$\mathcal{Q} = 1$
). To illustrate key flow behaviours, we show six representative simulations of (2.29): one (figure 4) from region
$\varPi _1$
in figure 3, one on the boundary
$\varPi _{12}$
(figure 5), three from region
$\varPi _2$
or its boundary with region
$\varPi _3$
(figures 6–8), and one from the ultra-low-viscosity regime near the
$\varPi _{3}/\varPi _4$
boundary (figure 9). In § 3.2 we compare results for increasing and decreasing injection rates.
3.1. Steady injection
To set the scene, we first illustrate spreading in the almost incompressible regime with low-viscosity ratio (
$\mathcal{M}=0.1$
,
$\zeta =10^{-4}$
; figure 4); as indicated in figure 3, this example lies in asymptotic region
$\varPi _1$
. As expected from previous studies of the incompressible problem with
$\mathcal{M}\lt 1$
(Pegler et al. Reference Pegler, Huppert and Neufeld2014), a finger of gas advances rapidly along the top of the channel (figure 4
b). Here
$X_u$
asymptotes promptly to
$t/\mathcal{M}$
(figure 4
d), advancing at a rate determined by the strength of the source, as predicted by (2.28b
); the lower contact line advances only marginally before the gas phase reaches the outlet. The term
$P(0,t)$
, which measures changes in gas pressure near the source relative to the hydrostatic scale
$\Delta \rho g H$
, is non-monotonic (figure 4
c): its early rise can be attributed to transient compressible effects, following (2.35); its later fall is due to the shortening length of the liquid-filled domain, across which there is the primary viscous pressure drop. The pressure field at any time is monotonic in
$x$
(figure 4
a) and is linear across
$X_u\lt x\lt \mathcal{L}$
. The pressure at the upper contact line rises rapidly before falling slowly (figure 4
a,c).
Figure 5 illustrates the impact of increasing the compressibility parameter to
$\zeta = 10^{-3}$
, a case lying on the line
$\mathcal{M} = \zeta \mathcal{L}$
in the parameter map shown in figure 3. The time taken for the pressure to rise is extended in comparison to figure 4; the early rise is still well captured by (2.35) (figure 5
c). Transient compression of the gas slows the rate of spreading relative to the incompressible case predicted by (2.28b
) (figure 5
$d$
). The spreading rate is further inhibited by the pressure reduction associated with increasing
$\zeta$
(compare figure 5
a to figure 4
a), which reduces the viscous pressure gradient along the liquid column. A more compressible gas (larger
$\zeta$
) has higher density, enabling it to sustain similar mass fluxes at lower pressures. The lower contact line in this example advances marginally (figure 5
b,d). This example corresponds to problem
$\varPi _b$
, discussed in Appendix B, so that buoyancy effects are expected to influence the motion of
$X_l$
. Otherwise, the pressure and interfacial profiles broadly resemble the incompressible example in figure 4.
A further increase in the compressibility parameter, to
$\zeta =1$
, brings more dramatic changes (figure 6). First, it significantly extends the proportion of spreading time in which
$P(0,t)$
rises as the gas is compressed (figure 6
c), almost until the breakthrough time. Spreading is substantially delayed relative to the incompressible case, giving time for buoyancy effects to drive the lower contact line backwards towards
$x=0$
(figure 6
b,d).
Figure 6. A simulation of (2.29) for
$\zeta = 1$
,
$\mathcal{M} = 0.1$
and
$\mathcal{L} = 100$
. Panels (a) and (b) show
$P(x,t)$
and
$F(x,t)$
at
$t = \{0,10,50,80,170,292.2\}$
. Solutions in (a–d) are plotted using the format described in figure 4. In (b, d), buoyancy forces drive the lower contact line backwards towards the origin and then up the wall at
$x = 0$
.
Reducing the viscosity ratio to
$\mathcal{M}=0.01$
with
$\zeta =10^{-3}$
also slows the rate at which the gas pressure initially rises (relative to the breakthrough time; compare figure 7
$c$
to figure 5
$c$
). In this example, with very low gas viscosity, spreading is sufficiently rapid for the initial interface profile to remain almost stationary as the gas film advances over the liquid (figure 7
$b$
), consistent with region-
$\varPi _3$
analysis in Appendix C. The thickness of the gas finger is reduced as a result of its low viscosity: the two layers have comparable fluxes driven by a shared pressure gradient, so that a thinner gas layer balances the smaller flux of the much more viscous liquid layer. Increasing
$\zeta$
, while retaining
$\mathcal{M}=0.01$
, again slows spreading substantially, giving buoyancy time to drive the lower contact line backwards (see figure 8
d for
$\zeta =1$
), as anticipated in Appendix B. Rather than equilibrating rapidly,
$P(0,t)$
rises continuously as the gas bubble elongates (figure 8
c), a signature of sustained compressive effects.
Figure 8. A simulation of (2.29) and the outer problem (2.44), for
$\zeta = 1$
,
$\mathcal{M} = 10^{-2}$
and
$\mathcal{L} = 100$
. The solid curves in (a) and (b) show the full solution of (2.29) plotted in the outer coordinates of problem
$\varPi _2$
,
$\check {P}(\bar {x},\check {t})$
and
$\check {F}(\bar {x},\check {t})$
, the dashed curves are numerical solutions of the outer equations (2.44). The solutions are shown at times
$\check {t} = \{0.8, 0.9, 1.1, 1.4, 1.7, 1.9\}$
. Solutions in (c) and (d) are from (2.29) and are plotted using the format described in figure 4.
A quantitative comparison with the predictions of problem
$\varPi _2$
is illustrated in figure 8, where solutions of (2.29) are compared with solutions of the reduced model (2.44). The close agreement illustrates that the reduced model captures the dominant processes. It is helpful at this point to review the elements of the reduced model. Equation (2.44a
) describes the evolution of the mass of gas, integrated across the thin film and represented by the density
$\check {F}\check {P}$
, advected by the pressure gradient
$\check {P}_{\bar {x}}$
. The gas layer and the thicker liquid beneath share the same pressure gradient. Equation (2.44b
) describes the evolution of the mass of liquid below the gas film. Across this region
$\check {P}_{\bar {x}\bar {x}}\lt 0$
, indicating stretching of the liquid column that allows the gas thickness to increase locally. Equation (2.44c
) captures the viscous pressure drop of the liquid column ahead of the gas bubble; as the column shortens, the pressure at the bubble tip falls. Equation (2.44d
) is the kinematic relationship between contact-line speed and pressure gradient: as the liquid column shortens, the pressure gradient across it rises and the contact line accelerates. Finally, (2.44e
) couples mass input from the source to compression of gas in the bulk of the bubble. This is awkward to implement numerically: we implemented a local expansion that fails to capture the precise shapes of
$\check {F}$
curves for smaller
$\bar {x}$
in figure 8.
Further reducing the viscosity ratio to
$\mathcal{M} = 10^{-3}$
with
$\zeta = 10^{-4}$
(figure 9) reveals the characteristic features of the ultra-low-viscosity regime described in Appendix C. In this regime, most of the interface remains stationary throughout the evolution, except for the thin film along the upper boundary; this is so narrow (figure 9
b) that the source supplies a negligible flux to this region. Consequently, the pressure rises within the main gas bubble at a rate that is approximately linear in time when
$\mathcal{Q}=1$
(figure 9
c); the slight deviation of
$P(0,t)$
from linearity reflects the slight expansion of the gas domain as the thin film elongates. An analytical solution of the reduced-order model (C12) for
$\mathcal{Q}=1$
implies that
\begin{align} {X_u}(t)=3\mathcal{L}\left [1-\left (1-\frac {t^2}{3L_0 \zeta \mathcal{L}^2}\right )^{1/2}\right ]\quad \left (0\lt t\lt \mathcal{L}\sqrt {5L_0\zeta /3}\right ). \end{align}
In this limit, spreading is driven primarily by the buildup of pressure in the gas bubble rather than directly from the mass flux delivered by the source. As a result, the position of the upper contact line
$X_u$
initially grows quadratically in time. The approximation (3.1) agrees well with the full numerical simulation of (2.29) across the entire evolution (figure 9
d).
Figure 9. A simulation of (2.29) in the ultra-low-viscosity regime, with
$\zeta = 10^{-4}$
,
$\mathcal{M} = 10^{-3}$
and
$\mathcal{L} = 100$
. Solutions are plotted using the same format as figure 4. (a,b) Evolution of the pressure
$P(x,t)$
and interface height
$F(x,t)$
at times
$t = \{0, 0.1, 0.5, 1, 1.5, 1.8\}$
. The inset in (b) shows
$F$
where it is close to unity. The dashed pink curve in (c) is the approximation of the linear pressure rise within the main gas bubble, given by (C12
$e$
). The dashed red line in (d) is the approximation given by (3.1).
Figure 10. Contour plots showing (a) breakthrough times
$t_b$
and (b) density
$\zeta P$
at the origin at the breakthrough time, in
$(\zeta , \mathcal{M} )$
-parameter space. Dots capture data from 2500 simulations conducted with
$\mathcal{L} = 100$
. Black contour lines, derived from interpolated data, are evenly spaced on a logarithmic scale. The magenta line illustrates the approximate location of the asymptotic boundary
$\varPi _{12}$
separating region
$\varPi _1$
to the left from region
$\varPi _2$
to the right (regions are illustrated in figure 3).The blue triangle in (a) indicates a slope of
$-1$
, illustrating the predicted scaling
$t_b \sim (\mathcal{L}^3 \mathcal{M} \zeta )^{1/2}$
. The asymptotic boundary
$\varPi _{23}$
lies along
$\mathcal{M}\sim \mathcal{L}^{-1}$
at the base of the map. Here
$t_b$
and
$\zeta P$
are predicted to be independent of
$\mathcal{M}$
in region
$\varPi _3$
, beneath this boundary (see Appendix C).
Computations across the parameter space show how the breakthrough time
$t_b$
increases with both
$\zeta$
and
$\mathcal{M}$
(figure 10
a). To the left of the line
$\zeta \mathcal{L} = \mathcal{M}$
, the flow approaches the incompressible limit (problem
$\varPi _1$
), making the breakthrough time independent of compressibility. In this regime,
$t_b\approx \mathcal{L} \mathcal{M}$
, which is reflected in the horizontal levelling off of the contours. To the right, problem
$\varPi _2$
predicts that compressible and viscous effects determine the breakthrough time, with
$t_b \sim (\mathcal{L} / \mathcal{M} \zeta )^{1/2}$
. This scaling implies that the contours should have a slope of
$-1$
, as demonstrated in the top right corner of the figure. Closer to the lower boundary of region
$\varPi _2$
(i.e. the lower boundary of figure 10
a), the contours become distorted due to the influence of very small viscosity effects. In the limit of vanishing gas viscosity (region
$\varPi _3$
in figure 3), the breakthrough time scale
$t_b \sim (\zeta \mathcal{L}^2)^{1/2}$
becomes independent of
$\mathcal{M}$
(Appendix C), manifesting as the contours becoming vertical.
Figure 10(b) shows how
$\zeta P$
, the gas density measured relative to
$\rho _{g0}$
, increases with
$\zeta$
and decreases with
$\mathcal{M}$
across region
$\varPi _2$
. The dominant density in the incompressible region (
$\varPi _1$
, to the left of the figure) is
$O(1)$
, with a correction of
$O(\zeta \mathcal{L}/\mathcal{M})$
(arising from the pressure
$\hat {P}_1$
in (2.40)). Accordingly, the density contours share the slope of the line
$\mathcal{L}\zeta =\mathcal{M}$
near this region. Across region
$\varPi _2$
,
$\zeta P$
is
$O(\zeta \mathcal{L}/\mathcal{M})^{1/2}$
(see (2.44)) and, therefore, the contours maintain a unit slope. However, near the lower boundary of the figure, the contours become vertical, because
$\zeta P=O(\mathcal{L}\zeta ^{1/2})$
in the underlying region
$\varPi _3$
(figure 3), which is independent of
$\mathcal{M}$
.
Figure 11. The solutions of (2.29). Panels (a) and (b) are for
$\mathcal{M} = 0.1$
, while panels (c), (d), (e), (f) corresponds to
$\mathcal{M} = 0.01$
. Panels (a) and (c) show the pressure at the origin and panels (b) and (d) show the upper contact line, all replotted in the outer coordinates of (2.44), given by
$P = (\mathcal{L}/\mathcal{M} \zeta )^{1/2} \check {P}$
,
$X_u = \mathcal{L} \check {X}_u$
and
$t = (\mathcal{L}^3 \mathcal{M} \zeta )^{1/2} \check {t}$
. The solid curves represent the full numerical solution of (2.29), while the dashed orange curves correspond to the outer approximation obtained by numerically solving (2.44). For both values of
$\mathcal{M}$
, the outer problem was initialised using the full solution for
$\zeta = 1$
at
$t = 80$
, by which time the late time inner–outer structure has fully developed. Panel (e) shows the speed of the upper contact line for
$\mathcal{M} = 0.01$
, which is compared with the speed of incompressible propagation (dashed black curve). In (f) we again replotted the curves presented in (e) in the outer coordinates, the dashed cyan curve
$\check {t}/2$
is plotted for comparison.
To investigate further the underlying scaling relationships, we replot simulation results using variables appropriate to region
$\varPi _2$
(the blue region in figure 3). Despite parameter variation over three orders of magnitude, the gas pressure at the source, the overall pressure distribution and the location of the upper contact line show only minor relative variation when plotted relative to
$\check {t}=t/(\mathcal{L}^3\mathcal{M}\zeta )^{1/2}$
(figure 11
a–d); the outlying case with
$\zeta =10^{-3}$
lies on the boundary with region
$\varPi _1$
. The speed of the upper contact line is shown in regular coordinates in figure 11(e) for multiple values of
$\zeta$
, with
$\mathcal{M} = 0.01$
, at the border between regions
$\varPi _2$
and
$\varPi _3$
. Initially, the contact line undergoes a brief deceleration, followed by a phase of approximately constant acceleration. As it approaches the outlet, a sharp increase in speed is observed, driven by the rapid decrease in viscous resistance from the draining liquid. For reference, the constant incompressible contact-line speed, given by
$1/\mathcal{M}$
, is indicated in figure 11(e). Across all values of
$\zeta$
considered, the contact-line speed remains below this incompressible limit throughout. When the same data are replotted in the outer coordinate of problem
$\varPi _2$
, the curves again collapse (figure 11
f) with
$\check {X}_{u,\check {t}}\approx ({1}/{2})\check {t}$
. The quadratic dependence of contact-line location on time for parameters at the interface of regions
$\varPi _2$
and
$\varPi _3$
(figure 11
d,e,f), which (following (3.1)) we attribute to dynamic pressure changes at the source, contrasts with the flux-driven linear dependence in the incompressible limit (2.28b
).
Figure 12. Numerical simulations of (2.29) for time-dependent injection rates, with
$\mathcal{M} = 10^{-2}$
,
$\zeta = 0.1$
,
$\mathcal{L} = 100$
and
$\varOmega = 10^{-2}$
. The corresponding steady injection case is shown for comparison. (a) Evolution of the pressure at the origin. (b) Evolution of the pressure at the upper contact line.
3.2. Time-dependent injection
Finally, we briefly consider two cases of time-dependent injection to examine how variations in the injection rate affect the total mass stored in the channel at the breakthrough time. Specifically, we investigate an increasing rate,
$\mathcal{Q}(\varOmega t) = 1 + \varOmega t$
, and a decreasing rate,
$\mathcal{Q}(\varOmega t) = 1 - \varOmega t$
, as defined in (2.24). Since the breakthrough time is not known a priori, we take
$\varOmega \ll 1$
to prevent mass withdrawal from the system.
Figure 12 compares the pressure evolution at the source and at the upper contact line for constant and time-dependent injection rates, with
$\mathcal{L} = 100$
,
$\mathcal{M} = 10^{-2}$
,
$\zeta = 10^{-1}$
and
$\varOmega = 10^{-2}$
(on the region
$\varPi _2/\varPi _3$
boundary). The initial pressure field (2.30) may be substituted into (2.33) to evaluate the initial mass of gas in the channel, giving
$M_{g0} = 3$
for
$\zeta = 10^{-1}$
. For steady injection, the breakthrough time is
$t_b \approx 60$
, corresponding to a final mass of
$M_g(t_b) \approx M_{g0} + 60$
. When the injection rate increases over time, the pressure rises more rapidly and the elevated pressure scale and steeper pressure gradient at
$X_u$
lead to faster spreading, resulting in an earlier breakthrough at
$t_b \approx 55$
with
$M_g(t_b) \approx M_{g0} + 71$
. Conversely, a decreasing injection rate shortens the pressure rise time, allowing the pressure to saturate near the source by the time that breakthrough occurs (figure 12
a). Quicker pressure saturation is consistent with the early time problem in (2.35), where a decreasing injection rate leads to a decline in
$P_{0,t}$
. Spreading is also slowed, so that the breakthrough time is delayed to
$t_b \approx 68$
, at which time the mass in the channel (
$M_g(t_b) \approx M_{g0} + 45$
) is less than the cases with steady or increasing injection rate. These findings suggest that an increasing injection rate enhances storage efficiency, since the total gas volume remains the same at breakthrough across all cases according to (2.45).
4. Discussion
4.1. Dimensional scales
We have used a simplified model of gas injection into a confined liquid-filled porous medium in order to investigate the role of gas compressibility. In a long-wave limit the spreading of the gas is modelled by coupled evolution equations (2.29) for the gas pressure and interface height. Focusing on the regime
$\mathcal{M}\ll 1$
(high gas mobility, leading to spreading being confined to a thin layer of gas that advances above the brine),
$\mathcal{L}\gg 1$
(a long domain length) and
$\zeta \ll 1$
(weak gas compressibility), we identified three dominant regions of parameter space (figure 3) in which spreading is regulated by different dominant balances: in region
$\varPi _1$
, compressible effects are transient and spreading is essentially incompressible, being regulated by the source strength; in
$\varPi _2$
, compressible effects regulate dynamics in the thin gas film and, hence, the overall spreading rate; in
$\varPi _3$
and
$\varPi _4$
, injection drives a steady increase of the pressure in the main gas bubble and spreading is confined to an ultra-thin film of gas. Figures 4–11 provide both qualitative and quantitative validation of the behaviour across all regimes, at parameter values indicated in figure 3. A table summarising these codimension-0 regions of
$(\zeta ,\mathcal{M})$
-parameter space, their codimension-1 boundaries and the codimension-2 points where the boundaries intersect is provided in table 3.
It is helpful to revisit the predicted time and pressure scales in these regions, expressing them using dimensional quantities. In the incompressible region
$\varPi _1$
, the scales for pressure, breakthrough time and rise time are
\begin{align} p_g^* - p_{g0} \sim \frac {q \mu _w L}{k_0 \rho _{g0} H} , \quad t_b^* \sim \frac { \rho _{g0} H L \mu _g}{q \mu _w}, \quad t_r^* \sim \frac {\mu _w L \rho _w g H^2}{q \mu _g c^2}, \end{align}
respectively. Here, the pressure scale (above the reference pressure
$p_{g0}$
at the downstream end of the aquifer) is primarily set by the viscous pressure drop along the liquid-filled column. The breakthrough time (4.1b
) is set by the source-driven transit time,
$\rho _{g0} H L / q$
, reduced by the mobility ratio,
$\mu _g / \mu _w$
, to account for the formation of a thin gas film over the liquid, which facilitates faster spreading. The much shorter scale for the pressure rise time (4.1c
) is governed by a balance between viscous dissipation, gas compressibility and buoyancy. The full-dimensional gas pressure field for the incompressible problem
$\varPi _1$
is, from (2.43),
\begin{align} p_g^* = p_{g0} + \frac {q \mu _w L \mathcal{Q} }{k_0 \rho _{g0} H} - \frac {q^2 \mu _w^2 \mathcal{V}^* \mathcal{Q}}{3\rho _{g0}^2 \mu _g k_0 H^2} \left [1 + 2 \left (\frac { \mu _g \rho _{g0} H}{\mu _w q} \frac {x^*}{\mathcal{V}^*} \right )^{3/2} \right ], \quad \mathcal{V}^* = \int _0^{t^*} \mathcal{Q} \, \mathrm{d} t^*, \end{align}
for
$x^*$
between the source and the upper contact-line location
$X_u^*(t^*)= q \mu _w \mathcal{V}^*/(\rho _{g0} H \mu _g)$
. The channel length
$L$
dictates the large uniform pressure component needed to drive the motion; this builds in the gas through rapid compression over the time scale
$t_r^*$
. However, this pressure component is decoupled from the spatially non-uniform spreading pressure field, which remains independent of the domain length.
The corresponding dimensional scales in region
$\varPi _2$
of figure 3, when compressible effects become significant, are
\begin{align} p_g^*\! -\! p_{g0}\! +\! \rho _{g0} c^2\! \sim \left [\frac {q \mu _w c^2L}{H k_0 } \right ]^{1/2}\!, \, t_b^* \sim \left [\frac {H \mu _g^2 L^3}{q k_0 \mu _w c^2} \right ]^{1/2}\!, \, t_r^* \sim \left [\frac {\rho _w^2 g^2 H^5 k_0 \rho _{g0}^2 \mu _w L }{\mu _g^2 q^3 c^2}\right ]^{1/2}. \end{align}
The pressure scale (4.3a
) is given by the geometric mean of the viscous pressure drop along the liquid-filled channel and the compressible pressure scale
$\rho _{g0}c^2$
. The breakthrough time scale (4.3b
) emerges from a leading-order balance between gas compressibility and viscous dissipation in both the liquid and gas. The shorter pressure rise time (4.3c
) is determined by all physical effects in the model. The scales (4.3a,b,c
) show how, for a more compressible gas (with smaller
$c^2$
), the pressure required to displace the liquid falls but the time taken to deliver the gas rises. As demonstrated in (2.33), the breakthrough time determines the mass of gas delivered when the source injection is steady, via
$q t_b^*$
.
In region
$\varPi _3$
(the green region in the parameter map, figure 3; Appendix C), the dimensional scales simplify to
\begin{align} p_g^* - p_{g0} + \rho _{g0} c^2\sim \frac {q \mu _g L}{k_0 \rho _{g0}H} \left ( \frac {\rho _{g0} c^2}{\rho _w g H} \right )^{1/2}, \qquad t_b^* \sim t_r^* \sim \frac {\rho _{g0} H L }{q} \left ( \frac {\rho _w g H }{\rho _{g0} c^2} \right )^{1/2}. \end{align}
The pressure scale is independent of the liquid viscosity but now depends on buoyancy: this is because the very thin gas layer displaces the liquid by slightly depressing it, overcoming hydrostatic pressure rather than displacing a large column of liquid along the channel against viscous resistance. The time scales for breakthrough and pressure rise are governed by the source-driven transit time, modulated by compressibility and buoyancy effects. Here we have taken the initial bubble length
$L_0$
to be order unity. If instead we let
$L_0 \rightarrow 0$
, the dimensionless scales become
$P \sim (\mathcal{L}^2/ \zeta L_0)^{1/2}$
,
$t_b \sim (\zeta \mathcal{L}^2 L_0)^{1/2}$
, showing how the pressure scale increases for a smaller initial bubble and the breakthrough time decreases; the scale for the breakthrough time is consistent with (3.1). The final region,
$\varPi _4$
in figure 3 (Appendix C), shares the pressure scale of region
$\varPi _1$
(4.1a
) and the time scale of region
$\varPi _3$
(4.4b
).
At the breakthrough time, the gas is predominantly stored in the film (the outer region) for
$\mathcal{L}^{-1}\ll \mathcal{M}\ll 1$
and in the bubble (the inner region) for
$\mathcal{ML}\ll 1$
. Scaled relative to
$\rho _{g0}HL$
, the mass of gas is
$O(\mathcal{M})$
in region
$\varPi _1$
, rising to
$(\mathcal{L} \mathcal{M} \zeta )^{1/2}$
in region
$\varPi _2$
. In the incompressible case (
$\varPi _1$
), the viscosity ratio determines the film thickness, however, the mass of gas is independent of other flow-related parameters. The mass of stored gas when
$\mathcal{ML}\ll 1$
is
$O(L_0\zeta ^{1/2})$
and the dimensional mass is of order
$LL_{g0}q\mu _g /[k_0 (\rho _{g0} \rho _w g H)^{1/2} c]$
. Here the pressure (and, therefore, the gas density) is elevated by the flow of gas in the film.
Having reviewed the primary regions, we comment briefly on the boundaries between them. These are specified by the thresholds
$\varPi _{12}$
,
$\varPi _{23}$
,
$\varPi _{34}$
,
$\varPi _{14}$
summarised in table 3. A primary condition defining when compressibility strongly influences spreading is the requirement
$\mathcal{M}\ll \mathcal{L}\zeta$
when
$\mathcal{L}^{-2}\ll \zeta \ll \mathcal{L}^{-1}$
, a threshold defined by the boundary between regions
$\varPi _1$
and
$\varPi _2$
in figure 3. Expressed in terms of input parameters, this requires that the pressure generated by viscous pressure drop along the channel is sufficiently large to cause appreciable density changes in the gas, i.e.
\begin{align} \frac {q \mu _w L}{k_0\rho _{g0}H} \gg \rho _{g0} c^2. \end{align}
Thus, compressibility effects are promoted by rapid injection into long, narrow, low-permeability channels. An alternative pathway arises at much lower viscosity ratios, involving a transition from flux-driven to pressure-driven spreading on crossing boundary
$\varPi _{14}$
, followed by a transition to flux-driven spreading on crossing boundary
$\varPi _{34}$
with
\begin{align} \left (\frac {q \mu _g L}{k_0 \rho _{g0} H}\right )^2 \gg (\rho _{g0} c^2) (\rho _w g H). \end{align}
Here viscous losses in the gas must exceed a threshold defined by compressive effects and buoyancy.
4.2. Practical implications
The results for time-dependent injection rates (figure 12) suggest that variations in injection rate may provide an effective means of controlling key aspects of flow and storage. Depending on the primary concerns at a given storage site, either an increasing or decreasing injection rate may offer advantages. An increasing injection rate can lead to elevated pressure buildup and more rapid spreading, potentially compromising reservoir integrity; however, if the reservoir is sufficiently robust, this strategy enables greater mass storage by the time of breakthrough. Conversely, a decreasing injection rate slows the spreading and results in less total mass delivered at breakthrough. Nevertheless, the associated reduction in pressure scale may reduce mechanical stress on the reservoir rock, making this approach preferable in settings where the risk of fracturing is a primary concern.
To illustrate the practical implications of compressibility, we report in table 4 the ranges of the dimensionless parameters in (2.29) for hydrogen storage, together with values calculated from the data provided in table 4 of Zheng & Neufeld (Reference Zheng and Neufeld2019) for two CO
$_2$
storage sites: the Sleipner aquifer in the North Sea (Boait et al. Reference Boait, White, Bickle, Chadwick, Neufeld and Huppert2012) and the In Salah aquifer in Algeria (Vasco et al. Reference Vasco, Rucci, Ferretti, Novali, Bissell, Ringrose, Mathieson and Wright2010). Estimates of the speed of sound for CO
$_2$
were obtained from Lemmon et al. (Reference Lemmon, Bell, Huber, McLinden, Linstrom and Mallard2023). For hydrogen, the ranges were computed from the typical dimensional parameters presented in table 1, yielding dimensionless parameters spanning several orders of magnitude. This reflects the diversity of potential storage sites in terms of geometry, thermodynamic conditions and pore-scale permeability. The lower end of
$\mathcal{L}$
, corresponding to extreme cases of very slow injection into a short, highly permeable and deep channel, lies outside the validity of the present model. Therefore, restricting attention to
$\mathcal{L} \gtrsim 100$
places hydrogen in the compressible regime
$\varPi _2$
(based on the threshold
$\varPi _{12}$
;
$\mathcal{M} \ll \theta$
) and potentially close to the boundary
$\varPi _{23}$
where ultra-low-viscosity effects may influence spreading. The parameter ranges in table 4 for Sleipner and In Salah data places the CO
$_2$
spreading in regions
$\varPi _1$
and
$\varPi _2$
, respectively. This suggests that spreading at Sleipner is effectively incompressible, whereas at In Salah compressibility may be significant, thereby slowing the rate of plume evolution.
Table 4. Estimated dimensionless parameter ranges for hydrogen and CO
$_2$
. The hydrogen ranges (column 2) are derived from table 1, reflecting the wide variability in geometry, thermodynamic conditions and permeability of potential underground storage sites. Column 3 lists estimates for CO
$_2$
storage at the Sleipner aquifer in the North Sea (Boait et al. Reference Boait, White, Bickle, Chadwick, Neufeld and Huppert2012), based on the data in table 4 of Zheng & Neufeld (Reference Zheng and Neufeld2019). Column 4 gives corresponding estimates for CO
$_2$
storage at the In Salah aquifer in Algeria, also from Zheng & Neufeld (Reference Zheng and Neufeld2019).
4.3. Model limitations and final remarks
The present model rests on numerous assumptions, which may limit its applicability to practical applications. In particular, inertial effects may become significant in coarse-grained porous media, especially for gases, which typically exhibit low viscosity. An extension of Darcy’s law (the Darcy–Forchheimer equation) accounts for the additional inertial drag in high-Reynolds-number flows (Farahani et al. Reference Farahani, Chiapponi, Longo and Federico2024). At the opposite end of the pore-size spectrum, when pore diameters approach the mean free path of gas molecules, the continuum assumption can break down. At these small scales, molecular interactions with solid boundaries become more frequent than intermolecular collisions, giving rise to non-continuum transport phenomena. This includes pressure-dependent permeability through the emergence of slip flow (Wu, Pruess & Persoff Reference Wu, Pruess and Persoff1998) and Knudsen diffusion, a mechanism driven by molecule–wall interactions that dominate over bulk-phase collisions (Darabi et al. Reference Darabi, Ettehad, Javadpour and Sepehrnoori2012). A simple model to account for slip flow, originally proposed by Klinkenberg (Reference Klinkenberg1941), introduces a nonlinear correction to the permeability that is inversely proportional to pressure. In the context of this study, lower pressures were observed for more compressible gases, which may enhance slip flow relative to weakly compressible gases.
We presented a model for a regime where capillary forces are subdominant to gravitational forces, in which case the interface will be sharp and the pressure across the interface continuous. This modelling assumption is consistent with recent experiments by Mortimer, Mingotti & Woods (Reference Mortimer, Mingotti and Woods2024), which showed that the interface remained relatively sharp under slow air injection into a water-saturated tank of glass beads. However, when capillary forces become significant, the sharp-interface approximation may no longer be valid, and it may become necessary to incorporate a macroscale capillary pressure jump into the model. Increased capillary forces also promote residual trapping of isolated gas bubbles through capillary snap-off (Lysyy, Ersland & Fernø Reference Lysyy, Ersland and Fernø2022). The extent of this trapping is governed by the system’s wettability and the heterogeneity of the pore structure. In CO
$_2$
storage, residual trapping is beneficial since storage is intended to be permanent; in contrast, this may complicate recovery in hydrogen storage. However, once isolated, trapped bubbles may undergo Ostwald ripening, whereby spatial variations in capillary pressure drive diffusive mass transfer from smaller to larger bubbles through the surrounding liquid (Zhang et al. Reference Zhang, Bijeljic, Gao, Goodarzi, Foroughi and Blunt2023). This redistribution could promote the formation of larger, more mobile gas clusters, potentially enhancing recoverability.
Given the potential importance of capillary effects, it is useful to identify conditions under which the sharp-interface approximation is valid. Zheng & Neufeld (Reference Zheng and Neufeld2019) extended the incompressible injection models of Pegler et al. (Reference Pegler, Huppert and Neufeld2014) and Zheng et al. (Reference Zheng, Guo, Christov, Celia and Stone2015) by incorporating partial saturation and a saturation-dependent capillary pressure. Their formulation employed a Brooks–Corey capillary–pressure saturation relation together with the multiphase extension of Darcy’s law (Bear Reference Bear1972), which captures reductions in permeability due to partial saturation. The resulting dimensionless formulation depends on three key parameters: the viscosity ratio scaled by the endpoint relative permeability of the non-wetting phase,
$k_{rn0}/\mathcal{M}$
; the Bond number,
$Bo = \Delta \rho g H / p_e$
, where
$p_e$
is the capillary entry pressure; and an empirical pore-size distribution parameter,
$\varLambda$
. In this framework, a sharp interface is recovered in two limits:
$Bo \to \infty$
(see the estimate in § 2), where buoyancy forces dominate over all capillary entry pressures, and
$\varLambda \to \infty$
, corresponding to monodisperse pores with constant identical entry pressures, so that displacement occurs as a single advancing front rather than a diffuse transition zone.
In addition to these capillary-driven interfacial phenomena, fluid displacement may also be destabilised dynamically through the onset of viscous fingering, which arises when a low-viscosity fluid displaces a more viscous one. A detailed discussion of how this instability impacts underground hydrogen storage can be found in Paterson (Reference Paterson1983). It is noteworthy, however, that certain features of liquid displacement by a compressible gas can act to delay the onset and mitigate the severity of such instabilities. In particular, the large density contrast between hydrogen and brine has a stabilising effect on the interface, as demonstrated by Saffman & Taylor (Reference Saffman and Taylor1958). More recently, Cuttle et al. (Reference Cuttle, Morrow and MacMinn2023) showed how higher gas compressibility can have an effect comparable to that of an elevated capillary number in delaying the onset of fingering.
While interfacial instabilities such as viscous fingering shape the large-scale displacement pattern, the rheological properties of the displacing fluid itself can also influence the flow dynamics. In addition to the shear viscosity of the fluids, the compressible Navier–Stokes equations reveal an additional viscous contribution associated with volumetric deformation, quantified by the bulk viscosity,
$\mu _b = \lambda + ({2}/{3}) \mu$
, where
$\lambda$
is the second (or volume) viscosity (Li, Hu & Jiang Reference Li, Hu and Jiang2017) that is not accounted for in Darcy’s law. As a fundamentally pore-scale effect, the role of bulk viscosity, and its potential to be systematically upscaled to reservoir-scale flow, remains unclear. However, when energy dissipation due to volumetric deformation is negligible, the Stokes hypothesis (
$\mu _b = 0$
) is often invoked (Graves & Argrow Reference Graves and Argrow1999). This approximation has been shown to hold exactly for monatomic gases such as helium, and approximately for weakly compressible polyatomic gases like hydrogen (Graves & Argrow Reference Graves and Argrow1999). However, for highly compressible gases such as CO
$_2$
, the bulk viscosity can exceed the shear viscosity by several orders of magnitude (Pan & Johnsen Reference Pan and Johnsen2017) and may therefore have a significant influence on flow.
Gas compression leads to localised temperature increases, resulting in both spatial and temporal temperature fluctuations within the injected gas. Capturing these effects would require a more sophisticated equation of state, coupled with conservation of energy. However, in porous media, the key question is whether heat exchange with the solid matrix is sufficient to suppress such variations. In this study we considered isothermal compression. When analysing temperature and pressure variations during cyclical injection and withdrawal in underground compressed-air storage, Kushnir et al. (Reference Kushnir, Ullmann and Dayan2012b
) identified the dominant factor distinguishing adiabatic from isothermal behaviour in a porous medium as the rate of heat transfer between the gas and the solid matrix. This rate can be quantified by the thermal effusivity,
$e = \sqrt {\kappa \rho c_p}$
, where
$\kappa$
is the thermal conductivity,
$\rho$
the density and
$c_p$
the specific heat capacity (Dante Reference Dante2015). In underground reservoirs the matrix is typically sandstone. Measurements by Labus & Labus (Reference Labus and Labus2018) give sandstone effusivity
$e_s \approx 1400$
–
$3700\,\mathrm{J\,m^{-2}\,K^{-1}\,s^{-1/2}}$
. For hydrogen, typical reservoir conditions yield order-of-magnitude estimates
$c_p$
$\sim 10\,\mathrm{kJ\,kg^{-1}\,K^{-1}}$
and
$\kappa \sim 0.1\,\mathrm{W\,m^{-1}\,K^{-1}}$
(Lemmon et al. Reference Lemmon, Bell, Huber, McLinden, Linstrom and Mallard2023), giving hydrogen effusivity
$e_g \sim 10$
–
$10^2\,\mathrm{J\,m^{-2}\,K^{-1}\,s^{-1/2}}$
. Since
$e_g \ll e_s$
, the solid acts as an efficient heat sink, absorbing heat rapidly and keeping the gas temperature nearly constant over the long compression time scale
$t_r$
. However, this may not hold for CO
$_2$
due to its relatively high density. Kushnir et al. (Reference Kushnir, Ullmann and Dayan2012b
) further noted that in high-effusivity rocks, temperature variations can generate thermal stresses that modify pore structure and permeability, while conductive heat loss into the matrix increases gas density and may thereby enhance storage efficiency.
Since shear viscosity and bulk viscosity are pressure- and temperature-dependent, the spreading dynamics may be influenced by thermodynamic effects arising from gas compression. However, the significance of these effects depends upon the gas under consideration. For the range of typical pressures and temperatures of underground reservoirs (
$T \lt 150\,\rm^{\circ }C$
,
$P \lt 50\,\mathrm{MPa}$
), hydrogen exhibits minimal variation in viscosity. We therefore expect that this neglected variation may affect the quantitative agreement with the model, but it is unlikely to substantially alter the overall dynamics or the specific flow regime. In contrast, CO
$_2$
can experience an order-of-magnitude increase in viscosity at elevated pressures (Heinemann et al. Reference Heinemann2021). For CO
$_2$
, this would generate a viscosity gradient along the channel, modifying both flow velocities and the pressure required to drive the fluid.
A distinguishing feature of hydrogen storage, relative to other gases, is the approximately seasonal cycle of injection and withdrawal. Such periodic injection behaviour has been explored in the context of underground thermal energy storage, where hot, dense fluids are cyclically injected and extracted (Dudfield & Woods Reference Dudfield and Woods2014), and more recently, adapted to hydrogen storage scenarios (Whelan & Woods Reference Whelan and Woods2025). Notably, Whelan & Woods (Reference Whelan and Woods2025) showed that when the injection velocity exceeds the characteristic buoyancy speed, the system approaches a state in which oscillations remain confined near the source. Conversely, when buoyancy dominates, oscillations propagate throughout the domain. Although the present study focused mainly on a steady injection rate, the model formulation readily accommodates periodic forcing. As injection cycles drive pressure fluctuations, the associated compression and decompression of the gas will also oscillate. Understanding how the phase relationship between these oscillatory components influences the large-scale flow may be a fruitful direction for future investigation.
In practical settings, stratification and curvature may also play important roles. Curvature of the confining boundaries due to tectonic deformation may be more favourable for trapping buoyant gas (Hagemann et al. Reference Hagemann, Rasoulzadeh, Panfilov, Ganzer and Reitenbach2015). While we have assumed permeability to be isotropic and uniform, natural reservoirs typically exhibit stratification, with porous channels forming layered structures that result in higher horizontal than vertical permeability. These directional variations are further compounded by spatial variation in the permeability arising from the varying pore structure, which can significantly influence large-scale flow behaviour. Additionally, the present model is readily generalisable to an axisymmetric geometry with a point source. A key qualitative difference from the planar case is that, depending on a viscous/buoyancy ratio (Guo et al. Reference Guo, Zheng, Celia and Stone2016; Hutchinson, Gusinow & Worster Reference Hutchinson, Gusinow and Worster2023), the injected gas plume may not extend to the lower domain boundary, a consequence of the point source (Guo et al. Reference Guo, Zheng, Celia and Stone2016).
As noted in the introduction, Cuttle et al. (Reference Cuttle, Morrow and MacMinn2023) studied a mathematically equivalent set-up involving liquid displacement in a radial Hele-Shaw cell, focusing on how gas compressibility influences the onset of interfacial instabilities. Our findings are in qualitative agreement with their observations of bulk flow dynamics in the absence of such instabilities. Specifically, they reported a non-monotonic pressure evolution, i.e. an initial rise due to gas compression, followed by a decline as viscous resistance diminished due to liquid drainage. They also found that the pressure scale decreased with increasing compressibility and that breakthrough times were correspondingly delayed. Notably, these bulk features were relatively insensitive to the presence of instabilities and exhibited only weak dependence on the capillary number, suggesting that the overall dynamics is primarily governed by a balance between gas compressibility and viscous dissipation in the displaced liquid. Building on these insights, we incorporated buoyancy effects and viscous dissipation in the gas. In contrast to their model, where the gas pressure is assumed spatially uniform and the dynamics is entirely coupled to viscous dissipation in the liquid, we revealed that, in the high-gas-mobility limit
$\mathcal{M}\ll 1$
, the flow develops an inner–outer structure. Near the source, for
$\mathcal{L}^{-1}\ll \mathcal{M}\ll 1$
, the gas pressure equilibrates and the flow is incompressible at late times, while in the outer region, compressibility and viscous dissipation govern the spreading; for smaller
$\mathcal{M}$
, the bubble pressure can continue to rise throughout spreading.
In summary, we have investigated the influence of compressibility on flow dynamics using a long-wave model, exploring a broad range of parameter space. We presented reduced models and scalings for the dominant balances across distinct regions of parameter space, testing reduced-order models against numerical results. Our findings show that gas compression slows the spreading rate and increases the mass-to-volume ratio of the gas plume by the end of the injection cycle, compared with the incompressible case. These effects are particularly relevant for underground gas-storage applications, such as hydrogen storage, where accurate plume evolution predictions are crucial. Practically, our results suggest that elongated aquifers (
$H/L \ll 1$
), with relatively low permeability, amplify compressibility effects, leading to more efficient storage. We have demonstrated that, even for weakly compressible gases, compressibility effects may not be purely transient, but can persist throughout the entire spreading process when the channel is sufficiently long.
Acknowledgements
LM would like to acknowledge the funding support from the EPSRC Industrial Decarbonisation Research and Innovation Centre (IDRIC MIPs 7.4 & 7.8; EP/V027050/1). For the purpose of open access, the authors have applied a Creative Commons Attribution (CC-BY) licence to any author accepted manuscript version arising. The authors acknowledge numerous helpful comments from reviewers.
Declaration of interests
The authors report no conflict of interest.
Data availability statement
The MATLAB code required to reproduce the numerical simulations of this study is available at https://github.com/PeterCastellucci/Compressible_Evolution_Equations.
Appendix A. Numerical methods
To solve (2.29), we employ a coordinate transformation, mapping
$(0,X_l(t))$
and
$(X_l(t), X_u(t))$
onto fixed domains. We use the method of lines, discretising space using a finite-volume method based on central differences and integrating the resulting system of ordinary differential equations (ODEs) in time using MathWorks MATLAB R2024b’s stiff solver ode15s. This solver employs variable-step, variable-order integration, resulting in a scheme that achieves second-order accuracy in both space and time. For each simulation, the error in global mass conservation is maintained below 1 %. To enhance computational efficiency, we numerically compute the Jacobian matrix and supply it explicitly to the solver. By avoiding the overhead of the solver recalculating the Jacobian at each iteration, this approach significantly reduces computation time.
In the early time formulation, where the compressed pressure and dynamic pressure are decoupled, we time step
$P_0(t)$
using (2.35), while the interface height
$F(x,t)$
is evolved independently using (2.34b
). At each time step, the solution of (2.35) is substituted into (2.34), allowing us to solve for
$P_{1x}$
via the resulting ODE.
When solving the outer problem
$\varPi _{2}$
, we solve on the domain
$\bar {x} = (\delta , \check {X}_u(\check {t}))$
, where
$\delta \ll 1$
to capture the singular behaviour from the boundary condition (2.44e
). In order to improve the resolution of the singularity at the origin, we use the transformation
$\bar {x} = \exp (\xi ) \check {X}_u$
. A uniform grid in
$\xi$
then clusters grid points near the origin. We discretise in space using a second-order-accurate finite-volume method and use ode15s to time step these equations.
Appendix B. Problem
$\boldsymbol{\varPi _b}$
: Buoyancy
A distinguished limit appears when
$\phi \equiv \mathcal{L} \mathcal{M}^2\sim 1$
and
$\theta \equiv \zeta \mathcal{L} \sim \mathcal{M} \ll 1$
(problem
$\varPi _b$
, figure 3). Here, buoyancy effects become dominant in the inner region of problem
$\varPi _{12}$
, with (2.38) replaced by
\begin{align} 0& = \left [ (1-\tilde {F}) \tilde {P}_0 \tilde {P}_{1\tilde {x}}\right ]_{\tilde {x}}, \\[-9pt] \nonumber \end{align}
\begin{align} \tilde {F}_{\tilde {t}} &=\left [ \tilde {F}(\tilde {P}_{1\tilde {x}}+\phi ^{-1}\tilde {F}_{\tilde {x}})\right ]_{\tilde {x}}, \quad \tilde {X}_l\lt \tilde {x}, \\[-9pt] \nonumber \end{align}
\begin{align} - \tilde {P}_0 \tilde {P}_{1\tilde {x}}&=\mathcal{Q}(\mathcal{M}/\theta ) \quad (\tilde {x}=0), \\[-9pt] \nonumber \end{align}
\begin{align} \tilde {X}_{l,\tilde {t}}&=- (\tilde {P}_{\tilde {x}}+\phi ^{-1}\tilde {F}_{\tilde {x}})\quad ( \tilde {x}={\tilde {X}_l+}). \\[9pt] \nonumber \end{align}
The interface evolution is governed by the mixed hyperbolic–parabolic evolution equation
\begin{align} \tilde {F}_{\tilde {t}}+\frac {\mathcal{Q}}{\tilde {P}_0}\frac {\mathcal{M}}{\theta }\frac {\tilde {F}_{\tilde {x}}}{(1-\tilde {F})^2}=\frac {1}{\phi }\left [\tilde {F}\tilde {F}_{\tilde {x}}\right ]_{\tilde {x}}. \end{align}
This corresponds to (2.38a
–e) with buoyancy effects restored. For
$\mathcal{L}^{-1/2}\ll \mathcal{M}\ll 1$
, buoyancy effects are confined to a boundary layer near the lower contact line, perturbing its location. For
$\phi \ll 1$
, the inner region splits into two zones: a short region in which
$\tilde {F}_{\tilde {t}}=\phi ^{-1}(\tilde {F}\tilde {F}_{\tilde {x}})_{\tilde {x}}$
to leading order, with the lower contact line receding towards the origin; and a longer unsteady region of mixed parabolic–hyperbolic form with
$F$
close to 1 that allows matching to the outer region. We write the latter region using
$\tilde {F}=1-\tilde {f}(\tilde {x},\tilde {t})$
with
$\vert \tilde {f}\vert \ll 1$
, so that to leading order
\begin{align} \tilde {f}_{\tilde {t}}+\frac {\mathcal{Q}}{\tilde {P}_0}\frac {\mathcal{M}}{\theta }\frac {\tilde {f}_{\tilde {x}}}{\tilde {f}^2}=\frac {1}{\phi }\tilde {f}_{\tilde {x}\tilde {x}}. \end{align}
The outer region matches via the first two terms, where
$\tilde {f}$
becomes very small over long length scales, while the nonlinear diffusion region matches to (B2) via the first and third terms where
$\tilde {f}$
is larger. Equation (B2) admits a similarity solution involving a balance of all three terms with
$\tilde {f}=(\mathcal{Q}\mathcal{M}/\tilde {P}_0\theta )^{1/2} (\phi t^{1/4})g(\xi )$
,
$\xi =(\phi /t)^{1/2} x$
provided the time-dependence of
$\mathcal{Q}$
and
$\tilde {P}_0$
is neglected. However, we use (B1a
) to argue that the gas flux through this region is uniform, allowing the source (B1c
) to communicate directly to the outer region, so that recession of the lower contact line does not have a prominent effect on spreading further downstream.
Appendix C. Very low viscosity ratio
Here, we discuss features that emerge at the base of regions
$\varPi _1$
and
$\varPi _2$
in figure 3, where the viscosity ratio is very low. When
$\mathcal{M}$
is very small in the early time problem (2.34),
$F_0$
and
$X_l$
are effectively frozen at their initial state (2.22) and (2.34a
,
c
) can be integrated to give
\begin{align} P_{1x}=\frac {\Delta _0(P_{0t}L_0-\mathcal{Q})}{P_0(X_u(0)-x)}-\frac {X_u(0)-x}{2\Delta _0}\frac {P_{0t}}{P_0}. \end{align}
This solution accommodates a mass flux
$P_0 P_{1x}(1-F)=P_{0t}L_0-\mathcal{Q}$
as
$x\rightarrow X_u(0)-$
, consistent with (2.35) for
$\mathcal{M}\ll 1$
. Here,
$P_1$
diverges logarithmically as
$x\rightarrow X_u(0)-$
and is regularised by viscous effects over short length scales as follows.
Set
$x=X_u(0)+\mathcal{M}^{1/3}\breve {x}$
,
$F=1-\mathcal{M}^{1/3}\breve {F}(\breve {x},t)$
. Then (2.34) gives, to leading order in
$\mathcal{M}\ll 1$
,
\begin{align} 0=[\breve {F}P_{1\breve {x}}]_{\breve {x}}, \quad \breve {F}_t=-P_{1\breve {x}\breve {x}}. \end{align}
This matches to the frozen region via
$\breve {F}\approx -\breve {x}/\Delta _0$
for
$-\breve {x}\gg 1$
. The mass flux of gas is uniform through this short region, with value
$\breve {F}P_0P_{1\breve {x}}\approx P_{0t}L_0-\mathcal{Q}$
; (C2) gives the local evolution equation
\begin{align} \breve {F}_t+\frac {\mathcal{Q}-P_{0t}L_0}{P_0}\frac {\breve {F}_{\breve {x}}}{\breve {F}^2}=0. \end{align}
This variation of the hyperbolic transport problem for the interface shape in the inner region (2.38e
) accounts for compressible effects through transient variations in
$P_0$
, arising over short time scales. The solution of (C3) satisfying
$\breve {F}(\breve {x},0)=-\breve {x}/\Delta _0$
in
$\breve {x}\lt 0$
can be evaluated (via characteristics) as
\begin{align} \breve {F}^2(\breve {x}-\breve {F}\Delta _0)=\int _0^t \frac {\mathcal{Q}}{P_0}\,\mathrm{d}t-L_0\log \left (\frac {P_0(t)}{P_0(0)}\right ). \end{align}
This links the upstream limit
$\breve {F}\approx -\breve {x}/\Delta _0\rightarrow \infty$
with the downstream limit
\begin{align} {F}\approx 1- \left (\frac {\mathcal{M}}{x-X_u(0)}\right )^{1/2} \left [\int _0^t \frac {\mathcal{Q}}{P_0}\,\mathrm{d}t-L_0\log \left (\frac {P_0(t)}{P_0(0)}\right ) \right ]^{1/2} \end{align}
(expressed in the original variables; recall
$P=P_0/\zeta +P_1+\ldots$
). Thus, over short time scales, the outer region is supplied directly by a flux that accommodates dynamic variations in
$P_0$
. This adjustment can be expected to have an influence along the line
$\mathcal{M}\zeta \mathcal{L}^3\sim 1$
, when the time scale of problem
$\varPi _2$
is no longer long in comparison to the time scale of the early evolution. This line radiates from the distinguished limit
$\zeta \sim 1/\mathcal{L}^2$
,
$\mathcal{M}\sim 1/\mathcal{L}$
.
Problem
$\varPi _c$
: we therefore reformulate (2.29) setting
$P=\mathcal{L}^2\ddot {P}(\bar {x},t)$
,
$F=1-\mathcal{M}\ddot {F}(\bar {x},t)$
,
$x = \mathcal{L} \bar {x}$
, to give (for
$\mathcal{L}\gg 1$
)
\begin{align} \left [\ddot {F}\ddot {P}\right ]_t &= \left [\ddot {F} \ddot {P} \ddot {P}_{\bar {x}}\right ]_{\bar {x}} \quad (0\lt \bar {x}\lt \bar {X}_u), \\[-9pt] \nonumber \end{align}
\begin{align} -\ddot {F}_t& = \ddot {P}_{\bar {x}\bar {x}}\quad (0\lt \bar {x}\lt \bar {X}_u), \\[-9pt] \nonumber\end{align}
\begin{align} \ddot {P}+\ddot {P}_{\bar {x}} (1 - \bar {X}_u)& =\frac {1}{\zeta \mathcal{L}^2}\quad (\bar {x}=\bar {X}_u), \\[-9pt] \nonumber \end{align}
\begin{align} \bar {X}_{u,t}&=- \ddot {P}_{\bar {x}}\quad (\bar {x}=\bar {X}_u-), \\[9pt] \nonumber \end{align}
with
\begin{align} \ddot {F}\approx \left [\frac {1}{\bar {x}} \left ( \frac {1}{\mathcal{M}\zeta \mathcal{L}^3}\int _0^t \frac {\mathcal{Q}}{\ddot {P}(0,t)}\,\mathrm{d}t-\frac {L_0}{\mathcal{M}\mathcal{L}}\log \left (\frac {\zeta \mathcal{L}^2 \ddot {P}(0,t)}{P_0(0)}\right ) \right )\right ]^{1/2} \quad (\bar {x}\rightarrow 0). \end{align}
This problem (
$\varPi _c$
in figure 3) has two independent parameters,
$\zeta \mathcal{L}^2$
and
$\mathcal{M}\mathcal{L}$
, which are both
$O(1)$
in the distinguished limit. We can state (C6e
) alternatively as
\begin{align} \mathcal{ML}(\ddot {F}^2\bar {x})_t \approx \frac {\mathcal{Q}/(\zeta \mathcal{L}^2)-{L_0} \ddot {P}_t}{\ddot {P}} \quad (\bar {x}\rightarrow 0). \end{align}
Here the flux entering the thin film (on the left) balances the source flux and dynamic (compressive) changes to the pressure in the main gas bubble. The parameter
$\mathcal{M}\mathcal{L}$
regulates the flux into the thin film while the parameter
$\zeta \mathcal{L}^2$
regulates whether spreading is driven by the mass flux at the source or by the evolving bubble pressure.
When
$1\ll \mathcal{ML}\sim \zeta \mathcal{L}^2\ll \mathcal{L}$
, we recover from (C6) the outer limit (2.39) of problem
$\varPi _{12}$
by setting
$\ddot {P}=\hat {P}/(\mathcal{LM})$
and
$t=\mathcal{LM}\tilde {t}$
. The log term in (C6e
) provides an
$O(1/(\mathcal{ML}))$
correction.
When
$\max ((\zeta \mathcal{L}^2)^{1/2},\zeta \mathcal{L}^2)\ll \mathcal{ML}\ll \mathcal{L}$
, we set
$\ddot {P}=(\zeta \mathcal{L}^2)^{-1}+(\mathcal{ML})^{-1}\hat {P}_1$
and
$t=\mathcal{ML}\tilde {t}$
recovers (2.40) from (C6), i.e. the outer limit of problem
$\varPi _1$
, with error
$O(L_0\zeta /\mathcal{M}^2)$
.
When
$1\ll \mathcal{ML}\ll \zeta \mathcal{L}^2$
, we recover from (C6) the outer problem (2.44) (with error
$1/(\zeta \mathcal{L}^2)$
) by setting
$\ddot {P}=\check {P}/(\mathcal{ML}^3\zeta )^{1/2}$
,
$t=(\mathcal{ML}^3\zeta )^{1/2}\check {t}$
, but with the revised source condition
\begin{align} \check {F}\approx \left [\frac {1}{\bar {x}} \left (\int _0^{\check {t}} \frac {\mathcal{Q}}{\check {P}(0,\check {t}^{\prime })}\,\mathrm{d}\check {t}^{\prime }-\frac {L_0}{\mathcal{M}\mathcal{L}}\log \left (\frac {(\mathcal{ML})^{-1}\check {P}(0,\check {t})}{P_0(0)}\right ) \right )\right ]^{1/2}. \end{align}
The log term identifies the lower boundary of problem
$\varPi _2$
as the line
$\mathcal{ML}\sim 1$
. Along this boundary, the two terms in (C8) can be expected to balance.
Below this boundary, we set
$t=(\zeta \mathcal{L}^2)^{1/2}\dot {t}$
,
$\ddot {P}=(\zeta \mathcal{L}^2)^{-1/2}\dot {P}$
,
$\bar {X}_u(t) = \dot {X}_u(\dot {t})$
. This recovers from (C6) the problem
\begin{align} \left [\dot {F}\dot {P}\right ]_{\dot {t}}& = \left [\dot {F} \dot {P} \dot {P}_{\bar {x}}\right ]_{\bar {x}}\quad\ (0\lt \bar {x}\lt \dot {X}_u), \\[-9pt] \nonumber \end{align}
\begin{align} -\dot {F}_{\dot {t}} &= \dot {P}_{\bar {x}\bar {x}}\qquad\qquad\! (0\lt \bar {x}\lt \dot {X}_u), \\[-9pt] \nonumber \end{align}
\begin{align} \dot {P}+\dot {P}_{\bar {x}} (1 - { \dot {X}_u})& =0\qquad\qquad\quad\kern-2.3pt (\bar {x}=\dot {X}_u), \\[-9pt] \nonumber \end{align}
\begin{align} {\dot {X}_{u,\dot {t}}}&=- \dot {P}_{\bar {x}}\qquad\qquad\kern-5.8pt (\bar {x}=\dot {X}_u-), \\[-9pt] \nonumber \end{align}
\begin{align} \dot {P}_{\dot {t}}&=\mathcal{Q}/{L_0}\qquad\quad\kern-3.7pt\ (\bar {x}=0), \\[9pt] \nonumber \end{align}
which holds for
$\mathcal{ML}\ll 1$
and
$1\ll \zeta \mathcal{L}^2$
(region
$\varPi _3$
in figure 3). The pressure at the origin rises independently of the degree of spreading via (C9e
), because there is very small leakage of gas from the bubble into the thin film. However, we have now lost a boundary condition on
$\dot {F}$
. This is resolved via retention of buoyancy, for example, in (C3), which becomes
\begin{align} \breve {F}_t+\frac {\mathcal{Q}-P_{0t}L_0}{P_0}\frac {\breve {F}_{\breve {x}}}{\breve {F}^2}=\mathcal{M}^{1/3}\breve {F}_{\breve {x}\breve {x}}. \end{align}
If the film becomes sufficiently thin, the flux it carries is too small to allow gas to escape the bubble and the supply is balanced by a continuous increase in pressure. In this state, the advection term in (C10) is then very weak, allowing buoyancy to regulate the transition between the bubble and the thin gas film further downstream. Equation (C10) suggests that, under (C9e
), buoyancy becomes dominant over a short length scale
$x\sim \mathcal{M}^{1/2}$
near
$X_u(0)$
and can be expected to modify the matching condition on
$\dot {F}$
.
On the boundary
$\zeta \mathcal{L}^2\sim 1$
with
$\mathcal{ML}\ll 1$
, (C6) becomes
\begin{align} \left [\ddot {F}\ddot {P}\right ]_t& = \left [\ddot {F} \ddot {P} \ddot {P}_{\bar {x}}\right ]_{\bar {x}}\quad\ (0\lt \bar {x}\lt \bar {X}_u), \\[-9pt] \nonumber \end{align}
\begin{align} -\ddot {F}_t& = \ddot {P}_{\bar {x}\bar {x}}\qquad\qquad\! (0\lt \bar {x}\lt \bar {X}_u), \\[-9pt] \nonumber \end{align}
\begin{align} \ddot {P}+\ddot {P}_{\bar {x}} (1 - \bar {X}_u)& ={1}/({\zeta \mathcal{L}^2})\qquad (\bar {x}=\bar {X}_u), \\[-9pt] \nonumber\end{align}
\begin{align} \bar {X}_{u,t}&=- \ddot {P}_{\bar {x}}\qquad\qquad\kern-5.8pt (\bar {x}=\bar {X}_u-), \\[-9pt] \nonumber \end{align}
\begin{align} \ddot {P}_t&={\mathcal{Q}}/(L_0 \zeta \mathcal{L}^2)\quad\kern-7.7pt\ (x\rightarrow 0). \\[9pt] \nonumber \end{align}
This is distinguished from (C9) by the additional term in (C11c ), elevating the pressure.
The final major sublimit of (C6) arises for
$\mathcal{M}\ll \zeta ^{1/2}$
and
$\zeta \mathcal{L}^2\ll 1$
(region
$\varPi _4$
in figure 3). Here, set
$\ddot {P}=(\zeta \mathcal{L}^2)^{-1}+\mathring {P}_1/(\zeta \mathcal{L}^2)^{1/2}$
,
$\ddot {F}=\mathring {F}(\bar {x},\mathring {t})$
,
$t=(\zeta \mathcal{L}^2)^{1/2}\mathring {t}$
and
$\bar {X}_u (t)= \mathring {X}_{u}(\mathring {t})$
, giving with error
$O(\mathcal{M}/\zeta ^{1/2})$
,
\begin{align} \mathring {F}_{\mathring {t}} &= \left [\mathring {F} \mathring {P}_{1\bar {x}}\right ]_{\mathring {x}}\quad\ (0\lt \bar {x}\lt \mathring {X}_u), \\[-9pt] \nonumber \end{align}
\begin{align} -\mathring {F}_{\mathring {t}}& = \mathring {P}_{1\bar {x}\bar {x}}\qquad\qquad\kern-11.2pt (0\lt \bar {x}\lt \mathring {X}_u), \\[-9pt] \nonumber \end{align}
\begin{align} \mathring {P}_1 + \mathring {P}_{1x} (1 - \mathring {X}_u) &= 0\qquad\qquad\kern3.7pt (\bar {x}=\mathring {X}_u-), \\[-9pt] \nonumber \end{align}
\begin{align} \mathring {X}_{u,\mathring {t}}&=- \mathring {P}_{1\bar {x}}\qquad\kern3.7pt\ (\bar {x}=\mathring {X}_u-), \\[-9pt] \nonumber \end{align}
\begin{align} \mathring {P}_{1\mathring {t}}& =\mathcal{Q}/{L_0}\qquad\qquad\kern-16.7pt (\bar {x}=0). \\[9pt] \nonumber\end{align}
This resembles the incompressible problem
$\varPi _1$
, but is driven by a pressure condition instead of an interfacial condition. Finally, along the line
$\mathcal{M}\sim \zeta ^2\ll 1/\mathcal{L}$
, the
$(F^2 \bar {x})_{\mathring {t}}$
term is restored to the boundary condition (C12); this term balances
$\mathcal{Q}$
on moving back into the incompressible region
$\varPi _1$
. Combining (C12a
,
b
), integrating up to
$\mathring {X}_u$
and applying the boundary conditions (C12d
,
e
) relates the pressure gradient to the velocity of the upper contact line
\begin{align} \mathring {P}_{1x} = - \frac {\mathring {X}_{u,\mathring {t}}}{(1+F)}. \end{align}
Substituting this expression into (C12b ) leads to the hyperbolic problem (2.41), which can be integrated along characteristics to give
\begin{align} \mathring {F} = \left (\frac {\mathring {X}_u(\mathring {t})}{\bar {x}} \right )^{1/2} - 1. \end{align}
This solution has the equivalent form as the incompressible solution (2.42), but here
$\mathring {X}_u$
is determined by the buildup of pressure in the gas bubble via (C12e
), rather than being driven directly by the source. Integrating (C13), applying (C12e
) and using (C14) to evaluate the integral of
$1/(1+\mathring {F})$
, yields
\begin{align} \mathring {P} = \frac {\mathring {\mathcal{V}}(\mathring {t})}{L_0} -\frac {2 \bar {x}^{3/2}}{3\mathring {X}_u^{1/2}} \mathring {X}_{u,\mathring {t}}, \quad \mathring {\mathcal{V}}(\mathring {t}) = \int _0^{\mathring {t}} \mathcal{Q} \,\mathrm{d}\mathring {t}^\prime . \end{align}
Combining the boundary condition (C12c
) with (C15) gives
$\mathring {X}_{u,\mathring {t}}(1-( {1}/{3})\mathring {X}_u)=\mathring {\mathcal{V}}/L_0$
, with the solution
\begin{align} \mathring {X}_u(\mathring {t}) = 3\left [ 1- \left (1 - \frac {2}{3 L_0} \int _0^{\mathring {t}} \mathring {\mathcal{V}}(\mathring {t}^\prime ) \,\mathrm{d} \mathring {t}^{\prime } \right )^{1/2} \right ]. \end{align}
Taking
$\mathcal{Q} = 1$
, (C16) reveals that, for small
$\mathring {t}$
,
$\mathring {X}_u \approx \mathring {t}^2/(2 L_0)$
or equivalently
$X_u\approx t^2/(2\zeta \mathcal{L}L_0)$
. Equation (C16) predicts a breakthrough time of
$t_b = ((5/3) L_0 \zeta \mathcal{L}^2)^{1/2}$
.
Open access
