2. Theoretical framework

2.1. Problem description

We consider a semi-infinite, unbounded, homogeneous porous medium characterised by its porosity $\phi$ and permeability $K$ , with a rigid and impermeable curvilinear lower surface $f(x)$ (see figure 1). In this system, the characteristic pore-scale length $l_{p}$ is assumed to be significantly smaller than the characteristic length of the curvilinear surface $l_{c}$ (i.e. $l_{p} \ll l_{c}$ ), ensuring a clear scale separation. Initially, the domain is saturated with an ambient fluid of density $\rho _{0}$ and viscosity $\mu _{0}$ . At the origin, a finite volume of a second fluid, which is denser (density $\rho$ ) and more viscous (viscosity $\mu$ ) than the ambient fluid, is released. This will create a porous media gravity current driven by the density difference $\Delta \rho =\rho -\rho _{0}$ . Under the assumption of negligible macroscopic capillary effects and the absence of Saffman–Taylor instability at the interface, the denser fluid maintains a sharp, well-defined interface with the ambient fluid as it propagates in the creeping flow regime.

Figure 1.

Schematic diagram of a gravity current in an homogeneous porous medium over an impermeable and rigid curvilinear surface.

2.2. Low-dimensional model for local evolution

2.2.1. Pressure and velocity fields

We first establish the hydrostatic pressure field in the local coordinates $(s,n)$ , where $s$ denotes the tangential direction along the curvilinear surface and $n$ is the direction normal to it. For convenience, we neglect the effect of the lighter ambient fluid and express the pressure as

(2.1) \begin{equation} P = p_0 + \Delta \rho g\,\textrm {cos}\,\theta (h_n + f_n - n), \end{equation}

where $p_0$ is a reference pressure and $\theta$ is the angle between the local normal direction $n$ and the global vertical direction $z$ . Here, $h_n$ denotes the fluid thickness measured along the normal direction and $f_n = f(s)/\textrm {cos}\,\theta$ represents the projection of the surface elevation (given by $f(x)$ in the global coordinate) onto the normal direction.

The local pressure gradient along the tangential direction is thus given by

(2.2) \begin{equation} \frac {\partial P}{\partial s} = \underbrace {\Delta \rho g\,\textrm {cos}\,\theta \frac {\partial h_n}{\partial s}}_{\text{Interface slope}} + \underbrace {\Delta \rho g\,\textrm {cos}\,\theta \frac {\partial f_n}{\partial s}}_{\text{Surface slope}} + \underbrace {\Delta \rho g (h_n + f_n - n)\frac {\partial }{\partial s}(\textrm {cos}\,\theta )}_{\text{Surface curvature}}, \end{equation}

showing that the local pressure gradient comprises contributions from the interface slope, the surface slope and the curvature of the surface.

To evaluate the surface slope term, we apply the chain rule

(2.3) \begin{equation} \frac {\partial f_n}{\partial s} = \frac {1}{\textrm {cos}\,\theta }\frac {\partial f}{\partial s} + f(s)\frac {\partial }{\partial s}\left (\frac {1}{\textrm {cos}\,\theta }\right ) = \frac {1}{\textrm {cos}\,\theta }\frac {\partial f}{\partial s} + f(s)\left (\frac {\sin \theta }{\textrm {cos}^2\,\theta }\right )\frac {\partial \theta }{\partial s}. \end{equation}

Introducing the local surface curvature defined as

(2.4) \begin{align} \kappa (s) = \frac {\partial \theta }{\partial s}, \end{align}

this expression becomes

(2.5) \begin{equation} \frac {\partial f_n}{\partial s} = \frac {1}{\textrm {cos}\,\theta }\frac {\partial f}{\partial s} + \frac {f(s) \sin \theta \kappa (s)}{\textrm {cos}^2\,\theta }. \end{equation}

For the surface curvature term, the derivative of $\textrm {cos}\,\theta$ with respect to $s$ is given by

(2.6) \begin{equation} \frac {\partial }{\partial s}(\textrm {cos}\,\theta ) = -\sin \theta \frac {\partial \theta }{\partial s} = -\sin \theta \kappa (s). \end{equation}

Thus, the final expression for the tangential pressure gradient is

(2.7) \begin{equation} \frac {\partial P}{\partial s} = \Delta \rho g\,\textrm {cos}\,\theta \left (\frac {\partial h_n}{\partial s} + \frac {1}{\textrm {cos}\,\theta }\frac {\partial f}{\partial s} + \frac {f\sin \theta \kappa (s)}{\textrm {cos}^2\,\theta }\right ) - \Delta \rho g(h_n + f_n - n)\sin \theta \kappa (s). \end{equation}

In the normal direction, the pressure is assumed to be in hydrostatic balance

(2.8) \begin{equation} \frac {\partial P}{\partial n} = -\Delta \rho g\,\textrm {cos}\,\theta . \end{equation}

This hydrostatic condition implies that the pressure gradient in the $n$ direction serves only to maintain local static equilibrium, rather than driving fluid motion with a significant velocity in the normal direction. Consequently, the flow is effectively confined to the tangential direction, allowing us to neglect the contribution of the normal pressure gradient.

Substituting the tangential pressure gradient into Darcy’s law for the creeping flow regime,

(2.9) \begin{equation} u_s = -\frac {K}{\mu }\frac {\partial P}{\partial s}, \end{equation}

we obtain the local velocity field

(2.10) \begin{align} u_s = -\frac {K\Delta \rho g}{\mu } \left [ \textrm {cos}\,\theta \left (\frac {\partial h_n}{\partial s} + \frac {1}{\textrm {cos}\,\theta }\frac {\partial f}{\partial s} + \frac {f \sin \theta \kappa (s)}{\textrm {cos}^2\,\theta }\right ) - \sin \theta \kappa (s)\left (h_n + \frac {f}{\textrm {cos}\,\theta } - n\right ) \right ].\nonumber\\[2pt] \end{align}

That is,

(2.11) \begin{align} u_s = -\frac {K\Delta \rho g}{\mu } \!\left [\! \underbrace {\textrm {cos}\,\theta \frac {\partial h_n}{\partial s}}_{\text{Interface slope}} \! + \underbrace {\frac {\partial f}{\partial s}}_{\text{Surface slope}} \! + \underbrace {\left (\!\frac {f\sin \theta \kappa (s)}{\textrm {cos}\,\theta } - \sin \theta \kappa (s)\left (\!h_n + \frac {f}{\textrm {cos}\,\theta } - n \!\right )\! \right )}_{\text{Curvature-driven terms}}\! \right ].\nonumber\\[2pt] \end{align}

2.2.2. Local evolution

We substitute the local velocity field into the local mass conservation equation

(2.12) \begin{equation} \phi \frac {\partial h_n}{\partial t} + \frac {\partial }{\partial s}\int _{f_n}^{h_n+f_n} u_s\, {\rm d}n = 0. \end{equation}

Evaluating the integral term, we obtain

(2.13) \begin{equation} \int _{f_n}^{h_n+f_n} \textrm {cos}\,\theta \frac {\partial h_n}{\partial s}\, {\rm d}n = h_n\, \textrm {cos}\,\theta \frac {\partial h_n}{\partial s}, \end{equation}

(2.14) \begin{equation} \int _{f_n}^{h_n+f_n} \frac {\partial f}{\partial s}\,{\rm d}n = h_n\frac {\partial f}{\partial s}, \end{equation}

(2.15) \begin{equation} \int _{f_n}^{h_n+f_n} \left [\frac {f \sin \theta \kappa (s)}{\textrm {cos}\,\theta } - \sin \theta \kappa (s)(h_n + f_n - n)\right ]\, {\rm d}n = \sin \theta \kappa (s)\left (\frac {fh_n}{\textrm {cos}\,\theta } - \frac {h_n^2}{2}\right ). \end{equation}

Thus, the local evolution equation is given by

(2.16) \begin{equation} \phi \frac {\partial h_n}{\partial t} = \frac {K\Delta \rho g}{\mu } \frac {\partial }{\partial s} \left [ h_n \textrm {cos}\,\theta \frac {\partial h_n}{\partial s} + h_n\frac {\partial f}{\partial s} + \sin \theta \kappa (s)\left ( \frac {fh_n}{\textrm {cos}\,\theta } - \frac {h_n^2}{2} \right ) \right ]. \end{equation}

This local evolution equation reveals that the fluid dynamics is governed by three distinct physical mechanisms. First, the term $h_n\textrm {cos}\,\theta \partial h_n/\partial s$ represents classical gravity-driven spreading, where gradients in fluid thickness induce pressure differences that drive the flow. Second, the term $h_n\partial f/\partial s$ embodies topographic forcing, whereby variations in the curvilinear surface elevation generate additional pressure gradients to drive the flow. Third, the curvature-driven terms, represented by $\sin \theta \kappa (s)(f h_n/\textrm {cos}\,\theta - h_n^2/2)$ , capture more effects arising from the interplay between fluid and surface geometry. Therein, the component $f h_n/\textrm {cos}\,\theta$ represents a centripetal effect due to fluid motion along curved streamlines. This effect tends to resist the flow in convex regions (where $\kappa \gt 0$ ) and to assist it in concave regions (where $\kappa \lt 0$ ). The term $-h_n^2/2$ accounts for hydrostatic pressure variations induced by thickness gradients in regions of curvature, tending to drive the flow towards areas of lower curvature. These curvature-induced effects can lead to complex flow behaviours near transitions between convex and concave sections.

2.3. Towards global evolution

2.3.1. Current height in global coordinates

While the local evolution equation effectively describes fluid motion relative to the curvilinear surface, practical applications are primarily concerned with the overall global performance. Therefore, it is essential to transform the governing equation from local to global coordinates.

We begin by relating the vertical fluid thickness, $h_z$ , measured vertically from the curvilinear surface to the fluid interface, to the thickness measured along the surface normal, $h_n$ . Since $h_z$ is naturally aligned with gravity and the global coordinate system, it is advantageous to express $h_n$ in terms of $h_z$ . Their relationship is given by

(2.17) \begin{equation} h_z = h_n\textrm {cos}\,\theta + \beta \kappa _z h_n^2 + {\mathcal{O}}(h_n^3), \end{equation}

where the first term, $h_n\textrm {cos}\,\theta$ , represents the geometric projection from the normal to the vertical direction, and the second term, $\beta \kappa _z h_n^2$ , is a curvature correction that accounts for the bending of the interface. Here, $\kappa _z$ denotes the local curvature of the interface, defined as

(2.18) \begin{equation} \kappa _z = \frac {\partial ^2 h_z/\partial x^2}{\left [1+\left (\partial h_z/\partial x\right )^2\right ]^{3/2}}. \end{equation}

Under the thin-film approximation, i.e. $\partial h_z/\partial x \sim h_n/L_c \ll 1$ , where $L_c$ is the characteristic horizontal length scale, the interface curvature simplifies to

(2.19) \begin{equation} \kappa _z \approx \frac {\partial ^2 h_z}{\partial x^2}\sim \frac {h_n}{L_c^2} \ll 1. \end{equation}

The dimensionless parameter $\beta$ (typically $\beta \sim {\mathcal{O}}(1)$ ) quantifies this first-order curvature correction. As a result, the second term is negligible compared with the leading-order term $h_n$ and we are justified in using the simplified relation

(2.20) \begin{equation} h_z = h_n\textrm {cos}\,\theta , \end{equation}

which preserves the essential geometric projection while neglecting higher-order effects that have minimal influence on the overall flow behaviour in the thin-film regime. Since the surface is assumed to be rigid (i.e. time-independent), the time derivative transforms as

(2.21) \begin{equation} \frac {\partial h_n}{\partial t} = \frac {1}{\textrm {cos}\,\theta }\frac {\partial h_z}{\partial t}. \end{equation}

2.3.2. Transform to a global evolution model

To assess the relative importance of surface slope versus curvature effects, we non-dimensionalise the two terms by defining the following dimensionless variables:

(2.22) \begin{equation} s^* = s/L_c, \quad f^*(x^*) = f(x)/L_c, \quad h_n^*(x^*) = h_n(x)/L_c, \quad \kappa ^*(s^*) = \kappa (s)L_c. \end{equation}

If the slope of a curvilinear surface is small overall, we introduce a small parameter $\epsilon$ , defined as the maximum value of the dimensionless surface slope,

(2.23) \begin{equation} \epsilon = \left |{\rm d}f^*/ {\rm d}x^*\right | \ll 1, \end{equation}

and for a small angle $\theta$ , we have

(2.24) \begin{equation} \tan \theta \approx \theta \approx |{\rm d}f^*/ {\rm d}x^*|, \end{equation}

which implies that

(2.25) \begin{equation} \theta \sim \epsilon , \end{equation}

and we further obtain

(2.26) \begin{equation} \textrm {cos}\,\theta \approx 1 - \frac {\theta ^2}{2} \approx 1, \quad \sin \theta \approx \theta \sim \epsilon . \end{equation}

Because the curvature is defined as

(2.27) \begin{equation} \kappa = \frac {{\rm d}\theta }{{\rm d}s}, \end{equation}

and over a characteristic length $L_c$ , the typical change in $\theta$ is of order $\epsilon$ , we can estimate

(2.28) \begin{equation} \kappa \sim \frac {\Delta \theta }{L_c} \sim \frac {\epsilon }{L_c}, \end{equation}

(2.29) \begin{equation} \kappa ^* = \kappa L_c \sim \frac {\epsilon }{L_c} \,{\boldsymbol\cdot}\, L_c = \epsilon . \end{equation}

Thus, the key geometric quantities scale as

(2.30) \begin{equation} \theta \sim \epsilon , \quad \textrm {cos}\,\theta \approx 1, \quad \sin \theta \sim \epsilon , \quad \kappa ^* \sim \epsilon . \end{equation}

For the spatial derivative, we obtain

(2.31) \begin{equation} \frac {\partial f}{\partial s} = \frac {\partial (L_c f^*)}{\partial (L_c s^*)} = \frac {\partial f^*}{\partial s^*} \approx \frac {\partial f^*}{\partial x^*}. \end{equation}

Consequently, the slope term transforms to

(2.32) \begin{equation} h_n\frac {\partial f}{\partial s} = L_c h_n^*\frac {\partial f^*}{\partial s^*} \end{equation}

and the curvature term becomes

(2.33) \begin{equation} \sin \theta \kappa (s)\left (\frac {f h_n}{\textrm {cos}\,\theta } - \frac {h_n^2}{2}\right ) = L_c \sin \theta \kappa ^*(s^*)\left (\frac {f^* h_n^*}{\textrm {cos}\,\theta } - \frac {(h_n^*)^2}{2}\right ). \end{equation}

To quantify the relative importance of two terms, we examine their ratio by considering the abovementioned geometric approximations,

(2.34) \begin{equation} \frac {\sin \theta \kappa ^*(s^*)\left (\dfrac {f^* h_n^*}{\textrm {cos}\,\theta } - \dfrac {(h_n^*)^2}{2}\right )}{h_n^*\dfrac {\partial f^*}{\partial s^*}}\sim \dfrac {\epsilon \,{\boldsymbol\cdot}\, \epsilon \,{\boldsymbol\cdot}\, \left(f^* h_n^* - \dfrac {(h_n^*)^2}{2}\right )}{h_n^* \,{\boldsymbol\cdot}\, \epsilon } \sim \epsilon \left(f^* - \frac {h_n^*}{2} \right) \sim \epsilon f^*, \end{equation}

where the final approximation follows the thin-film approximation $h_n^* \ll 1$ .

This analysis reveals that when $\epsilon f^* \ll 1$ , which is satisfied for moderate surface heights ( $f^* \lesssim {\mathcal{O}}(1)$ ) and small slopes ( $\epsilon \ll 1$ ), the curvature term represents an ${\mathcal{O}}(\epsilon ^2)$ correction to the slope term ${\mathcal{O}}(\epsilon )$ contribution. Therefore, the local evolution equation can be simplified by omitting curvature terms, resulting in the following global model:

(2.35) \begin{equation} \frac {\partial h_z}{\partial t} = \gamma \frac {\partial }{\partial x}\left [h_z\left (\frac {\partial h_z}{\partial x} + \frac {{\rm d}f}{{\rm d}x}\right )\right ], \end{equation}

where $\gamma = K\Delta \rho g /\phi \mu$ . This simplified equation captures the essential physics in global coordinates through two mechanisms: a nonlinear diffusion term $h_z \partial h_z/\partial x$ representing gravity-driven spreading, and a topographic forcing term $h_z df/dx$ describing the influence of surface gradients. This derivation approach can also be applied to the situation of viscous gravity currents, while we merely made an assumption about the limited slopes and neglected curvatures (Di & Huppert Reference Di and Huppert2024).

2.3.3. Finite-volume release

For finite-volume releases, the volume conservation for viscous fluid is expressed as

(2.36) \begin{equation} \phi \int _{0}^{x_f(t)} h_z \,{\rm d}x = \phi \int _{0}^{l} h_z \,{\rm d}x = q, \end{equation}

where $q$ represents the total fluid volume. Therein, considering that the front of the porous media gravity current denoted by $x_f(t)$ is unknown, we map the current domain $[0,x_f(t)]$ onto a fixed domain $[0,l]$ , with $x=l$ corresponding to the right end.

The right boundary condition is imposed as

(2.37) \begin{equation} h_z(l,t) = 0, \end{equation}

indicating that no fluid extends beyond $x=l$ .

The left boundary condition is given by

(2.38) \begin{equation} \left (\frac {\partial h_z}{\partial x}+\frac {{\rm d}f}{{\rm d}x}\right )\bigg |_{x=0} = 0, \end{equation}

which, derived from the governing equation, ensures a no-flux condition at the origin.

For the initial condition, the released fluid is commonly configured as a rectangular profile with height $h_0$ and length $x_0$ (where $x_0=q/(\phi h_0)$ ) at time $t=0$ , i.e.

(2.39) \begin{equation} h_z(x,0)=\begin{cases} h_0, & 0\leqslant x\leqslant x_0, \\[4pt] 0, & x \gt x_0. \end{cases} \end{equation}

2.4. Extension to axisymmetric propagation

For axisymmetric propagation over curvilinear surfaces, we can follow a derivation similar to that used for the two-dimensional propagation. The hydrostatic pressure and velocity fields in the local coordinate $(s,n)$ remain the same, but the mass conservation equation should be modified to account for axisymmetry. It becomes

(2.40) \begin{equation} \phi \frac {\partial h_n}{\partial t} + \frac {1}{r}\frac {\partial }{\partial s}\left (r\int _{f_n}^{h_n+f_n} u_s {\rm d}n\right ) = 0, \end{equation}

where the operator $1/r \,{\boldsymbol\cdot}\, \partial (r\,{\boldsymbol\cdot}\, )/\partial s$ reflects the variation of the annular cross-sectional area with radius.

Accordingly, the local evolution equation in global axisymmetric coordinates is

(2.41) \begin{equation} \begin{aligned} \phi \frac {\partial h_n}{\partial t} = \frac {K\Delta \rho g}{\mu }\frac {1}{r}\frac {\partial }{\partial s}\left [ r\,\textrm {cos}\,\theta \left (h_n\frac {\partial h_n}{\partial s} + h_n\frac {\partial f}{\partial s}\right ) + r\sin \theta \kappa (s)\left (\frac {f h_n}{\textrm {cos}\,\theta } - \frac {h_n^2}{2}\right )\right ], \end{aligned} \end{equation}

and the corresponding global equation becomes

(2.42) \begin{equation} \frac {\partial h_z}{\partial t} = \gamma \frac {1}{r}\frac {\partial }{\partial r}\left [rh_z\left (\frac {\partial h_z}{\partial r} + \frac {{\rm d}f}{{\rm d}r}\right )\right ], \end{equation}

where the term $h_z\partial h_z/\partial r$ captures radial spreading due to gravity, and $h_z {\rm d}f/{\rm d}r$ represents the leading-order topographic forcing.

We transform the domain from the current front $[0,r_f(t)]$ to the right end $r=l$ , and the conservation of volume requires

(2.43) \begin{equation} \phi \int _{0}^{r_f(t)}2\pi rh_z{\rm d}r = \phi \int _{0}^{l}2\pi rh_z{\rm d}r = q. \end{equation}

The boundary conditions include

(2.44) \begin{equation} h_z(l,t) = 0 \end{equation}

and

(2.45) \begin{equation} \left (\frac {\partial h_z}{\partial r} + \frac {{\rm d}f}{{\rm d}r}\right )\bigg |_{r=0} = 0. \end{equation}

The initial condition is taken to be a cylindrical volume of fluid with height $h_0$ and radius $r_0$ (where $r_0=\sqrt {q/\phi \pi h_0}$ ) at time $t=0$ ,

(2.46) \begin{equation} h_z(r,0)=\begin{cases} h_0, & 0\leqslant r\leqslant r_0, \\[2pt] 0, & r \gt r_0. \end{cases} \end{equation}

At this point, we have completed the theoretical framework for porous media gravity currents over rigid curvilinear surfaces, with the boundary and initial conditions specified for both two-dimensional and axisymmetric finite-volume releases.

3. Non-dimensionalisation

3.1. Dimensionless form

We consider sinusoidal surfaces as typical examples of the curvilinear topography commonly encountered in geological formations, which often represent gently undulating surfaces, defined by

(3.1) \begin{equation} f\left (x\right )=A\left [1-\cos(\lambda x)\right ] \ \text{(2D)},\quad f\left (r\right )=A\left [1-\cos(\lambda r)\right ]\ \text{(axisymmetric)}. \end{equation}

where the wavelength is given by $2\pi /\lambda$ and the total vertical variation is $2A$ .

To analyse the sinusoidal surfaces within our theoretical framework, we non-dimensionalise the system using scalings

(3.2) \begin{equation} S_x = 1/\lambda ,\quad S_h = A,\quad S_t = (\gamma \lambda ^2 A)^{-1}, \end{equation}

which leads to the dimensionless variables $X = x/S_x$ , $H = h_z/S_h$ and $T = t/S_t$ .

To ensure that the global model remains valid, we verify that our scalings are consistent with the fundamental approximations. Taking the characteristic length as $L_c = 1/\lambda$ , we require, first, that the dimensionaless surface height, defined as

(3.3) \begin{equation} f^* = f(x)/L_c = f\lambda = \lambda A\left [1-\textrm {cos}(\lambda x)\right ], \end{equation}

remains moderate (i.e. $f^* \lesssim {\mathcal{O}}(1)$ ), ensuring that the surface variations are small compared with the wavelength.

Second, the surface slope should be small,

(3.4) \begin{equation} \epsilon = \left |{\rm d}f/{\rm d}x\right | = \left |\lambda A\sin (\lambda x)\right |\ll 1. \end{equation}

Third, the fluid thickness remains small throughout the flow evolution,

(3.5) \begin{equation} h^* = h_n(x)/L_c \approx h_z\lambda \ll 1, \end{equation}

thereby justifying the thin-film approximation. When constraints are satisfied, the curvature effects contribute only high-order corrections, which can be omitted in the global model. Under these scalings, the sinusoidal surfaces are rendered in dimensionless form as

(3.6) \begin{equation} F(X) = 1-\textrm {cos}(X) \ \text{(2D)},\quad F(R) = 1-\textrm {cos}(R) \ \text{(axisymmetric)}. \end{equation}

The model for two-dimensional propagation of finite-volume releases then becomes

(3.7) \begin{equation} \frac {\partial H}{\partial T} = \frac {\partial }{\partial X}\left [H\left (\frac {\partial H}{\partial X} + \frac {{\rm d}F}{{\rm d}X}\right )\right ], \end{equation}

subject to the volume constraint

(3.8) \begin{equation} \int _{0}^{L} H \,{\rm d}X = N, \end{equation}

and boundary conditions

(3.9) \begin{equation} \left .\left (\frac {\partial H}{\partial X}+\frac {{\rm d}F}{{\rm d}X}\right )\right |_{X=0} = 0 \ \text{(no flux condition)},\quad H(X_{F},T) = H(L,T) = 0, \end{equation}

and the initial condition

(3.10) \begin{equation} H(X,0) = \begin{cases} H_0, & 0 \leqslant X \leqslant N/H_0, \\[4pt] 0, & X \gt N/H_0. \end{cases} \end{equation}

Similarly, the model for axisymmetric propagation of finite-volume releases becomes

(3.11) \begin{equation} \frac {\partial H}{\partial T} - \frac {1}{R}\frac {\partial }{\partial R}\left [RH\left (\frac {\partial H}{\partial R} + \frac {{\rm d}F}{{\rm d}R}\right )\right ] = 0, \end{equation}

subject to the volume constraint

(3.12) \begin{equation} 2\pi \int _{0}^{L} R H \,{\rm d}R = N, \end{equation}

the boundary conditions

(3.13) \begin{equation} \left .\left (\frac {\partial H}{\partial R}+\frac {{\rm d}F}{{\rm d}R}\right )\right |_{R=0} = 0,\quad H(R_{F},T) = H(L,T) = 0, \end{equation}

and the initial condition

(3.14) \begin{equation} H(R,0) = \begin{cases} H_0, & 0 \leqslant R \leqslant \sqrt {N/(\pi H_0)}, \\[6pt] 0, & R \gt \sqrt {N/(\pi H_0)}. \end{cases} \end{equation}

3.2. Physical interpretation

The dimensionless governing equations reveal two dimensionless parameters: the volume ratio $N$ and the initial height $H_0$ . Therein, the volume ratio $N$ quantifies the relative magnitude of the released fluid volume compared with the characteristic volume associated with one wavelength of the sinusoidal surface.

For the two-dimensional case, it is defined as

(3.15) \begin{equation} N = q/\phi A_s, \quad \text{with} \quad A_s = \int _{0}^{2\pi /\lambda } A\left [1-\textrm {cos}(\lambda x)\right ]\,{\rm d}x = 2\pi A/\lambda , \end{equation}

where $A_s$ represents the characteristic area of the topography.

For the axisymmetric case, we define

(3.16) \begin{equation} N = 4\pi ^3 q/\phi V_s, \quad \text{with} \quad V_s = 2\pi \int _{0}^{2\pi /\lambda } r A\left [1-\textrm {cos}(\lambda r)\right ]\,{\rm d}r = 4\pi ^3A/\lambda ^2, \end{equation}

where $V_s$ is the three-dimensional characteristic volume of the topography.

Physically, the dimensionless parameter $N$ governs the overall flow regime. When $N \ll 1$ with $H_0 \sim {\mathcal{O}}(1)$ , the fluid volume is too small to overcome the topography, leading to local topographic trapping that is largely influenced by the topography. As $N$ increases to order unity ( $N \sim 1$ ), a rich interplay between gravity-driven spreading and topographic forcing occurs. For even larger $N$ , the flow dynamics gradually resemble those over a flat surface.

For dimensionless parameter $H_0$ , a small value indicates that the initial release spreads over multiple wavelengths, whereas $H_0 \sim {\mathcal{O}}(1)$ suggests a more concentrated initial release that dominates the early-time dynamics. However, as gravity redistributes the fluid, the current thickness eventually decreases and the flow evolves into a thin-film regime.

4. Computational fluid dynamics approach

In this section, we present a more sophisticated CFD model for numerical validation of low-dimensional models. We first present two universal CFD models for two-phase flow in non-porous and porous media. By combining these two models, we develop a CFD model for macroscopic sharp interface flow in porous media. This CFD model is then used to simulate porous media gravity currents to provide validation of low-dimensional models.

4.1. The volume of fluid model

We begin with the volume of fluid model for classical two-phase flows with a sharp interface. Its incompressible version, as proposed by Hirt & Nichols (Reference Hirt and Nichols1981), is

(4.1) \begin{equation} {\boldsymbol\nabla} \,{\boldsymbol\cdot}\, \textbf{u}=0, \end{equation}

(4.2) \begin{equation} \frac {\partial \left (\rho \textbf{u}\right )}{\partial t}+{\boldsymbol\nabla} \,{\boldsymbol\cdot}\, \left (\rho \textbf{u}\textbf{u}\right )=-{\boldsymbol\nabla} p+\mu {\boldsymbol\nabla} ^{2}\textbf{u}+\rho \boldsymbol{\textbf{g}} +\sigma {\boldsymbol\nabla} \,{\boldsymbol\cdot}\, \left (\frac {{\boldsymbol\nabla} \alpha }{\left |{\boldsymbol\nabla} \alpha \right |}\right ){\boldsymbol\nabla} \alpha , \end{equation}

(4.3) \begin{equation} \frac {\partial \alpha }{\partial t}+{\boldsymbol\nabla} \,{\boldsymbol\cdot}\, \left (\boldsymbol\{\textbf{u} \alpha \right )+{\boldsymbol\nabla} \,{\boldsymbol\cdot}\, \left [ \alpha \left (1-\alpha \right )\left (C\left |{\textbf{u}}\right |\frac {{\boldsymbol\nabla} \alpha }{\left |{\boldsymbol\nabla} \alpha \right |}\right )\right ]=0. \end{equation}

This model adopts a single-field formulation where both phases share common velocity ( $\textbf{u}$ ) and pressure ( $p$ ) fields. The spatial distribution of the two fluids is characterised by a volume fraction $\alpha$ , where $\alpha = 1$ denotes regions fully occupied by fluid 1 and $\alpha = 0$ indicates the presence of only fluid 2. At the interface between fluids, surface tension $\sigma$ induces capillary effects, while the local density and viscosity are determined through linear interpolation, $\rho = \rho _{1} \alpha + \rho _{2} (1-\alpha )$ and $\mu = \mu _{1} \alpha + \mu _{2} (1-\alpha )$ . To maintain a numerically sharp interface between the fluids, the model incorporates a compression factor $C$ .

4.2. Macroscopic generalised two-phase flow model

Two-phase flows in porous media exhibit fundamentally different behaviour across multiple length scales. At the pore scale, the volume of fluid model, for instance, may be used to capture interfacial dynamics and local flow patterns. However, practical engineering applications, which typically span larger spatial domains, necessitate macroscopic models derived through upscaling. For instance, the volume-averaging method (Quintard & Whitaker Reference Quintard and Whitaker1988; Whitaker Reference Whitaker1998) establishes a framework for upscaling by volume-averaging pore-scale phenomena over representative elementary volumes.

Many different models have been proposed, starting from the classical generalised Darcy’s law to more complex models incorporating more dynamic and phase interaction effects (Davit & Quintard Reference Davit and Quintard2019). In this study, the assumed physical constraints lead to a two-phase flow regime characterised by a macroscopic sharp interface (i.e. a very small capillary fringe compared with the relevant characteristic lengths). As a consequence, the development of the model may start with a classical formulation of the macro-scale momentum equations incorporating a classical generalised Darcy’s law, since the details of the two-phase flow model will play a negligible role outside the macro-scale interface. Various arrangements of the variable field (saturation, pressures, velocities) may be used, e.g.: (i) a two-field formulation that resolves separate Darcy velocities for each phase (Horgue et al. Reference Horgue, Soulaine, Franc, Guibert and Debenest2015); (ii) a single-field formulation based on the total filtration velocity (Carrillo, Bourg & Soulaine Reference Carrillo, Bourg and Soulaine2020). The latter can be expressed as

(4.4) \begin{equation} {\boldsymbol\nabla} \,{\boldsymbol\cdot}\, \textbf{U}=0, \end{equation}

(4.5) \begin{align} \frac {\partial \left (\rho \textbf{U}\right )}{\partial t}+{\boldsymbol\nabla} \,{\boldsymbol\cdot}\, \left (\rho \textbf{U}\textbf{U}\right )&=-{\boldsymbol\nabla} P+\mu {\boldsymbol\nabla} ^{2}\textbf{U}+\rho \boldsymbol{\textbf{g}}-\phi k_{0}^{-1}\left (\frac {k_{r1}}{\mu _{1}}+\frac {k_{r2}}{\mu _{2}}\right )^{-1}\textbf{U}\nonumber\\[6pt] &\quad +\left [k_{0}^{-1}\left (\frac {\mu _{1}}{k_{r1}}+\frac {\mu _{2}}{k_{r2}}\right )^{-1}\left (M_{1}\alpha _{2}-M_{2}\alpha _{1}\right )\frac {\partial P_{c}}{\partial \alpha _{1}}-P_{c}\right ]{\boldsymbol\nabla} S_{1}, \end{align}

(4.6) \begin{equation} \quad \frac {\partial S_{1}}{\partial t}+{\boldsymbol\nabla} \,{\boldsymbol\cdot}\, \left (\textbf{U}S_{1}\right ) +{\boldsymbol\nabla} \,{\boldsymbol\cdot}\, (S_{1}S_{2}\textbf{U}_{r})=0, \end{equation}

(4.7) \begin{align} \textbf{U}_{r} & =-\phi ^{-1}\left [\left (\frac {M_{1}}{S_{1}}-\frac {M_{2}}{S_{2}}\right ){\boldsymbol\nabla} P+\left (\frac {\rho _{1}M_{1}}{S_{1}}-\frac {\rho _{2}M_{2}}{S_{2}}\right )\textbf{g}\right . \nonumber\\[6pt] & \left . \quad +\left (\frac {M_{1}S_{2}}{S_{1}}+\frac {M_{2}S_{1}}{S_{2}}\right ){\boldsymbol\nabla} P_{c}-\left (\frac {M_{1}}{S_{1}}-\frac {M_{2}}{S_{2}}\right )P_{c}{\boldsymbol\nabla} S_{1}\right ]. \end{align}

This single-field formulation characterises macroscopic flow by the intrinsic averaged velocity $\textbf{U}$ and pressure $P$ , while the fluid distribution is represented by saturations $S_i$ ( $i=1,2$ ) within the porous medium. The model incorporates both the intrinsic permeability $k_0$ and relative permeability $k_{ri}$ . The fluid mobility $M_i$ is defined as the ratio of these permeabilities, $M_i = k_0k_{ri}/\mu _{i}$ . The model requires closure relations for the macroscopic capillary pressure $P_c(S_1)$ and relative permeability functions $k_{ri}(S_1)$ .

4.3. Macroscopic sharp-interface flow model

Given the negligible macroscopic capillary effects and gravity segregation effects, the macroscopic interface between fluids in porous media gravity currents can be treated as sharp. Drawing parallels with the volume of fluid model, we modify the single-field Darcy formulation through strategic reconstruction of the macroscopic capillary term and relative velocity in the saturation equation. This yields the following sharp-interface model for macroscopic flows in porous media:

(4.8) \begin{equation} {\boldsymbol\nabla} \,{\boldsymbol\cdot}\, \textbf{U}=0, \end{equation}

(4.9) \begin{equation} \begin{aligned} \frac {\partial \left (\rho \textbf{U}\right )}{\partial t}+{\boldsymbol\nabla} \,{\boldsymbol\cdot}\, \left (\rho \textbf{U}\textbf{U}\right )=-{\boldsymbol\nabla} P +\mu {\boldsymbol\nabla} ^{2}\textbf{U} +\rho \boldsymbol{\textbf{g}} -\phi \mu K^{-1}\textbf{U}+\sigma _{eff}{\boldsymbol\nabla} \,{\boldsymbol\cdot}\, \left (\frac {{\boldsymbol\nabla} S}{\left |{\boldsymbol\nabla} S\right |}\right ){\boldsymbol\nabla} S, \end{aligned} \end{equation}

(4.10) \begin{equation} \frac {\partial S}{\partial t}+{\boldsymbol\nabla} \,{\boldsymbol\cdot}\, \left (\textbf{U}S\right )+{\boldsymbol\nabla} \,{\boldsymbol\cdot}\, \left [\phi S\left (1-S\right )\left (C\left |\textbf{U}\right |\frac {{\boldsymbol\nabla} S}{\left |{\boldsymbol\nabla} S\right |}\right )\right ]=0. \end{equation}

This formulation offers a significant advantage in that it eliminates the need for empirical closure relations describing macroscopic capillary effects and relative permeability, making it particularly well suited for macroscopic sharp-interface flows in porous media. The effective surface tension at the macro scale is denoted by $\sigma _{eff}$ , which is neglected in the assumption of porous media gravity currents, and $K$ is the intrinsic permeability scalar in this study. Because porous media gravity currents characteristically propagate at low Reynolds numbers, inertial corrections are neglected. We implement this model within the open-source CFD software OpenFOAM $^{\circledR }$ for numerical simulation of porous media gravity currents.

5. Numerical validation

We perform CFD simulations using the macroscopic sharp-interface flow model and compare the results with predictions from low-dimensional models. For the physical parameters of the porous medium, we adopt values from the literature (Acton, Huppert & Worster Reference Acton, Huppert and Worster2001; Lyle et al. Reference Lyle, Huppert, Hallworth, Bickle and Chadwick2005), where an homogeneous porous medium was generated in a bead-filled tank, yielding a porosity of 0.37 and a permeability of $6.8\times 10^{-9}$ m $^{2}$ . Glycerol is used as the working fluid, with a dynamic viscosity of $\mu = 1.5$ Pa $\,{\boldsymbol\cdot}\,$ s and density $\rho = 1.26\times 10^3$ kg m $^{-3}$ .

Our numerical domain is designed to be 0.15 m in the horizontal direction and 0.05 m vertically. For the bottom sinusoidal surfaces, we employ the same parameters as used in our previous study on viscous gravity currents (Di & Huppert Reference Di and Huppert2024), with a wavelength $2\pi /\lambda = 0.05$ m and amplitude $A = 0.002$ m. This configuration results in $A\lambda = 0.25 \lt 1$ , which satisfies the small-slope and limited surface amplitude conditions required by the low-dimensional models. The laboratory-scale settings ensure that the flow remains in the creeping flow regime where Darcy’s law applies, thereby justifying the use of the CFD model for numerical validation.

We impose atmospheric conditions at the top and right boundaries, a symmetry condition at the left boundary and a no-slip condition at the bottom boundary. Local mesh refinement near the bottom is employed to the wall regions. The initial condition consists of a glycerol volume released from the origin, with its upper surface conforming to the bottom profile to ensure a consistent initial rectangular configuration relative to the surface. Here, we examine three initial configurations with a release length of $x_0 (r_0)= 0.05$ m (i.e. $X_0 (R_0)= 1$ ) and heights of $h_0 = 0.01$ , 0.02 and 0.03 m (i.e. $H_0 = 5$ , 10 and 15).

Numerical accuracy is ensured via grid independence studies, where simulations were conducted using four mesh resolutions (120 000, 240 000, 360 000 and 480 000 cells). Convergence is assessed by comparing the current profile and the front position, using MATLAB $^{\circledR }$ image processing for post-processing. An intermediate resolution of 240 000 cells was found to offer an optimal balance between accuracy and computational efficiency. All simulations were performed on a workstation equipped with Intel Xeon E5-2678 processors using 12 CPU cores, with typical runs requiring approximately 30 hours to simulate 1000 s of physical time.

Figures 2, 3 and 4 illustrate the evolution of two-dimensional porous media gravity currents for three cases with increasing initial heights, while keeping the initial release length constant. This corresponds to increasing the released fluid volume, which not only drives the propagation of the current further, but also enhances the gravitational potential energy available during the initial deformation.

Figure 2.

Numerical comparison of two approaches to current profiles for a two-dimensional porous media gravity current over a sinusoidal surface (Case 1). Dotted yellow line, numerical solution.

Figure 3.

Numerical comparison of two approaches to current profiles for a two-dimensional porous media gravity current over a sinusoidal surface (Case 2). Dotted yellow line, numerical solution.

Figure 4.

Numerical comparison of two approaches to current profiles for a two-dimensional porous media gravity current over a sinusoidal surface (Case 3). Dotted yellow line, numerical solution.

In Case 1, with the lowest initial height, the CFD results and the low-dimensional model predictions agree very well over time. The low initial height ensures that the initial configuration closely satisfies the thin-film approximation. Consequently, the fluid quickly reaches a regime in which the low-dimensional model is valid. Moreover, owing to the finite volume, once the fluid fills the first two troughs of the topography, the interface becomes nearly flat relative to the global coordinate, which limits further propagation.

In Case 2, the increased initial height leads to noticeable discrepancies at early times (e.g. at $t = 50$ s). Here, the CFD model still exhibits an influence of the initial fluid shape, whereas the low-dimensional model rapidly develops a monotonically decreasing interface that propagates further to the right. As the flow continues to evolve, the predictions of both models gradually converge.

In Case 3, where the initial height is further increased, the differences between the two models become even more pronounced at early times. For instance, at $t = 50$ s, the CFD results clearly reflect the initial configuration, and even at $t = 100$ s, the influence of the initial shape persists. Only after a sufficiently long time does the flow evolve towards the thin-film regime, at which point, the predictions of the low-dimensional model begin to match those of the CFD simulations.

Overall, these examples demonstrate that although a higher initial fluid height intensifies the dynamic influence of the initial conditions, as time increases, the released fluid eventually leads to a state where the thin-film approximation is well satisfied. The results of axisymmetric propagation of the same initial released dimension are displayed in Appendix A.

Figure 5 quantitatively compares the evolution of the current front as predicted by low-dimensional models and captured by CFD simulations. For the smallest release volume, both approaches show excellent agreement. However, as the initial fluid height increases, discrepancies begin to emerge, with the low-dimensional model slightly overestimating the propagation rate. This divergence is primarily due to significant vertical deformation in the early-stage dynamics leads to enhanced viscous dissipation, but cannot be fully captured by the horizontal Darcy’s law in low-dimensional models, whereas the Navier–Stokes formulation in the CFD simulations resolves these effects more accurately.

It should be noted that this study focuses on finite-volume releases. In this regime, both the low-dimensional and CFD models predict that the fluid ultimately reaches an equilibrium state, consistent with energy minimisation principles, in which the finite volume is confined within the troughs of the sinusoidal surface. In contrast, under constant-flux injection, the continuous input of fluid sustains a dynamic, non-equilibrium state as the fluid persistently overcomes topographic variations. Therefore, it would be interesting to explore the applicability of the low-dimensional approach to this flow pattern.

Figure 5.

Numerical comparison of two approaches to dimensionless current fronts for a porous media gravity current over a sinusoidal surface in both two-dimensional and axisymmetric propagations.

7. Relevance to CO $_{\textbf{2}}$ geological sequestration

Although our numerical simulations were performed at laboratory scales, the fundamental dynamics captured by low-dimensional models can extend naturally to engineering-scale scenarios, such as CO $_2$ geological sequestration. In practice, CO $_2$ is injected under constant-flux conditions. Once injection ceases, the system behaves as a finite-volume release. Under such conditions, the flow is governed by the same topographic trapping mechanisms that the finite-volume framework describes. Moreover, even over kilometre-scale distances, the balance between viscous forces and gravity, central to the thin-film approximation, remains valid, ensuring that the key physical processes are effectively captured. This provides a foundation for applying low-dimensional models to study supercritical CO $_2$ flow over wavy cap rocks in realistic geological formations. Such conditions are discussed in detail by Huppert & Neufeld (Reference Huppert and Neufeld2014).

For quantitative analysis, we consider supercritical CO $_2$ flow with physical properties from Pegler et al. (Reference Pegler, Huppert and Neufeld2013): a density of $\rho = 700$ kg m $^{-3}$ , a dynamic viscosity of $\mu = 6 \times 10^{-5}$ Pa $\,{\boldsymbol\cdot}\,$ s and an ambient brine density of $\rho _0 = 1000$ kg m⁻ $^3$ , within a geological formation characterised by a porosity $\phi = 0.4$ and a permeability $K = 1 \times 10^{-12}$ m $^2$ . Following Woods (Reference Woods2015), we represent cap rocks, distorted by compaction and tectonic stresses, as sinusoidal surfaces with wavelengths $2\pi /\lambda = 10, 15$ m and amplitudes $A = 0.5, 1$ m. To assess the applicability of low-dimensional models, we compute the constraint $A\lambda$ . For a 10 m wavelength, $A\lambda$ is 0.31 and 0.63 for amplitudes of 0.5 m and 1 m, respectively; for a 15 m wavelength, $A\lambda$ takes on lower values of 0.21 and 0.42. These results ensure that low-dimensional models remain valid for predicting the supercritical CO $_2$ propagation. For supercritical CO $_2$ volume, we set the initial released dimensions at 10 m $\times$ 20 m, corresponding to an area of 200 m $^2$ in a two-dimensional configuration and a volume of 283.2 m $^3$ in an axisymmetric configuration.

Figure 6 illustrates the propagation of supercritical CO $_2$ over both flat and wavy cap rocks under finite-volume release conditions. Unlike the continuous, monotonically increasing propagation observed over a flat surface, the flow over wavy cap rocks exhibits a step-like pattern. This behaviour arises because the fluid gradually fills the troughs as it propagates to the right, then rapidly flows downhill upon reaching the crests of the sinusoidal surface, and subsequently accumulates slowly once more in the troughs. Over time, the step region expands, indicating that the current front decelerates and requires increasingly longer periods to fully fill the troughs until eventually all the fluid becomes trapped.

Figure 6.

Comparison of supercritical CO $_{2}$ over horizontal flat and wavy cap rocks in both two-dimensional and axisymmetric propagations.

A comparative analysis of different cap rock geometries reveals that a wavy cap rock with a smaller amplitude (0.5 m) and a longer wavelength (15 m) permits greater propagation, whereas cap rocks with a larger amplitude (1 m) and a shorter wavelength (10 m) more effectively restrict the flow. These findings suggest an inverse relationship between the surface characteristic volume and the dimensionless release volume, $N$ ; a larger characteristic surface volume corresponds to a smaller $N$ , thereby reducing propagation. Notably, for the 15 m wavelength, the propagation eventually stagnates, indicating complete topographic trapping within the troughs, in contrast to the unbounded spreading observed over flat surfaces.

Further, our results indicate that in two-dimensional propagation, the current is more sensitive to variations in amplitude than to wavelength, suggesting the dominant influence of amplitude. In axisymmetric propagation, however, while amplitude remains important, the effect of wavelength becomes more pronounced, i.e. longer wavelengths facilitate greater radial propagation. This distinction arises because the dimensionless volume $N$ scales linearly with wavelength in two-dimensional flow (3.15), but quadratically in axisymmetric configurations (3.16). In radially spreading currents, annular troughs with increasing circumference impose progressively larger volumetric demands with radius, unlike the constant cross-sectional area in two-dimensional propagation, thereby accentuating the wavelength sensitivity in axisymmetric propagation.

In the context of CO $_2$ geological sequestration, where millions of cubic metres of supercritical CO $_2$ are injected for long-term storage over thousands of years and may migrate several kilometres within subsurface strata, two primary mechanisms enhance storage security: capillary trapping within the porous medium and dissolution into the formation brine. In this light, wavy cap rocks not only serve as local structural traps by confining supercritical CO $_2$ within their troughs, but also extend the residence time of supercritical CO $_2$ , thereby enhancing the effectiveness of these stabilising processes. This suggests that naturally occurring wavy cap rock geometries can improve long-term storage security in practical CO $_2$ sequestration.

8. Future outlook

Our theoretical framework for porous media gravity currents over curvilinear surfaces extends the seminal contributions of Huppert & Woods (Reference Huppert and Woods1995) and Lyle et al. (Reference Lyle, Huppert, Hallworth, Bickle and Chadwick2005). By capturing essential topographic effects, it lays a foundation for understanding the complex physics of subsurface flows, with particular relevance to CO $_{2}$ sequestration. However, in practical applications, there are other factors that need to be studied.

One critical challenge is the migration of supercritical CO $_2$ through fractures or low-permeability zones in cap rocks. Although previous studies (Neufeld, Vella & Huppert Reference Neufeld, Vella and Huppert2009; Neufeld et al. Reference Neufeld, Vella, Huppert and Lister2011) have explored fissure leakage over flat surfaces, the coupling between leakage pathways and wavy cap rocks remains poorly understood. In such settings, the local geometry of the sealing formation can alter pressure distributions and divert flow towards or away from leakage points. A deeper understanding of these interactions is vital for developing robust risk assessment strategies.

Confinement effects introduce another complexity when combined with curvilinear boundaries (Pegler, Huppert & Neufeld Reference Pegler, Huppert and Neufeld2014; Zheng et al. Reference Zheng, Guo, Christov, Celia and Stone2015). In contrast to flat confined systems, the curvilinear boundaries generate variable effective geometries. This alteration influences both the pressure gradients that drive horizontal flow and the maximum hydrostatic head that limits propagation. Future studies should quantify how the geometric modifications affect flow dynamics under confinement.

With permeability varying between horizontal and vertical directions (Hinton & Woods Reference Hinton and Woods2018; Benham, Bickle & Neufeld Reference Benham, Bickle and Neufeld2021), the alignment of preferential flow paths with the local topography can change spatially. This directional dependence makes the effective permeability a function of both position and orientation, potentially amplifying or mitigating the influence of surface slope and curvature on pressure distributions. Addressing this anisotropy is essential for further accurately predicting flow patterns.

Moreover, while low-dimensional models assume a sharp, macroscopic interface, realistic geological settings often involve multiphase flow with complex interfacial dynamics (Golding et al. Reference Golding, Neufeld, Hesse and Huppert2011; Golding, Huppert & Neufeld Reference Golding, Huppert and Neufeld2017). Incorporating factors such as relative permeability and capillary pressure into low-dimensional models with curvilinear boundaries is necessary to capture the interplay between gravity and capillary forces, as well as processes like mixing, dissolution and instabilities that are critical for long-term CO $_{2}$ trapping.

In summary, future research should focus on integrating these coupled mechanisms into low-dimensional models that remain computationally tractable for large-scale geological applications. Investigating how cap rock geometries interact with features such as leakage pathways, confinement and anisotropy will be key to advancing our understanding and predictive capability in CO $_{2}$ sequestration scenarios.