# Low-dimensional modelling of viscous fluids in porous media over curvilinear surfaces

Xiaoning Di, Yi Huang, Michel Quintard, Herbert Huppert

► To cite this version:

Xiaoning Di, Yi Huang, Michel Quintard, Herbert Huppert. Low-dimensional modelling of viscous fluids in porous media over curvilinear surfaces. Journal of Fluid Mechanics, 2025, 1016, pp.A16. 10.1017/jfm.2025.10168. hal-05221362

HAL Id: hal-05221362

https://hal.science/hal-05221362v1

Submitted on 25 Aug 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Low-dimensional modelling of viscous fluids in porous media over curvilinear surfaces

Xiaoning Di<fnref n="1" />,<fnref n="2" />,<fnref n="3" /> ID, Yi Huang<fnref n="1" />, Michel Quintard<fnref n="2" /> ID and Herbert E. Huppert<fnref n="3" /> ID

**Corresponding authors:** Herbert E. Huppert, heh1@cam.ac.uk; Xiaoning Di, xiaoning0802@gmail.com

(Received 1 October 2024; revised 18 February 2025; accepted 20 March 2025)

We present a theoretical framework for porous media gravity currents propagating over rigid curvilinear surfaces. By reducing the flow dynamics to low-dimensional models applicable on surfaces where curvature effects are negligible, we demonstrate that, for finite-volume releases, the flow behaviour in both two-dimensional and axisymmetric configurations is primarily governed by the ratio of the released viscous fluid volume to the characteristic volume of the curvilinear surface. Our theoretical predictions are validated using computational fluid dynamics simulations based on a sharp-interface model for macroscopic flow in porous media. In the context of carbon dioxide geological sequestration, our findings suggest that wavy cap rock geometries can enhance trapping capacity compared with traditional flat-surface assumptions, highlighting the importance of incorporating realistic topographic features into subsurface flow models.

**Key words:** porous media, gravity currents, topographic effects

## 1. Introduction

The dynamics of porous media gravity currents is a fundamental issue in fluid mechanics, with applications ranging from groundwater hydrology to carbon dioxide (CO₂) geological sequestration. Despite their inherent complexity, these flows can be effectively described by low-dimensional models, which offer both predictive capacity and valuable insights into their spatial and temporal evolution (Philip 1970; Huppert & Neufeld 2014).

<footnote n="1">¹School of Naval Architecture and Ocean Engineering, Dalian University of Technology, Dalian 116024, PR China</footnote>

<footnote n="2">²Institut de Mécanique des Fluides de Toulouse (IMFT), Université de Toulouse, CNRS, Toulouse 31400, France</footnote>

<footnote n="3">³Institute of Theoretical Geophysics, King's College, University of Cambridge, Cambridge CB2 1ST, UK</footnote>
A theoretical framework for gravity currents propagating through porous media was established by Huppert & Woods (1995), who developed low-dimensional models for two-dimensional propagation over both horizontal and inclined flat surfaces, and proposed self-similar solutions. This framework was subsequently extended to axisymmetric configurations by Lyle et al. (2005) and further generalised to inclined situations by Vella & Huppert (2006). Later studies enriched the theoretical framework by incorporating non-Newtonian effects for power-law fluids (Di Federico, Archetti & Longo 2012), investigated experimentally by Longo et al. (2013), thereby underscoring the importance of rheological characterisation.

Despite these theoretical advances, a major limitation of existing models is their common assumption of flat, horizontal surfaces, which fails to represent the complex topography found in practical situations, such as in CO₂ geological sequestration, where cap rocks often exhibit irregular, wavy geometries (Woods 2015). While Pegler, Huppert & Neufeld (2013) extended the analysis to porous media gravity currents over power-law surfaces, these monotonic curvilinear approximations do not capture the undulating nature of cap rocks. Therefore, the disconnect between idealised models and realistic geological formations constitutes a critical gap in our understanding of porous media gravity currents over complex topography.

To bridge this gap, we develop a theoretical framework for porous media gravity currents propagating over non-monotonic curvilinear surfaces. Building on the principles of hydrostatic pressure distribution and Darcy's law, we first establish a model for the local evolution of viscous fluids over curvilinear surfaces. Then, under the assumptions of small slopes and negligible curvature effects, we transform the local model into a global model that describes the overall evolution of the current. Through non-dimensionalisation, we identify a key dimensionless parameter that governs the flow dynamics in finite-volume releases, offering a straightforward physical interpretation. Finally, we validate the global model using computational fluid dynamics (CFD) simulations.

This paper is organised as follows. In §2, we develop low-dimensional models for porous media gravity currents propagating over curvilinear surfaces. Section 3 employs non-dimensionalisation to identify the key parameters that govern topographic effects. In §4, we present our CFD approach for capturing macroscopic sharp interfaces of porous media gravity currents. Section 5 provides numerical validations of the theoretical predictions from our low-dimensional models. Section 6 offers further discussions of the proposed model, while §7 explores the application to CO₂ geological sequestration. Section 8 outlines the future outlook on this fundamental issue, and finally, §9 concludes the study.

## 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
<figure><img src="image_3.png" /><figcaption>Figure 1. Schematic diagram of a gravity current in an homogeneous porous medium over an impermeable and rigid curvilinear surface.</figcaption></figure>

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.

## 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

$$P = p_0 + \Delta\rho g \cos \theta (h_n + f_n - n), \quad (2.1)$$

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)/\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

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

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

$$\frac{\partial f_n}{\partial s} = \frac{1}{\cos\theta} \frac{\partial f}{\partial s} + f(s) \frac{\partial}{\partial s} \left( \frac{1}{\cos\theta} \right) = \frac{1}{\cos\theta} \frac{\partial f}{\partial s} + f(s) \left( \frac{\sin\theta}{\cos^2\theta} \right) \frac{\partial\theta}{\partial s}. \quad (2.3)$$

Introducing the local surface curvature defined as

$$\kappa(s) = \frac{\partial\theta}{\partial s}, \quad (2.4)$$

this expression becomes

$$\frac{\partial f_n}{\partial s} = \frac{1}{\cos\theta} \frac{\partial f}{\partial s} + \frac{f(s) \sin\theta \kappa(s)}{\cos^2\theta}. \quad (2.5)$$

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

$$\frac{\partial}{\partial s}(\cos\theta) = -\sin\theta \frac{\partial\theta}{\partial s} = -\sin\theta \kappa(s). \quad (2.6)$$
Thus, the final expression for the tangential pressure gradient is

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

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

$$ \frac{\partial P}{\partial n} = -\Delta \rho g \cos \theta. \quad (2.8) $$

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,

$$ u_s = -\frac{K}{\mu} \frac{\partial P}{\partial s}, \quad (2.9) $$

we obtain the local velocity field

$$ u_s = -\frac{K \Delta \rho g}{\mu} \left[ \cos \theta \left( \frac{\partial h_n}{\partial s} + \frac{1}{\cos \theta} \frac{\partial f}{\partial s} + \frac{f \sin \theta \kappa(s)}{\cos^2 \theta} \right) - \sin \theta \kappa(s) \left( h_n + \frac{f}{\cos \theta} - n \right) \right]. \quad (2.10) $$

That is,

$$ u_s = -\frac{K \Delta \rho g}{\mu} \left[ \underbrace{\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)}{\cos \theta} - \sin \theta \kappa(s) \left( h_n + \frac{f}{\cos \theta} - n \right) \right)}_{\text{Curvature-driven terms}} \right]. \quad (2.11) $$

### 2.2.2. Local evolution

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

$$ \phi \frac{\partial h_n}{\partial t} + \frac{\partial}{\partial s} \int_{f_n}^{h_n+f_n} u_s \, dn = 0. \quad (2.12) $$

Evaluating the integral term, we obtain

$$ \int_{f_n}^{h_n+f_n} \cos \theta \frac{\partial h_n}{\partial s} \, dn = h_n \cos \theta \frac{\partial h_n}{\partial s}, \quad (2.13) $$

$$ \int_{f_n}^{h_n+f_n} \frac{\partial f}{\partial s} \, dn = h_n \frac{\partial f}{\partial s}, \quad (2.14) $$

$$ \int_{f_n}^{h_n+f_n} \left[ \frac{f \sin \theta \kappa(s)}{\cos \theta} - \sin \theta \kappa(s)(h_n + f_n - n) \right] dn = \sin \theta \kappa(s) \left( \frac{fh_n}{\cos \theta} - \frac{h_n^2}{2} \right). \quad (2.15) $$
Thus, the local evolution equation is given by

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

This local evolution equation reveals that the fluid dynamics is governed by three distinct physical mechanisms. First, the term $h_n \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) (fh_n / \cos \theta - h_n^2 / 2)$, capture more effects arising from the interplay between fluid and surface geometry. Therein, the component $fh_n / \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 > 0$) and to assist it in concave regions (where $\kappa < 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

$$ h_z = h_n \cos \theta + \beta \kappa_z h_n^2 + \mathcal{O}(h_n^3), \quad (2.17) $$

where the first term, $h_n \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

$$ \kappa_z = \frac{\partial^2 h_z / \partial x^2}{\left[1 + \left(\partial h_z / \partial x\right)^2\right]^{3/2}}. \quad (2.18) $$

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

$$ \kappa_z \approx \frac{\partial^2 h_z}{\partial x^2} \sim \frac{h_n}{L_c^2} \ll 1. \quad (2.19) $$

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

$$ h_z = h_n \cos \theta, \quad (2.20) $$

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

$$ \frac{\partial h_n}{\partial t} = \frac{1}{\cos \theta} \frac{\partial h_z}{\partial t}. \qquad (2.21) $$

### 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:

$$ 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. \qquad (2.22) $$

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,

$$ \epsilon = |\mathrm{d}f^*/\mathrm{d}x^*| \ll 1, \qquad (2.23) $$

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

$$ \tan \theta \approx \theta \approx |\mathrm{d}f^*/\mathrm{d}x^*|, \qquad (2.24) $$

which implies that

$$ \theta \sim \epsilon, \qquad (2.25) $$

and we further obtain

$$ \cos \theta \approx 1 - \frac{\theta^2}{2} \approx 1, \quad \sin \theta \approx \theta \sim \epsilon. \qquad (2.26) $$

Because the curvature is defined as

$$ \kappa = \frac{d\theta}{ds}, \qquad (2.27) $$

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

$$ \kappa \sim \frac{\Delta\theta}{L_c} \sim \frac{\epsilon}{L_c}, \qquad (2.28) $$

$$ \kappa^* = \kappa L_c \sim \frac{\epsilon}{L_c} \cdot L_c = \epsilon. \qquad (2.29) $$

Thus, the key geometric quantities scale as

$$ \theta \sim \epsilon, \quad \cos \theta \approx 1, \quad \sin \theta \sim \epsilon, \quad \kappa^* \sim \epsilon. \qquad (2.30) $$

For the spatial derivative, we obtain

$$ \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^*}. \qquad (2.31) $$

Consequently, the slope term transforms to

$$ h_n \frac{\partial f}{\partial s} = L_c h_n^* \frac{\partial f^*}{\partial s^*} \qquad (2.32) $$

and the curvature term becomes

$$ \sin \theta \kappa(s) \left( \frac{fh_n}{\cos \theta} - \frac{h_n^2}{2} \right) = L_c \sin \theta \kappa^*(s^*) \left( \frac{f^*h_n^*}{\cos \theta} - \frac{(h_n^*)^2}{2} \right). \qquad (2.33) $$
To quantify the relative importance of two terms, we examine their ratio by considering the abovementioned geometric approximations,

$$ \frac{\sin \theta \kappa^*(s^*) \left( \frac{f^* h_n^*}{\cos \theta} - \frac{(h_n^*)^2}{2} \right)}{h_n^* \frac{\partial f^*}{\partial s^*}} \sim \frac{\epsilon \cdot \epsilon \cdot \left( f^* h_n^* - \frac{(h_n^*)^2}{2} \right)}{h_n^* \cdot \epsilon} \sim \epsilon \left( f^* - \frac{h_n^*}{2} \right) \sim \epsilon f^*, \quad (2.34) $$

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:

$$ \frac{\partial h_z}{\partial t} = \gamma \frac{\partial}{\partial x} \left[ h_z \left( \frac{\partial h_z}{\partial x} + \frac{df}{dx} \right) \right], \quad (2.35) $$

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_t 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 2024).

### 2.3.3. Finite-volume release

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

$$ \phi \int_{0}^{x_f(t)} h_z \, dx = \phi \int_{0}^{l} h_z \, dx = q, \quad (2.36) $$

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

$$ h_z(l, t) = 0, \quad (2.37) $$

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

The left boundary condition is given by

$$ \left. \left( \frac{\partial h_z}{\partial x} + \frac{df}{dx} \right) \right|_{x=0} = 0, \quad (2.38) $$

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.

$$ h_z(x, 0) = \begin{cases} h_0, & 0 \le x \le x_0, \\ 0, & x > x_0. \end{cases} \quad (2.39) $$

## 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

$$ \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 \, dn \right) = 0, \qquad (2.40) $$

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

Accordingly, the local evolution equation in global axisymmetric coordinates is

$$ \phi \frac{\partial h_n}{\partial t} = \frac{K \Delta \rho g}{\mu} \frac{1}{r} \frac{\partial}{\partial s} \left[ r \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}{\cos \theta} - \frac{h_n^2}{2} \right) \right], \qquad (2.41) $$

and the corresponding global equation becomes

$$ \frac{\partial h_z}{\partial t} = \gamma \frac{1}{r} \frac{\partial}{\partial r} \left[ r h_z \left( \frac{\partial h_z}{\partial r} + \frac{\mathrm{d}f}{\mathrm{d}r} \right) \right], \qquad (2.42) $$

where the term $h_z \partial h_z / \partial r$ captures radial spreading due to gravity, and $h_z \mathrm{d}f / \mathrm{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

$$ \phi \int_{0}^{r_f(t)} 2\pi r h_z \, dr = \phi \int_{0}^{l} 2\pi r h_z \, dr = q. \qquad (2.43) $$

The boundary conditions include

$$ h_z(l, t) = 0 \qquad (2.44) $$

and

$$ \left. \left( \frac{\partial h_z}{\partial r} + \frac{\mathrm{d}f}{\mathrm{d}r} \right) \right|_{r=0} = 0. \qquad (2.45) $$

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$,

$$ h_z(r, 0) = \begin{cases} h_0, & 0 \le r \le r_0, \\ 0, & r > r_0. \end{cases} \qquad (2.46) $$

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

$$ f(x) = A [1 - \cos(\lambda x)] \quad (2D), \quad f(r) = A [1 - \cos(\lambda r)] \quad (\text{axisymmetric}). \qquad (3.1) $$

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
$$S_x = 1/\lambda, \quad S_h = A, \quad S_t = (\gamma\lambda^2 A)^{-1}, \qquad (3.2)$$

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 dimensionless surface height, defined as

$$f^* = f(x)/L_c = f\lambda = \lambda A [1 - \cos(\lambda x)], \qquad (3.3)$$

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,

$$\epsilon = |df/dx| = |\lambda A \sin(\lambda x)| \ll 1. \qquad (3.4)$$

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

$$h^* = h_n(x)/L_c \approx h_z \lambda \ll 1, \qquad (3.5)$$

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

$$F(X) = 1 - \cos(X) (2D), \quad F(R) = 1 - \cos(R) \text{ (axisymmetric).} \qquad (3.6)$$

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

$$\frac{\partial H}{\partial T} = \frac{\partial}{\partial X} \left[ H \left( \frac{\partial H}{\partial X} + \frac{dF}{dX} \right) \right], \qquad (3.7)$$

subject to the volume constraint

$$\int_{0}^{L} H \, dX = N, \qquad (3.8)$$

and boundary conditions

$$\left(\frac{\partial H}{\partial X} + \frac{dF}{dX}\right)\bigg|_{X=0} = 0 \text{ (no flux condition),} \quad H(X_F, T) = H(L, T) = 0, \qquad (3.9)$$

and the initial condition

$$H(X, 0) = \begin{cases} H_0, & 0 \le X \le N/H_0, \\ 0, & X > N/H_0. \end{cases} \qquad (3.10)$$

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

$$\frac{\partial H}{\partial T} - \frac{1}{R} \frac{\partial}{\partial R} \left[ R H \left( \frac{\partial H}{\partial R} + \frac{dF}{dR} \right) \right] = 0, \qquad (3.11)$$

subject to the volume constraint

$$2\pi \int_{0}^{L} RH \, dR = N, \qquad (3.12)$$

the boundary conditions

$$\left(\frac{\partial H}{\partial R} + \frac{dF}{dR}\right)\bigg|_{R=0} = 0, \quad H(R_F, T) = H(L, T) = 0, \qquad (3.13)$$
and the initial condition

$$H(R, 0) = \begin{cases} H_0, & 0 \le R \le \sqrt{N/(\pi H_0)}, \\ 0, & R > \sqrt{N/(\pi H_0)}. \end{cases} \qquad (3.14)$$

## 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

$$N = q/\phi A_s, \quad \text{with} \quad A_s = \int_{0}^{2\pi/\lambda} A [1 - \cos(\lambda x)] \, dx = 2\pi A/\lambda, \qquad (3.15)$$

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

For the axisymmetric case, we define

$$N = 4\pi^3 q/\phi V_s, \quad \text{with} \quad V_s = 2\pi \int_{0}^{2\pi/\lambda} r A [1 - \cos(\lambda r)] \, dr = 4\pi^3 A/\lambda^2, \qquad (3.16)$$

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 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 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 (1981), is

$$\nabla \cdot \mathbf{u} = 0, \qquad (4.1)$$

$$\frac{\partial (\rho \mathbf{u})}{\partial t} + \nabla \cdot (\rho \mathbf{u}\mathbf{u}) = -\nabla p + \mu \nabla^2 \mathbf{u} + \rho \mathbf{g} + \sigma \nabla \cdot \left( \frac{\nabla \alpha}{|\nabla \alpha|} \right) \nabla \alpha, \qquad (4.2)$$

$$\frac{\partial \alpha}{\partial t} + \nabla \cdot (\mathbf{u}\alpha) + \nabla \cdot \left[ \alpha (1-\alpha) \left( C |\mathbf{u}| \frac{\nabla \alpha}{|\nabla \alpha|} \right) \right] = 0. \qquad (4.3)$$
This model adopts a single-field formulation where both phases share common velocity (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 1988; Whitaker 1998) 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 2019). 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. 2015); (ii) a single-field formulation based on the total filtration velocity (Carrillo, Bourg & Soulaine 2020). The latter can be expressed as

$$ \nabla \cdot \mathbf{U} = 0, \qquad (4.4) $$

$$ \frac{\partial (\rho \mathbf{U})}{\partial t} + \nabla \cdot (\rho \mathbf{U} \mathbf{U}) = -\nabla P + \mu \nabla^2 \mathbf{U} + \rho \mathbf{g} - \phi k_0^{-1} \left( \frac{k_{r1}}{\mu_1} + \frac{k_{r2}}{\mu_2} \right)^{-1} \mathbf{U} \\ + \left[ k_0^{-1} \left( \frac{\mu_1}{k_{r1}} + \frac{\mu_2}{k_{r2}} \right)^{-1} (M_1\alpha_2 - M_2\alpha_1) \frac{\partial P_c}{\partial \alpha_1} - P_c \right] \nabla S_1, \qquad (4.5) $$

$$ \frac{\partial S_1}{\partial t} + \nabla \cdot (\mathbf{U}S_1) + \nabla \cdot (S_1 S_2 \mathbf{U}_r) = 0, \qquad (4.6) $$

$$ \mathbf{U}_r = -\phi^{-1} \left[ \left( \frac{M_1}{S_1} - \frac{M_2}{S_2} \right) \nabla P + \left( \frac{\rho_1 M_1}{S_1} - \frac{\rho_2 M_2}{S_2} \right) \mathbf{g} \right. \\ \left. + \left( \frac{M_1 S_2}{S_1} + \frac{M_2 S_1}{S_2} \right) \nabla P_c - \left( \frac{M_1}{S_1} - \frac{M_2}{S_2} \right) P_c \nabla S_1 \right]. \qquad (4.7) $$

This single-field formulation characterises macroscopic flow by the intrinsic averaged velocity $\mathbf{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:

$$ \nabla \cdot \mathbf{U} = 0, \qquad (4.8) $$

$$ \frac{\partial (\rho \mathbf{U})}{\partial t} + \nabla \cdot (\rho \mathbf{U}\mathbf{U}) = -\nabla P + \mu \nabla^2 \mathbf{U} + \rho \mathbf{g} - \phi \mu K^{-1} \mathbf{U} + \sigma_{eff} \nabla \cdot \left( \frac{\nabla S}{|\nabla S|} \right) \nabla S, \qquad (4.9) $$

$$ \frac{\partial S}{\partial t} + \nabla \cdot (\mathbf{U}S) + \nabla \cdot \left[ \phi S (1-S) \left( C |\mathbf{U}| \frac{\nabla S}{|\nabla S|} \right) \right] = 0. \qquad (4.10) $$

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® 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 2001; Lyle et al. 2005), 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} \text{ m}^2$. Glycerol is used as the working fluid, with a dynamic viscosity of $\mu = 1.5 \text{ Pa} \cdot \text{s}$ and density $\rho = 1.26 \times 10^3 \text{ 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 2024), with a wavelength $2\pi/\lambda = 0.05 \text{ m}$ and amplitude $A = 0.002 \text{ m}$. This configuration results in $A\lambda = 0.25 < 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
<figure><img src="image_4.png" /><figcaption>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.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>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.</figcaption></figure>

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® 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><img src="image_6.png" /><figcaption>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.</figcaption></figure>

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
<figure><img src="image_7.png" /><figcaption>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.</figcaption></figure>

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.

## 6. Discussion

The validity of low-dimensional models relies on three fundamental constraints. First, the thin-film approximation ensures that vertical pressure gradients dominate. If the fluid thickness becomes comparable to the characteristic horizontal length, the neglected dynamic terms would need to be incorporated. Second, the small-slope assumption underpins the simplifications made during the transformation between local and global coordinates. Third, the assumption of a small ratio between the surface amplitude and the characteristic length (i.e. for a sinusoidal surface with amplitude A and wavelength $L_c = 1/\lambda$, the condition $A\lambda \ll 1$ must hold) justifies neglecting curvature effects. Together, these constraints are key to the derivations presented earlier.

When the surface slopes or amplitudes become large, the underlying geometric approximations no longer hold. In such cases, the complete relationship between the vertical and normal interface heights, as well as the full curvature effects, requires the inclusion of higher-order terms. In fact, retaining the full curvature effects in the local equation introduces significant mathematical complexity when transforming the formulation to global coordinates. This complexity can hinder the consistent coupling of boundary conditions in the global framework, suggesting that a rigorous treatment of strong curvature effects is best confined to the local coordinate formulation. Nevertheless, in many practical situations, although local surface slopes or curvature may be significant, their effects decay rapidly with distance; thus, the violations of the underlying approximations remain localised, and the global model can still yield accurate predictions over the majority of the domain (e.g. the linear-exponential surfaces discussed in Appendix C).

In addition to the geometrical approximations, low-dimensional models are built upon several fundamental physical assumptions. First, a clear scale separation between the pore size and the curvilinear surface is required, which justifies treating the porous medium as a continuum and applying Darcy's law. Second, the use of a macroscopic sharp-interface approximation effectively reduces the complex multiphase system to a single-phase flow problem. This simplification is supported by the assumption that the released fluid is
significantly more viscous than the ambient fluid, thereby suppressing mixing, dissolution and diffusive processes at the interface. Although this means that the model cannot capture such interfacial phenomena, the resulting low-dimensional framework offers substantial computational efficiency compared with full multiphase CFD simulations. This efficiency is particularly advantageous for rapid assessments of large-scale scenarios where detailed numerical simulations would be prohibitively expensive.

We suggest that experimental validation of the low-dimensional models could be pursued using a Hele-Shaw cell, which has long served as an effective analogue for two-dimensional porous flows. In such an experiment, a precisely engineered curved bottom would be required, formed by two narrowly spaced parallel plates to create the flow cell. To maintain the creeping flow regime central to Darcy's law, the working fluid should be sufficiently viscous (e.g. glycerol or silicone oil) while still providing a substantial density contrast to drive gravity currents. The initial configuration must be carefully designed: the top gate should conform to the curved bottom profile, while the side walls remain vertical, to ensure the proper release of the initially rectangular fluid volume onto the curvilinear surface. The evolution of the interface can be monitored using imaging processing for quantitative comparison with the predictions of low-dimensional models. It should be noted, however, that experimental measurements may also capture effects not accounted for in the theoretical framework, such as surface tension, variations in the gap width and transient effects during rapid gate removal which could influence the early-stage released dynamics.

## 7. Relevance to CO₂ 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₂ geological sequestration. In practice, CO₂ 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₂ flow over wavy cap rocks in realistic geological formations. Such conditions are discussed in detail by Huppert & Neufeld (2014).

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

Figure 6 illustrates the propagation of supercritical CO₂ over both flat and wavy cap rocks under finite-volume release conditions. Unlike the continuous, monotonically
<figure><img src="image_8.png" /><figcaption>Figure 6. Comparison of supercritical CO₂ over horizontal flat and wavy cap rocks in both two-dimensional and axisymmetric propagations.</figcaption></figure>

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.

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₂ geological sequestration, where millions of cubic metres of supercritical CO₂ 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₂ within their troughs, but also extend the residence time of
supercritical CO₂, 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₂ sequestration.

## 8. Future outlook

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

One critical challenge is the migration of supercritical CO₂ through fractures or low-permeability zones in cap rocks. Although previous studies (Neufeld, Vella & Huppert 2009; Neufeld *et al.* 2011) 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 2014; Zheng *et al.* 2015). 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 2018; Benham, Bickle & Neufeld 2021), 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.* 2011; Golding, Huppert & Neufeld 2017). 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₂ 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₂ sequestration scenarios.

## 9. Conclusions

In this work, we have developed a comprehensive theoretical framework for porous media gravity currents propagating over curvilinear surfaces. By analysing the flow dynamics from local perspectives, we derived a local evolution equation that explicitly incorporates the effects of surface geometry through both slope and curvature terms. This local model reveals the interplay between the fluid and the underlying curvilinear topography.
Building on the local model, we performed a scaling analysis to transform into a global model. In doing so, we identified the key constraints, namely the thin-film approximation, small surface slopes and moderate surface amplitudes, that justify the neglect of higher-order geometric and curvature effects. As a consequence, the global model is predominantly governed by leading-order slope effects, thereby capturing the essential physics of the flow.

Our non-dimensionlisation further demonstrates that the evolution of finite-volume porous media gravity currents is primarily controlled by a key dimensionless parameter. This parameter, defined as the ratio of the released fluid volume to the characteristic volume of the curvilinear surface, delineates distinct flow regimes ranging from topography-dominated behaviour to volume-dominated dynamics.

The validity of low-dimensional models has been confirmed by numerical comparisons with computational fluid dynamics simulations based on a more sophisticated macroscopic sharp-interface flow model. The practical utility of low-dimensional models is applied to CO₂ geological sequestration scenarios. The results reveal that for finite-volume releases, wavy cap rock geometries can enhance trapping capacity compared with conventional flat-surface approximations, highlighting the critical role of topographic features in subsurface flow modelling. In the future, we will investigate the constant-flux condition, because this pattern more closely resembles those encountered in CO₂ sequestration to reflect the role of wavy cap rocks, thereby enhancing understanding on long-term CO₂ sequestration performance.

**Acknowledgements.** We are grateful to Professor M. Ungarish for thoughtful suggestions. We also thank the referees, who supplied a series of helpful insights.

**Declaration of interests.** The authors report no conflict of interest.

## Appendix A. Numerical validation on axisymmetric porous media gravity currents over sinusoidal surfaces

To numerically validate the axisymmetric propagation of porous media gravity currents, we employ an efficient approach that models only a thin wedge-shaped sector of the full three-dimensional domain, thereby substantially reducing computational overhead. In our simulations, a mesh comprising 480 000 cells was used. As shown in figures 7, 8 and 9, low-dimensional models and CFD simulations exhibit excellent agreement in axisymmetric propagation. Notably, under identical release criteria, the axisymmetric current propagates more slowly than in the two-dimensional case. These CFD results corroborate our analysis in the context of CO₂ geological sequestration, attributing the slower propagation to the increased volumetric demand required for radially spreading flows, which highlights the pronounced geometric influence on axisymmetric flow dynamics.

## Appendix B. Porous media gravity currents over horizontal flat surfaces

Before investigating porous media gravity currents over curvilinear surfaces, we consider both propagation over a horizontal flat surface, corresponding to the classical low-dimensional models (Huppert & Woods 1995; Lyle et al. 2005),

$$ \frac{\partial h}{\partial t} - \gamma \frac{\partial}{\partial x} \left( h \frac{\partial h}{\partial x} \right) = 0, \qquad (B1) $$

$$ \frac{\partial h}{\partial t} - \gamma \frac{1}{r} \frac{\partial}{\partial r} \left( r h \frac{\partial h}{\partial r} \right) = 0. \qquad (B2) $$
<figure><img src="image_9.png" /><figcaption>Figure 7. Numerical comparison of two approaches to current profiles for an axisymmetric porous media gravity current over a sinusoidal surface (Case 1). Dotted yellow line, numerical solution.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Figure 8. Numerical comparison of two approaches to current profiles for an axisymmetric porous media gravity current over a sinusoidal surface (Case 2). Dotted yellow line, numerical solution.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Figure 9. Numerical comparison of two approaches to current profiles for an axisymmetric porous media gravity current over a sinusoidal surface (Case 3). Dotted yellow line, numerical solution.</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>Figure 10. Numerical comparison of two approaches to current profiles for a porous media gravity current over a horizontal flat surface in both two-dimensional and axisymmetric propagations (t = 600 s). Dotted yellow line, numerical solution.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Figure 11. Numerical comparison of two approaches to dimensionless current fronts for a porous media gravity current over a horizontal flat surface in both two-dimensional and axisymmetric propagations.</figcaption></figure>

We conduct numerical simulations of a finite-volume glycerol release (0.05 m × 0.02 m) using both low-dimensional and CFD models. Figure 10 presents the current profiles obtained from both approaches at t = 600 s, illustrating the evolution from an initial square configuration to a porous media gravity current and demonstrating strong agreement between the models. Moreover, the close correspondence in the time-dependent dimensionless front position, as shown in figure 11, further validates our numerical framework prior to its extension to curvilinear surfaces. Minor discrepancies observed across three mesh resolutions can be attributed to the thin-film approximation employed in the low-dimensional models, which does not fully capture the vertical viscous dissipation during the initial fluid collapse, which is effectively addressed by the fully volume-averaged Navier-Stokes formulation in the CFD simulations.

## Appendix C. Porous media gravity currents over linear-exponential surfaces

We also explore the propagation over linear-exponential surfaces, as in our previous study on viscous gravity currents (Di & Huppert 2024). These surfaces are defined as

$$f(x) = axe^{-bx} \quad (2D), \quad f(r) = re^{-br} \quad (axisymmetric). \qquad (C1)$$

We identify the following scalings that emerge from the surface geometry:

$$S_x = 1/b, \quad S_h = a/b, \quad S_t = (\gamma ab)^{-1}. \qquad (C2)$$

Here, $1/b$ represents the characteristic horizontal length scale where the exponential decay becomes significant, while $a/b$ captures the characteristic height of the surface.
<figure><img src="image_14.png" /><figcaption>Figure 12. Numerical comparison of two approaches to current profiles for a porous media gravity current over a linear-exponential surface (t = 600 s). Dotted yellow line, numerical solution. (a) Two-dimensional propagation; (b) axisymmetric propagation.</figcaption></figure>

Accordingly, the dimensionless forms for these surfaces are

$$F(X) = Xe^{-X} \quad (2D), \quad F(R) = Re^{-R} \quad (axisymmetric), \quad (C3)$$

and the dimensionless released volumes are defined as

$$N = q/\phi A_e, \quad \text{with} \quad A_e = \int_{0}^{\infty} a x e^{-bx} dx = a/b^2 \quad (2D), \quad (C4)$$

$$N = 4\pi q/\phi V_e, \quad \text{with} \quad V_e = 2\pi \int_{0}^{\infty} ar^2 e^{-br} dr = 4\pi a/b^3 \quad (axisymmetric). \quad (C5)$$

For instance, by assigning the surface parameters $a = 1$ and $b = 10 \text{ m}^{-1}$, we examine three cases of fluid evolution over linear-exponential surfaces. Figure 12 qualitatively illustrates the current profiles in both two-dimensional and axisymmetric configurations at $t = 600$ s. In Case 1, a glycerol volume of $0.03 \text{ m} \times 0.02 \text{ m}$ ($X_0(R_0) = 0.3, H_0 = 0.2$) is released on the left side of the peak, forming a horizontal macroscopic sharp interface that accumulates in the left vacancy of the linear-exponential surface. Without a pressure source from the origin, the current eventually stops evolving. In Case 2, the release volume increases to $0.08 \text{ m} \times 0.02 \text{ m}$ ($X_0(R_0) = 0.8, H_0 = 0.2$), enabling the fluid to bypass the peak and begin flowing downslope. In Case 3, with a larger release volume of $0.15 \text{ m} \times 0.02 \text{ m}$ ($X_0(R_0) = 1.5, H_0 = 0.2$), the initial fluid extent exceeds the peak, resulting in a current that collects near the origin, while its front advances downslope.

Figure 13 quantitatively compares the current front between the low-dimensional and CFD models. The strong agreement between the two approaches across all cases further confirms the validity of our framework for porous media gravity currents over curvilinear surfaces. In Case 1, the current front remains nearly stationary over time, reflecting the stagnation within the left vacancy; in Cases 2 and 3, the advancing front indicates continuous propagation driven by the interplay between current thickness and surface slope. In summary, propagation over the linear-exponential surface provides a local view of the relationship between the released volume and the peak, leading to various evolutions.

It is worth noting that although the linear-exponential surface exhibits steep slopes near the origin, i.e. potentially violating the small slope assumption of our theoretical framework, the numerical validation confirms that the low-dimensional models remain
<figure><img src="image_15.png" /><figcaption>Figure 13. Numerical comparison of two approaches to dimensionless current fronts for a porous media gravity current over a linear-exponential surface in both two-dimensional and axisymmetric propagations.</figcaption></figure>

robust. This robustness arises because the region of steep slopes is limited to the origin, which is smaller than the overall horizontal extent of the current. Consequently, the primary assumption remains valid at the global scale, ensuring that the main flow dynamics is captured and that the low-dimensional models yield reliable predictions.

## Appendix D. Extension to generalised surfaces

### D.1. Axisymmetric power-exponential surfaces

We define a generalised axisymmetric power-exponential surface as

$$f(r) = ar^m e^{-br^n}, \qquad (D1)$$

where the parameters $m$ ($m > 1$) and $n$ ($n > 1$) are constants that determine the specific shape of the surface. Using the scalings

$$S_x = (1/b)^{1/n}, \quad S_h = a(S_r)^m, \quad S_t = (S_r)^{2-m}/(\beta a), \qquad (D2)$$

the dimensionless model is written as

$$\frac{\partial H}{\partial T} - \frac{1}{R} \frac{\partial}{\partial R} \left[ R H \left( \frac{\partial H}{\partial R} + \frac{dF}{dR} \right) \right] = 0, \qquad (D3)$$

$$2\pi \int_{0}^{L} R H \, dR = N, \qquad (D4)$$

$$\left( \frac{\partial H}{\partial R} + \frac{dF}{dR} \right) |_{R=0} = 0, \qquad (D5)$$

$$H(R_F, T) = H(L, T) = 0, \qquad (D6)$$

$$H(R, 0) = \begin{cases} H_0 & 0 \le R \le \sqrt{N/(\pi H_0)}, \\ 0 & R > \sqrt{N/(\pi H_0)}, \end{cases} \qquad (D7)$$

where the dimensionless power-exponential surface is

$$F(R) = R^m e^{-R^n}. \qquad (D8)$$
The definition of N is

$$N = \frac{2\pi \Gamma \left[ \frac{(m+2)}{n} \right] q}{\phi n V_{ge}}, \qquad (D9)$$

with the volume of the power-exponential surface

$$V_{ge} = 2\pi \int_{0}^{\infty} r a^{m} e^{-br^{n}} dr = \frac{2\pi a}{nb^{\left(m+2\right)/n}} \Gamma \left[ \frac{(m+2)}{n} \right]. \qquad (D10)$$

## D.2. Arbitrary curvilinear surfaces

Extending sinusoidal surfaces to arbitrary axisymmetric continuous surfaces

$$f(r) = \sum_{i=1}^{j} A_i [1 - \cos(\lambda_i r)]. \qquad (D11)$$

Letting $j = 3$, the curvilinear surface is written as

$$f(r) = A_1 [1 - \cos(\lambda_1 r)] + A_2 [1 - \cos(\lambda_2 r)] + A_3 [1 - \cos(\lambda_3 r)], \qquad (D12)$$

where $\lambda_3 > \lambda_2 > \lambda_1$. Using the scalings

$$S_x = 1/\lambda_1, \quad S_h = A_1, \quad S_t = 1/(\beta\lambda_1^2 A_1), \qquad (D13)$$

we write the dimensionless model as

$$\frac{\partial H}{\partial T} - \frac{1}{R} \frac{\partial}{\partial R} \left[ R H \left( \frac{\partial H}{\partial R} + \frac{dF}{dR} \right) \right] = 0, \qquad (D14)$$

$$2\pi \int_{0}^{L} R H \, dR = N, \qquad (D15)$$

$$\left(\frac{\partial H}{\partial R} + \frac{dF}{dR}\right) |_{R=0} = 0, \qquad (D16)$$

$$H(R_F, T) = H(L, T) = 0, \qquad (D17)$$

$$H(R, 0) = \begin{cases} H_0 & 0 \le R \le \sqrt{N/(\pi H_0)}, \\ 0 & R > \sqrt{N/(\pi H_0)}, \end{cases} \qquad (D18)$$

where the dimensionless curvilinear surface is

$$F(R) = 1 - \cos(R) + P_2 [1 - \cos(Q_2 R)] + P_3 [1 - \cos(Q_3 R)], \qquad (D19)$$

$$P_2 = A_2/A_1, \quad Q_2 = \lambda_2/\lambda_1, \qquad (D20)$$

$$P_3 = A_3/A_1, \quad Q_3 = \lambda_3/\lambda_1. \qquad (D21)$$

The definition of N is

$$N = 4\pi^3 aq/\phi V_{s1}, \qquad (D22)$$
where $V_{s1}$ is the volume given by the first wavelength of the primary term, $A_1[1 - \cos(\lambda_1 r)]$, calculated as

$$V_{s1} = 2\pi A_1 \int_{0}^{2\pi/\lambda_1} r[1 - \cos(\lambda_1 r)] dr = 4\pi^3 A_1 / \lambda_1^2. \quad (D23)$$

Similarly to our previous study on viscous gravity currents over curvilinear surfaces (Di & Huppert 2024), we reveal that the flow system for a finite-volume viscous fluid in porous media over a curvilinear surface is also jointly governed by the dimensionless surface shape $F(R)$ and the crucial dimensionless number $N$.

REFERENCES

ACTON, J.M., HUPPERT, H.E. & WORSTER, M.G. 2001 Two-dimensional viscous gravity currents flowing over a deep porous medium. *J. Fluid Mech.* **440**, 359–380.

BENHAM, G.P., BICKLE, M.J. & NEUFELD, J.A. 2021 Two-phase gravity currents in layered porous media. *J. Fluid Mech.* **922**, A7.

CARRILLO, F.J., BOURG, I.C. & SOULAINE, C. 2020 Multiphase flow modeling in multiscale porous media: an open-source micro-continuum approach. *J. Comput. Phys.* **8**, 100073.

DAVIT, Y. & QUINTARD, M. 2019 One-phase and two-phase flow in highly permeable porous media. *Heat. Transfer Engng* **40** (5-6), 391–409.

DI, X. & HUPPERT, H.E. 2024 Release of a finite volume of viscous fluid over rigid curvilinear surfaces. *J. Fluid Mech.* **1001**, A19.

DI FEDERICO, V., ARCHETTI, R. & LONGO, S. 2012 Similarity solutions for spreading of a two-dimensional non-Newtonian gravity current in a porous layer. *J. Non-Newtonian Fluid Mech.* **177**, 46–53.

GOLDING, M.J., HUPPERT, H.E. & NEUFELD, J.A. 2017 Two-phase gravity currents resulting from the release of a fixed volume of fluid in a porous medium. *J. Fluid Mech.* **832**, 550–577.

GOLDING, M.J., NEUFELD, J.A., HESSE, M.A. & HUPPERT, H.E. 2011 Two-phase gravity currents in porous media. *J. Fluid Mech.* **678**, 248–270.

HINTON, E.M. & WOODS, A.W. 2018 Buoyancy-driven flow in a confined aquifer with a vertical gradient of permeability. *J. Fluid Mech.* **848**, 411–429.

HIRT, C.W. & NICHOLS, B.D. 1981 Volume of fluid (VOF) method for the dynamics of free boundaries. *J. Comput. Phys.* **39** (1), 201–225.

HORGUE, P., SOULAINE, C., FRANC, H., GUIBERT, R. & DEBENEST, G. 2015 An open-source toolbox for multiphase flow in porous media. *Comput. Phys. Commun.* **187**, 217–226.

HUPPERT, H.E. & NEUFELD, J.A. 2014 The fluid mechanics of carbon dioxide sequestration. *Annu. Rev. Fluid Mech.* **46** (1), 255–272.

HUPPERT, H.E. & WOODS, A.W. 1995 Gravity-driven flows in porous layers. *J. Fluid Mech.* **292**, 55–69.

LONGO, S., DI FEDERICO, V., CHIAPPONI, L. & ARCHETTI, R. 2013 Experimental verification of power-law non-newtonian axisymmetric porous gravity currents. *J. Fluid Mech.* **731**, R2.

LYLE, S., HUPPERT, H.E., HALLWORTH, M., BICKLE, M. & CHADWICK, A. 2005 Axisymmetric gravity currents in a porous medium. *J. Fluid Mech.* **543**, 293–302.

NEUFELD, J.A., VELLA, D. & HUPPERT, H.E. 2009 The effect of a fissure on storage in a porous medium. *J. Fluid Mech.* **639**, 239–259.

NEUFELD, J.A., VELLA, D., HUPPERT, H.E. & LISTER, J.R. 2011 Leakage from gravity currents in a porous medium. Part 1. A localized sink. *J. Fluid Mech.* **666**, 391–413.

PEGLER, S.S., HUPPERT, H.E. & NEUFELD, J.A. 2013 Topographic controls on gravity currents in porous media. *J. Fluid Mech.* **734**, 317–337.

PEGLER, S.S., HUPPERT, H.E. & NEUFELD, J.A. 2014 Fluid migration between confined aquifers. *J. Fluid Mech.* **757**, 330–353.

PHILIP, J.R. 1970 Flow in porous media. *Annu. Rev. Fluid Mech.* **2** (1), 177–204.

QUINTARD, M. & WHITAKER, S. 1988 Two-phase flow in heterogeneous porous media: the method of large-scale averaging. *Transport Porous Media* **3** (4), 357–413.

VELLA, D. & HUPPERT, H.E. 2006 Gravity currents in a porous medium at an inclined plane. *J. Fluid Mech.* **555**, 353–362.

WHITAKER, S. 1998 *the Method of Volume Averaging*. Springer Science & Business Media.

WOODS, A.W. 2015 *Flow in Porous Rocks*. Cambridge University Press.

ZHENG, Z., GUO, B., CHRISTOV, I.C., CELIA, M.A. & STONE, H.A. 2015 Flow regimes for fluid injection into a confined porous medium. *J. Fluid Mech.* **767**, 881–909.