Reservoir Geometry: Riemannian Manifolds in Oil and Gas

Beneath the surface of every oil and gas field lies a porous rock formation, and inside that rock, fluid moves. The simplest model for that motion, Darcy's law, was measured empirically by Henry Darcy in 1856 from sand filtration experiments on the Dijon fountain supply. He found that volumetric flow rate through a porous column is proportional to the pressure drop and inversely proportional to the column length, a linear constitutive relation whose proportionality constant, the permeability, characterises the rock. For more than a century, reservoir engineers treated permeability as a scalar, a single number per rock sample, and built analytic solutions for pressure propagation on that assumption.

Real reservoirs are three-dimensional and anisotropic. Sedimentary layering, fracture networks, and depositional fabrics all impose directional preferences on fluid flow: the rock is more permeable along bedding planes than across them, more permeable along open fractures than perpendicular to them, more permeable horizontally than vertically by factors routinely exceeding one hundred. A scalar permeability misrepresents all of this. The correct object is a permeability tensor: a symmetric positive definite contravariant tensor of type $(2,0)$ whose coordinate matrix in a given chart is a $3 \times 3$ symmetric positive definite array assigning a different resistance to flow in each direction.

That tensor is the entry point for Riemannian geometry. A symmetric positive definite matrix is, simultaneously, a linear map and a Riemannian metric on the ambient space. Once the tensor field is read as a metric, the entire machinery of differential geometry becomes available: geodesics, curvature, the Laplace–Beltrami operator, the heat kernel. The pressure equation is a heat equation on a Riemannian manifold. Streamlines are geodesics. Permeability contrasts are Riemannian refractive index discontinuities, bending streamlines by an exact analogue of Snell's law. The geometry is the exact mathematical content of the equations.

This post develops that content in three dimensions from first principles. We begin by separating the permeability tensor from its coordinate matrix representation, introducing index notation, the Levi-Civita connection, and the fibre bundle structure of material property fields (§2). We construct the 3D permeability metric and its diffusion ellipsoid, and show that Darcy flux is a Riemannian gradient (§3). We derive the 3D Laplace-Beltrami heat equation, the 3D spherical pressure kernel, and the spherical-to-radial flow transition that governs early-time well test signatures (§4). We compute the 3D Christoffel symbols for a diagonal tensor, derive the Darcy-Snell refraction law at permeability interfaces, and develop the Jacobi field equation for streamline spreading (§5). We prove the 3D scalar curvature formula $R = 2\Delta k - \tfrac{5}{2}|\nabla k|^2/k$, show that the Weyl curvature vanishes identically in 3D, and connect curvature to pressure transient anomalies (§6). We solve the 3D spectral problem on a rectangular parallelepiped and extract the pseudosteady-state time constant for transversely isotropic formations (§7). We treat the space of $3 \times 3$ permeability tensor coordinate matrices as the manifold $\mathrm{SPD}(3)$ and derive the geodesic interpolation between measurements, including the fibre bundle and frame-bundle structure (§8). We extend the framework to a deforming reservoir, deriving the Biot poroelastic equations in abstract index form, establishing the Kozeny-Carman permeability update as a conformal metric flow, and proving that the full coupled Biot-Darcy system is a PDE on a time-parametrised Riemannian manifold (§9). We close by showing how all of this executes in no_std Rust on a bare-metal microcontroller, using the cartan geometry library, for real-time well-site monitoring (§10).

1. Darcy's Law and the 3D Permeability Tensor

1.1. The Empirical Law

Darcy's original experiment measured the steady-state volumetric flow rate $Q$ through a vertical sand column of cross-sectional area $A$ and length $L$, subject to a pressure head difference $\Delta h$:

where $\kappa$ is the hydraulic conductivity of the sand. Dividing both sides by $A$ and taking the continuum limit gives the volumetric flux (velocity per unit area), the Darcy velocity or specific discharge:

where $p$ is the fluid pressure, $\mu$ is the dynamic viscosity, and $k$ is the intrinsic permeability. This isotropic form holds only when the rock has no preferred flow direction, a condition satisfied in clean, well-sorted sandstone but almost nowhere else.

1.2. The Full 3D Permeability Tensor

In an anisotropic three-dimensional formation, the Darcy velocity in each direction depends on the pressure gradient in all three directions. The correct constitutive law is:

where $\mathbf{K}$ is the permeability tensor. Symmetry ($k_{ij} = k_{ji}$) follows from the Onsager reciprocal relations: the cross-coupling between a pressure gradient in direction $i$ and the resulting flux in direction $j$ must equal the coupling in the reverse direction, a consequence of microscopic reversibility. Positive definiteness is a physical requirement: the viscous dissipation rate per unit volume, $\mathbf{u}\cdot\nabla p = -(1/\mu)(\nabla p)^\top\mathbf{K}\nabla p$, must be non-negative for any pressure gradient.

The positive definiteness conditions for the $3 \times 3$ tensor are:

These are the three Sylvester conditions (leading principal minors all positive), which characterise symmetric positive definite matrices.

1.3. Transverse Isotropy: The Standard Reservoir Model

The dominant anisotropy type in clastic and carbonate reservoirs is transverse isotropy (TI), in which the rock has equal permeability in all horizontal directions but a different permeability vertically. With $z$ as the vertical axis:

where $k_h$ is the horizontal permeability and $k_v$ is the vertical permeability. All off-diagonal entries vanish when the coordinate axes are aligned with the bedding plane. The TI model reduces the full six-parameter tensor to two independent parameters, which is sufficient to capture the dominant flow behaviour in most conventional reservoirs.

1.4. Continuity and the 3D Pressure Equation

Combining Darcy's law with conservation of mass for a slightly compressible fluid (constant total compressibility $c_t$) in a porous medium with porosity $\phi$:

Substituting [eq:darcy-tensor] into [eq:continuity] gives the 3D pressure diffusion equation:

For constant $\mathbf{K}$, this is a constant-coefficient parabolic PDE. For heterogeneous $\mathbf{K}(\mathbf{x})$, it is a variable-coefficient elliptic-parabolic operator whose structure is governed by the Riemannian geometry of the tensor field.

2. Tensors, Coordinates, and the Permeability Object

2.1. The Permeability Tensor as a Coordinate-Free Object

The word "tensor" is often used interchangeably with "matrix" in engineering literature. In what follows that identification is not merely imprecise: it is wrong in a way that produces incorrect equations on a deforming manifold. The distinction is therefore introduced carefully here, before any geometry is built on top of it.

A tensor field of type $(r, s)$ on a smooth manifold $\mathcal{R}$ assigns to each point $x \in \mathcal{R}$ a multilinear map

where $T_x\mathcal{R}$ is the tangent space at $x$ and $T_x^*\mathcal{R}$ is the cotangent space. The tensor itself is this map. A matrix is the array of numbers obtained by evaluating the map on a chosen basis. Choosing a different basis changes the matrix; the tensor does not change.

The permeability tensor $K$ is a symmetric contravariant tensor of type $(2, 0)$: at each $x$ it is a symmetric bilinear map $K_x \colon T_x^*\mathcal{R} \times T_x^*\mathcal{R} \to \mathbb{R}$. Given a coordinate chart $\{x^i\}_{i=1}^{3}$ on $\mathcal{R}$ with coordinate basis $\{\partial_i\}$ for $T_x\mathcal{R}$ and dual basis $\{dx^i\}$ for $T_x^*\mathcal{R}$, the components of $K$ in that chart are

The $3 \times 3$ array $[K^{ij}(x)]$ is the permeability matrix in the chart $\{x^i\}$. Darcy's law in this chart is

where $u^i$ are the contravariant components of the Darcy velocity vector and $\partial_j p$ are the covariant components of the pressure gradient covector. The repeated index $j$ is summed from 1 to 3 (Einstein convention). Throughout this post, upper indices denote contravariant components and lower indices denote covariant components.

2.2. Transformation Law and the Primacy of the Tensor

Under a coordinate change $x^i \to \tilde{x}^k$ with Jacobian $J^k{}_i = \partial\tilde{x}^k/\partial x^i$, the components of a $(2,0)$ tensor transform as

Both indices transform with one factor of $J$ each, reflecting the contravariant character. The tensor $K$ is unchanged: it is the same physical object, described by different numbers in the new chart. A rotation from geographic coordinates (north, east, down) to borehole coordinates (along-hole, cross-hole, vertical) multiplies the matrix $[K^{ij}]$ by two factors of the rotation matrix, but the permeability of the rock at that point is unaltered.

2.3. The Metric as a Covariant Tensor

The Riemannian metric $g$ derived from $K$ is a $(0, 2)$ covariant symmetric tensor: at each $x$, it is a symmetric positive definite bilinear map $g_x \colon T_x\mathcal{R} \times T_x\mathcal{R} \to \mathbb{R}$. Its components in a chart are

The relationship $g = K^{-1}$ means: the matrix $[g_{ij}]$ is the matrix inverse of $[K^{ij}]$. More invariantly, the contraction $g_{ik}\,K^{kj} = \delta^j_i$ holds in every chart, where $\delta^j_i$ is the Kronecker delta. Indices are raised and lowered by $K^{ij}$ and $g_{ij}$:

The Darcy velocity has contravariant components $u^i$ (a vector, pointing in the flow direction); the pressure gradient has covariant components $\partial_i p$ (a covector, pointing in the direction of steepest pressure increase). The permeability tensor $K^{ij}$ converts covectors to vectors: it maps the pressure gradient covector to the flux vector.

2.4. The Levi-Civita Connection

A connection on $(\mathcal{R}, g)$ is a rule for differentiating tensor fields that is consistent with the manifold structure. The Levi-Civita connection $\nabla$ is the unique connection that is (i) metric-compatible, $\nabla_k g_{ij} = 0$, and (ii) torsion-free, $\Gamma^k{}_{ij} = \Gamma^k{}_{ji}$. Its Christoffel symbols in a coordinate chart are

Note: $\Gamma^k{}_{ij}$ are not the components of a tensor; they transform inhomogeneously under coordinate changes. They are the correction terms that account for the curvature of the coordinate grid, and they vanish in a geodesic normal chart centred at any given point. The covariant derivative of a vector field $v^i$ is

and its covariant divergence is

The last equality follows from the identity $\Gamma^i{}_{ij} = \partial_j\ln\sqrt{\det g}$ and is the form used in the Laplace-Beltrami operator of §4.

2.5. Why the Distinction Is Critical on a Deforming Manifold

On the static reservoir of §3–§8, the tensor/matrix distinction is important for correctness but does not change the form of the equations: the domain $\mathcal{R}$ is fixed, and a single coordinate chart suffices. On a deforming reservoir (§9), the solid skeleton displaces with a smooth map $\varphi_t \colon \mathcal{R}_0 \to \mathcal{R}(t)$ from the reference configuration $\mathcal{R}_0 = \mathcal{R}(0)$ to the current configuration $\mathcal{R}(t)$. The deformation gradient is

where $X^J$ are material (Lagrangian) coordinates and $x^i$ are spatial (Eulerian) coordinates. Under this map, the coordinate chart on $\mathcal{R}(t)$ in which the permeability matrix is measured changes at every time step.

The permeability tensor $K$ is a material property of the rock: it does not change when coordinates are relabelled. Its Eulerian matrix representation $[K^{ij}(x, t)]$ does change, because the coordinate axes $\{\partial_i\}$ at a point $x \in \mathcal{R}(t)$ are not the same objects as the axes $\{\partial_J\}$ at the corresponding material point $X \in \mathcal{R}_0$. The correct object to transport forward in time is not the matrix $[K^{ij}]$ but the push-forward of the tensor:

where $K_0^{JL}$ are the components of the permeability tensor in the reference configuration. This is the contravariant transformation law applied to the deformation map. Confusing the matrix $[K^{ij}]$ with the tensor $K$ on a deforming domain is equivalent to forgetting to account for the Jacobian $F$, leading to an error of order $|F - I|$ in the permeability (which equals the volumetric strain, a measurable and sometimes large quantity in compaction reservoirs). Section §9 develops this in full.

3. The 3D Reservoir as a Riemannian Manifold

3.1. The Permeability Metric

A symmetric positive definite $(2,0)$ permeability tensor $K$ defines a Riemannian metric $g$ on $\mathcal{R}$ by declaring $g = K^{-1}$: the components of the covariant metric tensor in any chart are $g_{ij} = (K^{-1})_{ij}$, the entries of the matrix inverse of $[K^{ij}]$. The following theorem shows why this choice is canonical.

Proof. The Riemannian gradient is defined as the unique vector field satisfying

for all tangent vectors $\mathbf{v}$. In coordinates, $g(\mathbf{w}, \mathbf{v}) = g_{ij} w^i v^j = (K^{-1})_{ij} w^i v^j$, so the condition reads $g_{ij}\,(\mathrm{grad}_g\,p)^i = \partial_j p$, i.e., $(K^{-1})_{ij}\,(\mathrm{grad}_g\,p)^i = \partial_j p$ for all $j$. Multiplying both sides by $K^{jk}$ and summing over $j$ gives $(\mathrm{grad}_g\, p)^k = K^{kj}\partial_j p$, i.e., $\mathrm{grad}_g\, p = \mathbf{K}\nabla p$. Therefore

The proof uses no properties specific to $n = 3$; it holds in any dimension. $\square$

3.2. Transverse Isotropy and the Diffusion Ellipsoid

For the transversely isotropic tensor $\mathbf{K}_\mathrm{TI} = \mathrm{diag}(k_h, k_h, k_v)$, the metric is

The Riemannian unit ball $\{(x,y,z) : g_{ij}x^i x^j \leq 1\}$ is the ellipsoid

with horizontal semi-axis $\sqrt{k_h}$ and vertical semi-axis $\sqrt{k_v}$. This is the diffusion ellipsoid: it describes the shape of the region reached by a unit-energy pressure pulse in unit Riemannian time. Since $k_h > k_v$ in most reservoirs, the ellipsoid is oblate (flattened vertically), reflecting the greater ease of horizontal flow.

3.3. Volume Form and Riemannian Divergence

The Riemannian volume form on $(\mathcal{R}, g)$ is

Here $\det g = \det[g_{ij}]$ is the determinant of the matrix of covariant metric components, which equals $1/\det[K^{ij}]$ since $[g_{ij}]$ is the matrix inverse of $[K^{ij}]$. For the TI tensor, $\det\mathbf{K}_\mathrm{TI} = k_h^2 k_v$, so $\mathrm{d}V_g = (k_h^2 k_v)^{-1/2}\mathrm{d}^3 x$. High-permeability regions have a smaller Riemannian volume element: they are geometrically compressed. The Riemannian divergence of a vector field $\mathbf{F}$ is

4. Pressure Diffusion as the 3D Laplace–Beltrami Equation

4.1. The 3D Laplace–Beltrami Operator

4.2. Constant vs. Heterogeneous Permeability

Proof. When $\mathbf{K}$ is constant, all derivatives of $K^{ij}$ and $\det\mathbf{K}$ vanish, so [eq:lb-explicit] reduces to $\Delta_g p = K^{ij}\,\partial_i\partial_j p = \nabla\cdot(\mathbf{K}\nabla p)$. The pressure equation then takes the form [eq:heat-beltrami]. $\square$

4.3. Heat Kernels in 2D and 3D

The fundamental solution to the heat equation $\partial_t p = D\Delta p$ depends critically on the spatial dimension. For isotropic permeability $k$ and diffusivity $D = k/(\phi\mu c_t)$, the heat kernels in two and three dimensions are:

where $r = |\mathbf{x}|$. The 2D kernel applies to an infinite vertical well producing from a formation of uniform thickness (the pressure diffusion is planar, independent of the $z$-coordinate). The 3D kernel applies to a point source in a three-dimensional infinite medium.

Proof. The result follows by integrating the 3D heat kernel against the source. The key identities are: the Laplacian of $1/r$ in 3D is $\Delta(1/r) = -4\pi\delta(\mathbf{x})$, and the convolution of the heat kernel with $1/(4\pi r)$ gives [eq:3d-drawdown]. Explicitly, verify the two limits: as $t \to 0$, $\mathrm{erfc}(r/\sqrt{4Dt}) \to 0$ since the argument goes to infinity, so $\Delta p \to 0$ (no drawdown before production). As $t \to \infty$, $\mathrm{erfc}(0) = 1$, giving the steady-state Newtonian potential $\Delta p \to q\mu/(4\pi k r)$, which satisfies the time-independent Laplace equation $\Delta(\Delta p) = 0$ at $r > 0$, confirming convergence to steady state. $\square$

4.4. The Spherical-to-Radial Flow Transition

The transition between 3D spherical and 2D radial flow is a concrete Riemannian geometry result: the dimension of the heat kernel changes as the pressure diffusion ellipsoid transitions from three-dimensional to quasi-two-dimensional.

Proof sketch. The TI metric maps the reservoir to a rescaled space with coordinates $X = x/\sqrt{k_h}$, $Y = y/\sqrt{k_h}$, $Z = z/\sqrt{k_v}$, in which the metric is Euclidean (shown in §5.2). In this space, the Riemannian distance from the well axis is $(X^2+Y^2+Z^2)^{1/2} = r_\mathrm{eff}/\sqrt{k_h}$. At early time, before the diffusion front is constrained by the finite well length, the 3D heat kernel applies in the rescaled space, giving [eq:spherical-regime]. At late time, the $Z$-direction has equilibrated over the well length and the problem reduces to a 2D Laplacian in the $(X, Y)$-plane, recovering the Theis solution [eq:radial-regime]. $\square$

5. 3D Streamlines as Geodesics

5.1. Christoffel Symbols for Diagonal 3D K

For a diagonal permeability tensor $\mathbf{K} = \mathrm{diag}(k_1, k_2, k_3)$ in three dimensions (with coordinates $x^1, x^2, x^3$), the metric is $g_{ij} = \delta_{ij}/k_i$ (no sum). The Levi-Civita Christoffel symbols are

where the outer index $\ell$ does not participate in a sum (since $g$ is diagonal, $g^{\ell m} = k_\ell \delta^{\ell m}$).

Proof. Compute each case from [eq:christoffel] using $g_{\alpha\beta} = \delta_{\alpha\beta}/k_\alpha$ and $g^{\alpha\beta} = k_\alpha\delta^{\alpha\beta}$.

For $\ell = i = j = \alpha$: $\Gamma^\alpha_{\alpha\alpha} = (k_\alpha/2)(2\partial_\alpha(1/k_\alpha) - \partial_\alpha(1/k_\alpha)) = (k_\alpha/2)\partial_\alpha(1/k_\alpha) = -\partial_\alpha k_\alpha/(2k_\alpha)$.

For $\ell = j = \beta$, $i = \alpha$, $\alpha \neq \beta$: $\Gamma^\beta_{\alpha\beta} = (k_\beta/2)(\partial_\alpha(1/k_\beta) + \partial_\beta(1/k_\alpha) - \partial_\beta(1/k_\alpha)) = (k_\beta/2)\partial_\alpha(1/k_\beta) = -\partial_\alpha k_\beta/(2k_\beta)$.

For $\ell = \beta$, $i = j = \alpha$, $\alpha \neq \beta$: $\Gamma^\beta_{\alpha\alpha} = (k_\beta/2)(0 + 0 - \partial_\beta(1/k_\alpha)) = (k_\beta/2)(\partial_\beta k_\alpha/k_\alpha^2) = k_\beta\partial_\beta k_\alpha/(2k_\alpha^2)$.

For constant $\mathbf{K}$, all partial derivatives of $k_\alpha$ vanish and so all Christoffel symbols are zero. The geodesic equation $\ddot\gamma^\ell + \Gamma^\ell_{ij}\dot\gamma^i\dot\gamma^j = 0$ reduces to $\ddot\gamma^\ell = 0$, giving straight-line solutions. $\square$

5.2. The 3D Rescaled Coordinate Transformation

For a constant diagonal tensor $\mathbf{K} = \mathrm{diag}(k_1, k_2, k_3)$, the coordinate change

(no summation) transforms the metric to the standard Euclidean form in three dimensions:

For the TI reservoir with $k_1 = k_2 = k_h$, $k_3 = k_v$, the rescaled domain is $[0, L_x/\sqrt{k_h}] \times [0, L_y/\sqrt{k_h}] \times [0, L_z/\sqrt{k_v}]$. The vertical dimension is stretched by the factor $\sqrt{\kappa}$, so a reservoir that is geometrically thin vertically may be effectively isotropic in the rescaled coordinates when $\kappa = k_h/k_v$ is large enough.

5.3. The Darcy–Snell Law at Permeability Interfaces

When a streamline crosses a planar boundary between two isotropic regions with permeabilities $k_1$ and $k_2$, it does not continue in the same direction. The angle of incidence and the angle of refraction are related by an exact analogue of Snell's law from optics.

Proof. By Fermat's principle applied to the Riemannian metric, geodesics minimise the arc-length functional $\int\mathrm{d}s_g = \int |\mathrm{d}\mathbf{x}|/\sqrt{k}$. This is formally identical to optical path minimisation with refractive index $n = 1/\sqrt{k}$: Snell's law states $n_1\sin\theta_1 = n_2\sin\theta_2$, which translates to [eq:darcy-snell]. More directly: since $k = k(z)$ depends only on $z$, the Lagrangian $L = (1/2k)(|\dot{\mathbf{x}}|^2)$ is independent of the tangential coordinates $(x,y)$. The Euler–Lagrange equations give the conserved quantities $\dot{x}/k = p_x$ and $\dot{y}/k = p_y$. Writing the tangential speed $|v_T|^2 = \dot{x}^2 + \dot{y}^2 = |\dot{\mathbf{x}}|^2\sin^2\theta$ and the Riemannian arc-length speed $|\dot{\mathbf{x}}|/\sqrt{k} = C$ (a constant for arc-length parametrisation), the conserved tangential co-momentum is $p_T = |\dot{\mathbf{x}}|\sin\theta/k = C\sin\theta/\sqrt{k}$. Setting this equal across the interface gives [eq:darcy-snell]. Total reflection occurs when no real $\theta_2$ satisfies [eq:darcy-snell], i.e., when $\sin\theta_1 > \sqrt{k_2/k_1}$. $\square$

5.4. Jacobi Fields and Streamline Spreading

In three dimensions, the transverse neighbourhood of a geodesic is two-dimensional. The Jacobi equation governs how nearby geodesics (streamlines) spread apart:

where $D/ds = \nabla_{\dot\gamma}$ is the covariant derivative along the geodesic $\gamma$, $J$ is a Jacobi field (transverse deviation vector, $g(J, \dot\gamma) = 0$), and $R$ is the Riemann curvature tensor. In 3D, there are two independent transverse directions, and two independent Jacobi fields $J_1, J_2$ govern the spreading of the flow tube cross-section.

The Jacobi field equation is the 3D geometric description of the phenomenon that practitioners call flow channelling: preferential flow through high-permeability streaks or fractures, resulting in early breakthrough in secondary recovery operations. The curvature of the permeability field determines whether and how rapidly channelling develops along a given flow path.

6. Curvature of the 3D Permeability Field

6.1. The 3D Conformal Isotropic Metric

For an isotropic but spatially variable permeability $k(\mathbf{x})$, the permeability tensor is $\mathbf{K} = k(\mathbf{x})\,\mathbf{I}_3$ and the metric is

This is a conformally flat metric in $\mathbb{R}^3$: it is related to the Euclidean metric by the scalar conformal factor $1/k(\mathbf{x})$. In index notation, $g_{ij} = (1/k)\,\delta_{ij}$ and $g^{ij} = k\,\delta^{ij}$, where $\delta_{ij}$ is the Euclidean metric tensor (components: 1 on the diagonal, 0 off it, in every orthonormal chart). Writing $g = e^{2\varphi}\delta$, the conformal factor is $e^{2\varphi} = 1/k$ and hence $\varphi = -\tfrac{1}{2}\ln k$.

The general formula for the scalar curvature of a conformally flat metric $g = e^{2\varphi}\delta$ on $\mathbb{R}^n$ is

where $\Delta$ and $\nabla$ are the ordinary Euclidean operators. For $n = 2$ this gives the Liouville formula $R = -2e^{-2\varphi}\Delta\varphi$. For $n = 3$, the formula acquires a second term.

6.2. The 3D Scalar Curvature Formula

Proof. Set $\varphi = -\tfrac{1}{2}\ln k$ and compute the ingredients:

For the Euclidean Laplacian of $\varphi$: