Porous-media flow

Porous-media flow is fluid motion through the interconnected voids of a porous medium, such as a sand bed, a rock or a filter. Individual passages can be tiny and tortuous, yet their combined discharge can be predicted without resolving every passage. Darcy law supplies this bulk description when viscous resistance dominates pore-scale inertia; porosity measures how much space stores fluid, while permeability of a porous medium measures how readily that space transmits it.

The basic fully saturated model assumes one Newtonian fluid fills the connected pore space of a stationary solid framework, with incompressible flow, a small pore-scale Reynolds number, and averaging volumes large compared with individual pores but small compared with macroscopic variations. These are assumptions of the model developed below, rather than a claim that Darcy law describes every porous-media flow. Two-phase porous-media flow and poroelastic aquifers require additional constitutive and storage information.

For a first calculation, see Darcy flow through a homogeneous slab: a pressure drop determines the discharge, but porosity determines the corresponding pore velocity. Layered porous permeability then shows how geometry alone can make the bulk response anisotropic.

An aquifer is a permeable geological formation containing and transmitting groundwater. A confined aquifer lies beneath a confining layer; an unconfined aquifer instead has a free water table.

An unconfined aquifer has a free groundwater table. In a shallow aquifer of thickness $h(x,t)$, hydrostatic pressure and Darcy's law give horizontal flux proportional to $-h{h}_x$.

For pore fraction $n$, permeability $k$, liquid density $\rho$, viscosity $\mu$, rainfall recharge $R$, and a shallow groundwater thickness $h$,This is a nonlinear diffusion equation with diffusivity proportional to the local thickness.

For $h_t=\kappa (h{h}_x{)}_x$ with a drained boundary $h(0,t)=0$ and a finite current vanishing at its nose, integration by parts gives $\dot{D}=-\kappa [h^2{]}_0^{\infty }/2=0$. The outlet can have finite volume flux because its contribution to the first-moment flux is multiplied by $x=0$. A positive prescribed boundary height instead gives $\dot{D}=\kappa h(0,t{)}^2/2$, so drainage is essential to this invariant.

The conserved first moment of a draining porous current fixes $\alpha +2\beta =0$ in a self-similar solution $h=a{t}^{\alpha }f(x/(d{t}^{\beta }))$. The Boussinesq equation for an unconfined aquifer supplies $\alpha =2\beta -1$, giving $\alpha =-1/2$, $\beta =1/4$. Normalize $\kappa a=d^2$ and the nose at $\eta =1$. Then $-f/2-\eta f^{\prime }/4=(f{f}^{\prime }{)}^{\prime }$ is solved byThe square-root outlet has finite drainage flux and the nose has zero flux. This is a particular similarity solution and its scaling family, not an exact description of arbitrary initial data.

For the dipole similarity solution of a draining porous current, $V=\varphi \int hdx=\varphi ad{t}^{-1/4}/18$. Thus $\dot{V}/V=-1/(4t)$. A shifted similarity origin replaces $t$ by $t+t_0$. Neither a positive rate nor a rate $-1/(2t)$ is compatible with this two-dimensional first-moment-conserving profile.

After recharge stops on a finite aquifer connected to a river, the late solution of the Boussinesq equation for an unconfined aquifer has the form $h\sim t^{-1}F(x/L)$. The profile solves a nonlinear eigenvalue problem, and river discharge decays as $t^{-2}$.

The Dupuit approximation treats flow in a shallow unconfined aquifer as predominantly horizontal, with hydrostatic pressure $p=\rho g(h-z)$ and horizontal pressure gradient independent of depth. Integrating Darcy's law over the saturated depth then gives a volume flux per unit width depending only on the groundwater height and its slope.

For saturated Darcy flow through a long dam occupying a channel of width $\alpha z^2$, let hydraulic conductivity be $K=k\rho g/\mu$ and free-surface height $h(x)$. The Dupuit approximation gives $Q=-K\alpha h^3{h}_x/3$. Fixed heads $h(0)=H$ and $h(L)=\eta$ give $h^4=H^4-(H^4-{\eta }^4)x/L$ and the displayed steady discharge. Porosity relates pore velocity to superficial Darcy velocity but does not independently multiply this discharge when $k$ is intrinsic permeability. For depth-dependent permeability, integrate the horizontal conductivity over the saturated cross-section first.

With permeability of a porous medium $k(z)=k_0(1+\beta z)$, the Dupuit approximation givesHere $\varphi$ is porosity and $R$ is recharge. The low-height limit has mobility proportional to $h$, whereas large heights have mobility proportional to $h^2$. This changes both the forced filling similarity for power-law diffusion and the separable draining profile for power-law diffusion.

Even when most of an unconfined aquifer with depth-dependent permeability is deep, its height vanishes at an absorbing outlet. Near that outlet, finite outward volume flux per unit width $Q$ gives $u_{b}h{h}_x\simeq Q$ and hence $h\simeq (2Qx/u_b{)}^{1/2}$. The high-depth outer approximation instead has $h\simeq [6Qx/(u_b\beta ){]}^{1/3}$. Their crossover lies at height of order $1/\beta$ and distance of order $u_b/({\beta }^{2}Q)$; the thin inner boundary layer transmits the same leading discharge.

A poroelastic aquifer couples pore-fluid storage and flow to elastic compaction of its solid skeleton. Changes in gravity, surface loading, or pore pressure can therefore change groundwater discharge even without changing rainfall.

For a skeleton with Young modulus $E$, undeformed solid fraction ${\varphi }_0$, grain-fluid density difference ${\rho }_s-\rho$, and reference gravity $g_0$, the hydrostatic poroelastic compaction length isIt is the depth over which self-weight changes the solid fraction appreciably.

For a fixed pressure difference across two noncommunicating layers, an activated viscous polymer segment of length ${\ell }_i$ adds hydraulic resistance in series with the remaining fluid. The Darcy flux is the displayed expression when polymer raises dynamic viscosity by a factor $\lambda$. Tracking particle entry time determines activation and exit, while ${\dot{X}}_i=q_i/\varphi$ gives each displacement front. Preferentially occupying the more permeable layer shifts the relative inflow toward the other layer, but cannot increase that layer's absolute flux at fixed pressure. It can delay its breakthrough and leave original reservoir fluid flowing after the unmodified reservoir would be fully swept.

A sharp advected temperature transition in a fluid-filled porous medium. If fluid and solid equilibrate locally in temperature, their effective heat capacity per total volume is $C_{\mathrm{eff}}=\varphi C_f+(1-\varphi )C_s$, while the advected heat flux is $C_{f}UT$ with $U$ the Darcy velocity. Neglecting conduction gives $C_{\mathrm{eff}}{T}_t+C_{f}U{T}_x=0$, so the thermal front speed is $\Gamma U$ with $\Gamma =C_f/C_{\mathrm{eff}}$. The material pore velocity is $U/\varphi$, so thermal retardation relative to that speed is $\Gamma \varphi$. Temperature-dependent viscosity couples thermal-front perturbations back to the Darcy flux.

A porous medium is a solid framework containing connected voids in which a fluid can reside and flow. Its porosity is the void fraction; its permeability of a porous medium measures the ease of pressure-driven flow through the connected pores. The averaged Darcy flux is measured per total cross-sectional area, whereas the pore velocity is that flux divided by porosity. Not every material with voids has connected permeable pores.

Without transverse exchange, a particle chosen uniformly from pore volume samples each streamline's transit time with resident-volume weight. This differs from the flux-weighted residence time in a porous layer, which samples at the inlet proportionally to $v$. For constant porosity and a linear positive velocity $v=v_0+(v_1-v_0)y/H$, the two means are $L\log (v_1/v_0)/(v_1-v_0)$ and $2L/(v_0+v_1)$, respectively. Their limits agree for uniform velocity.

A wetting liquid penetrates pore space under a pressure difference that includes capillary suction. With a planar front, fixed pressure drop and negligible gravity/inertia, Darcy's law and pore-volume storage give $z\propto t^{1/2}$. The front pressure must be defined relative to the gas phase consistently; adding a positive suction to the inlet pressure is different from prescribing a positive liquid front pressure.

Radial capillary imbibition into a half-space has hemispherical flow area $2\pi r^2$. Its radial hydraulic resistance is proportional to $1/R_s-1/R$, so its inlet volume flux approaches $2\pi k\Delta P{R}_s/\mu$. The growing front area gives late-time $R\propto t^{1/3}$, unlike the planar $t^{1/2}$ law. Both results assume capillary-dominated geometry with negligible directional gravity effects.

When a hemispherical wetted front loses liquid at fixed volumetric flux per unit area $F_e$, equilibrium balances inlet flow with $2\pi R^2{F}_e$. The displayed radius exceeds $R_s$ and is stable because inlet flow decreases while evaporating area increases with radius. Porosity affects the adjustment time but cancels from the equilibrium radius.

Immiscible fluids share pore space, with each phase's Darcy flux determined by its pressure, viscosity and relative permeability. Fluid saturation specifies how much pore space each phase occupies. A capillary pressure couples their pressures, while their individual mass conservation laws determine displacement and saturation evolution.

When a nonwetting phase recedes from pores, surface tension and pore geometry can retain a residual saturation $s$. If $h$ is the mobile invaded thickness and $M$ its maximum past value, stored phase volume per plan area is $\varphi [(1-s)h+sM]$. Newly invading pores have $M=h$; during recession $M$ remains fixed. With mobile Darcy velocity $u$, conservation gives advance speed $u/\varphi$ and recession speed $u/[\varphi (1-s)]$. The latter is faster for $0<s<1$.

When a thinning porous gravity current leaves fraction $s$ behind, receding storage changes by only $\varphi (1-s)dh$ after accounting for immobile fluid. Thus its advective thinning equation is $\varphi (1-s)h_t+u{h}_x=0$. Advancing into previously unwetted pores still requires the full storage $\varphi h$, so the nose moves at $u/\varphi$. With exponentially decreasing injection the mobile volume per unit well length is $V=Q\tau (e^{-st/\tau }-e^{-t/\tau })/(1-s)$. It obeys $V^{\prime }=Q{e}^{-t/\tau }-sV/\tau$; the second term becomes immobile retained volume.

For residual fraction $0<s<1$, differentiation of $V=Q\tau (e^{-st/\tau }-e^{-t/\tau })/(1-s)$ gives its unique positive maximum at the displayed time. The maximum volume is $Q\tau s^{s/(1-s)}$. When $s=0$ the mobile volume approaches $Q\tau$ monotonically, with no finite maximum. The limiting $s\to 1$ expression is $V=Qt{e}^{-t/\tau }$ with peak at $t=\tau$.

For capillary residual trapping of an initial triangular current, the method of characteristics translates its leading and trailing faces at $v_A=u/\varphi$ and $v_R=u/[\varphi (1-s)]$. Their intersection has height $a[L-(v_R-v_A)t/2]$. Extinction occurs at $t_{\ast }=2L/(v_R-v_A)$ and position $x_{\ast }=2L/s$. These formulas require $0<s<1$; at zero residual saturation the triangle translates indefinitely.

For the triangular current with capillary retention, the initial trailing face is already the maximum for $0<x<L$. Further upslope, the maximum occurs when the moving crest passes. HenceThe final trapped phase has saturation $s$ below this envelope. Its volume is $\varphi s\int M_{\infty }dx=\varphi a{L}^2$, exactly the initial mobile volume. Maximum invaded thickness is distinct from residual saturation: $s$ multiplies occupied pore volume rather than shrinking the geometrical envelope.

Injecting water to displace oil through a porous reservoir. Spatial permeability contrasts cause preferential flow and early water breakthrough, while unswept slower regions retain oil. The recovery also depends on fractional flow, capillary pressure and phase mixing; a passive streamline model alone does not capture all two-phase effects.

The one-dimensional mass conservation equation for saturation when total Darcy flux is constant and capillary pressure and gravity are neglected. It is a scalar conservation law with dimensional characteristic speed $(Q/\varphi )F^{\prime }(s)$. Shock speeds obey a Rankine-Hugoniot condition; physically admissible solutions satisfy an entropy selection. Capillary diffusion regularizes the saturation transition.

Eliminating phase pressures with $p_c=p_n-p_w$ gives $\varphi s_t+QF(s{)}_x=(D_c{s}_x{)}_x$. For decreasing capillary pressure, $D_c\ge 0$. This nonlinear diffusion broadens a saturation shock; it can degenerate where a phase mobility vanishes. Endpoint conservation retains the limiting Rankine-Hugoniot condition.

When a rarefaction wave attaches to a saturation shock, its terminal characteristic speed must match the Rankine-Hugoniot condition speed. This gives a tangent from the initial state $(s_0,F(s_0))$ to the fractional-flow curve at $s_s$. It selects the upstream shock saturation; the actual dimensional shock speed still includes total flux divided by porosity.

The fraction of total advective Darcy flux carried by one phase when capillary and gravity corrections are omitted. For two phases it is ${\lambda }_w/({\lambda }_w+{\lambda }_n)$. This quantity lies between zero and one for nonnegative phase mobilities. Its nonlinear dependence on fluid saturation determines characteristic speeds and shocks in the Buckley-Leverett equation.

The fraction of pore volume occupied by a fluid phase. For a completely fluid-filled two-phase pore space the saturations sum to one. Residual saturation refers to fluid remaining immobile under the modeled displacement; it can impose zero relative mobility at the endpoint.

For steady incompressible flow with no stagnant volume, partition the layer into streamtubes of volume flux $dQ$. The pore volume in a streamtube is $\tau dQ$, so $\int \tau dQ=V_{\mathrm{pore}}$. This gives the mean residence time $V_{\mathrm{pore}}/Q$ even when individual path times vary widely. Streamfunction labels measure flux in two-dimensional flow, so they supply the correct weighting.

A liquid layer of depth $h$ above a dry deep porous substrate drives a Darcy flux $w=(gk/\nu )(1+h/l)$ through its saturated depth $l$, assuming negligible capillary pressure and negligible resistance from the displaced phase. The pressure head contributes $h/l$ and gravity contributes one. Volume conservation gives $\varphi \dot{l}=w$. For a finite column, $h+\varphi l=h_0$ and initiallyThe singular initial flux is an idealization; the continuum model requires penetration beyond the pore scale.

Darcy law is a linear relation between the bulk Darcy flux and the net driving force per unit volume. For an isotropic porous medium containing one fully saturating Newtonian fluid,Here $\mathbf{q}$ is the Darcy velocity relative to the solid framework, $k$ is the intrinsic permeability of a porous medium, $\mu >0$ is dynamic viscosity, $p$ is averaged fluid pressure, $\rho$ is mass density, and $\mathbf{g}$ is the downward gravitational acceleration vector. The pressure gradient exerts force density $-\nabla p$, while gravity contributes $\rho \mathbf{g}$; the minus sign makes the discharge follow their sum. The symbols $\mathbf{u}$ and $\mathbf{q}$ are both commonly used for Darcy velocity; neither denotes pore velocity here.

With upward coordinate $z$ and positive gravitational acceleration magnitude $g$, $\mathbf{g}=-g{\mathbf{e}}_z$. Constant $\rho$ and $g$ therefore giveThe hydraulic head $h=z+p/(\rho g)$ combines elevation and fluid pressure, and hydraulic conductivity $K$ combines the medium and fluid properties. Hydrostatic pressure satisfies $\nabla p=\rho \mathbf{g}$, so it produces zero Darcy flux even though fluid pressure increases downward. Conversely, equal pressures at the two ends of a vertical column do not imply zero flow: their elevations give different hydraulic heads.

The appeal of Darcy law is that a complicated connected pore geometry enters a macroscopic calculation through a measurable coefficient. It predicts filter discharge and aquifer flow; combined with mass conservation, it gives an equation for the evolving fluid pressure. Its simple linear form requires a small pore-scale Reynolds number, slow adjustment compared with pore-scale inertial relaxation, and a valid continuum approximation. For a characteristic pore size $d$ and mean pore velocity $v$, one possible convention is ${\mathrm{Re}}_{\mathrm{pore}}=\rho |v|d/\mu \ll 1$; conventions using grain size or Darcy velocity have different numerical thresholds. Merely calling flow laminar does not guarantee negligible inertia. Unsaturated pores, several fluid phases, appreciable inertia, or a fluid whose stress is not Newtonian require further modeling. Woessner and Poeter discuss these limits in their authored hydrology text.^[ref]

Henry Darcy's 1856 book Les fontaines publiques de la ville de Dijon was published by Victor Dalmont; the University Grenoble Alpes digital collection verifies this bibliographic record and supplies a digitized copy.^[ref] Woessner and Poeter's account connects Darcy law to experiments on water flowing through sand filters, relating discharge to cross-sectional area, hydraulic head difference and sample length.^[ref] The modern intrinsic permeability of a porous medium notation above should not be read as a quotation of Darcy's original notation.

Hydraulic head is the sum of elevation and pressure head, measured in metres relative to a chosen height datum. For constant mass density $\rho$ and gravitational acceleration magnitude $g$, the pressure head $p/(\rho g)$ is the height of a column of the same fluid supported by fluid pressure $p$. A water-filled open measuring tube records hydraulic head as the elevation of its water level when $p$ is gauge fluid pressure. Changing the common pressure reference or height datum adds a constant to hydraulic head and leaves its driving gradient unchanged.

For water with $\rho =1000\mathrm{kg}{\mathrm{m}}^{-3}$ and $g=9.81\mathrm{m}{\mathrm{s}}^{-2}$, a point at $z=2\mathrm{m}$ with gauge fluid pressure $p=9810\mathrm{Pa}$ has $h=3\mathrm{m}$. A second point at $z=0$ with gauge fluid pressure $29430\mathrm{Pa}$ also has $h=3\mathrm{m}$: these values are compatible with a stationary column under hydrostatic pressure. It is a difference in hydraulic head, rather than a difference in fluid pressure alone, that drives saturated Darcy flow. This form uses constant mass density; one cannot generally pull spatially varying $\rho$ through the gradient.^[ref]

Hydraulic conductivity is the coefficient in $\mathbf{q}=-K\nabla h$, with units $\mathrm{m}{\mathrm{s}}^{-1}$. Because the gradient of hydraulic head is dimensionless, hydraulic conductivity has the dimensions of velocity, but it is not the actual Darcy velocity unless the driving gradient has magnitude one. Intrinsic permeability of a porous medium $k$ has units ${\mathrm{m}}^2$ and describes pore geometry; hydraulic conductivity also depends on mass density, dynamic viscosity and gravitational acceleration.^[ref]

For $k=10^{-12}{\mathrm{m}}^2$, $\rho =1000\mathrm{kg}{\mathrm{m}}^{-3}$, $g=9.81\mathrm{m}{\mathrm{s}}^{-2}$ and $\mu =10^{-3}\mathrm{Pa}\mathrm{s}$, hydraulic conductivity is $K=9.81\times 10^{-6}\mathrm{m}{\mathrm{s}}^{-1}$. Doubling dynamic viscosity while keeping pore geometry, mass density and gravitational acceleration fixed halves hydraulic conductivity without changing intrinsic permeability of a porous medium. Thus a value quoted as a material's hydraulic conductivity needs its fluid conditions.

Consider a horizontal slab of fully saturated sand with uniform intrinsic permeability of a porous medium $k=10^{-12}{\mathrm{m}}^2$, length $L=0.50\mathrm{m}$, total cross-sectional area $A=0.020{\mathrm{m}}^2$, and porosity $\varphi =0.25$. Use one Newtonian fluid with mass density $\rho =1000\mathrm{kg}{\mathrm{m}}^{-3}$ and dynamic viscosity $\mu =10^{-3}\mathrm{Pa}\mathrm{s}$. Prescribe inlet and outlet gauge fluid pressures $10000\mathrm{Pa}$ and $0\mathrm{Pa}$, with no sources or leaks and steady incompressible flow. The two ends have the same elevation, so gravity has no component along the slab. Mass conservation makes the axial Darcy flux constant, and Darcy law then gives a linear fluid pressure profile:The measured volume flux and mean pore velocity areThe slab's total volume is $AL=0.010{\mathrm{m}}^3$, of which $\varphi AL=0.0025{\mathrm{m}}^3=2.5\mathrm{L}$ is pore volume. In the ideal uniform advective model, a nonreacting tracer has transit time $L/v=6250\mathrm{s}\simeq 1.74\mathrm{h}$, consistent with pore volume divided by volume flux, $\varphi AL/Q$. The incorrect substitution of Darcy velocity for pore velocity would give $25000\mathrm{s}$, four times too long. Real passages have different path lengths and local speeds; these numbers describe a bulk mean rather than every particle's journey.

As an independent check, gravitational acceleration magnitude $g=9.81\mathrm{m}{\mathrm{s}}^{-2}$ gives hydraulic conductivity $K=9.81\times 10^{-6}\mathrm{m}{\mathrm{s}}^{-1}$, hydraulic head drop $\Delta h=10000/(1000\times 9.81)\mathrm{m}\simeq 1.01937\mathrm{m}$, and $q=K\Delta h/L=2.0\times 10^{-5}\mathrm{m}{\mathrm{s}}^{-1}$ again. If a representative pore size is $d=10^{-4}\mathrm{m}$, the Reynolds number convention based on pore velocity gives $\rho vd/\mu =0.008$, supporting neglect of pore inertia. This size is an illustrative extra measurement, not something determined by porosity and intrinsic permeability of a porous medium alone.

Regenerate the diagram by running python ../darcy-flux-and-pore-velocity.py from _media; it writes its PNG into the current directory.

Darcy flow is slow flow through a porous medium whose averaged flux obeys the Darcy law. The Darcy velocity is the fluid volume flux per total cross-sectional area; it differs from the velocity averaged only over pore space.

In a thin fluid-filled crack, lubrication theory gives flux per unit depth $q=-h^3{p}_x/(12\mu )$. The surrounding Darcy flow has normal velocity $-k{p}_y/\mu$. Conservation in the crack gives $q_x=-[u_y{]}_{-}^{+}$, yielding the displayed condition. Pressure is continuous across the crack. This reduced model neglects transverse pressure variation and assumes wall slip lengths are small compared with $h$.

For $h=H_0(1-x^2/a^2{)}^{1/6}$, elliptic coordinates make $h^3/\sin \eta$ constant on the upper crack face. Pressure continuity and harmonicity reduce the boundary condition to $F^{\prime }(0)=\alpha F(0)$ for the cosine mode, while uniform far flow gives $F\sim -e^{\xi }/2$. Solving $F^{\prime \prime }=F$ yields the displayed pressure. The midpoint crack flux is $2aU\alpha /(1+\alpha )$; small $\alpha$ leaves the background pressure gradient nearly unchanged, while large $\alpha$ makes the crack nearly equipotential.

Subtracting the uniform-flow pressure $-\mu Ux/k$ leaves $(\mu Ua/k)\alpha e^{-\xi }\cos \eta /(1+\alpha )$. The coordinate map has $r\sim a{e}^{\xi }/2$, so $e^{-\xi }\sim a/(2r)$, giving the displayed dipole coefficient. The one-half is essential. In the perfect-conductor limit the independent expression $p=-(\mu U/k)\mathrm{Re}\sqrt{z^2-a^2}$ has the same far-field expansion.

For dipole pressure $\delta p\sim Cx/(x^2+y^2)$, choose a rectangular control volume whose transverse size tends to infinity faster than its streamwise size. Each end face contributes $\pi C$ to $\int \delta p{n}_{x}ds$, so fixed-throughflow boundary work gives $\Delta \mathcal{D}=-2\pi UC$, the displayed result. This control-surface choice removes the integrated far-field velocity correction. A circular cutoff of the pressure-only functional gives $-\pi UC$ instead and does not implement the same throughflow normalization.

For constant wall heat-flux density $F>0$, effective thermal conductivity $k_e$, buoyancy mobility $B={\rho }_{0}g\alpha K/\mu$ and thermal transport coefficient $\kappa =k_e/({\rho }_0{c}_f)$, a slender Darcy flow boundary layer has width $\delta =(\kappa k_{e}z/(BF){)}^{1/3}$. Writing $\eta =y/\delta$, $\psi =\kappa zf(\eta )/\delta$ and $T-T_a=(F\delta /k_e)f^{\prime }(\eta )$ gives the displayed equation with $f(0)=0$, $f^{\prime \prime }(0)=-1$ and $f^{\prime }(\infty )=0$. The exact integral is ${\int }_0^{\infty }(f^{\prime }{)}^{2}d\eta =1$. Consequently its kinematic Darcy momentum flux is $\int w^{2}dy=B\kappa Fz/k_e$ and its buoyancy flux is $g\alpha \kappa Fz/k_e$. These are per span on one heated face.

Assume an exponentially decreasing temperature profile while imposing wall heat flux, impermeability and the exact integrated heat budget. For $f=C(1-e^{-\eta /d})$, the wall condition gives $C=d^2$ and $\int (f^{\prime }{)}^{2}d\eta =1$ gives $d^3=2$. This yields the displayed integral approximation and discharge $\int wdy\simeq 2^{2/3}\kappa z/\delta$. It satisfies the boundary conditions and exact flux budget, but not the differential equation pointwise; its mass-flux coefficient is approximate.

For two incompressible fluids separated by a planar material front $X(t)$, the Darcy flux is uniform and the pressure drop sums the resistances of the two layers: $U=k\Delta P/[\mu X+{\mu }_h(L-X)]$. The front moves at the pore velocity $X^{\prime }=U/\varphi$. Integrating from $X(0)=0$ yields the displayed relation. The formula holds up to breakthrough $X=L$, with either sign of the viscosity difference; equal viscosities give constant speed.

If a thermal front has speed $Y^{\prime }=\Gamma U$ with $U$ the Darcy velocity, while a material front has speed $X^{\prime }=U/\varphi$, their position ratio from coincident initial locations is $Y/X=\Gamma \varphi$. Calling a retardation factor the ratio to pore velocity instead would give $Y/X=\Gamma$; these conventions must not be mixed. With cold viscosity $\mu$ and heated viscosity $b\mu$, the displacement's effective viscosity is ${\mu }_{\mathrm{eff}}=\mu [b-(b-1)\Gamma \varphi ]$. In the pressure resistance the cold layer has length $Y$, the hot injected layer has length $X-Y$, and the displaced layer has length $L-X$.

A slender planar plume in a saturated porous medium has buoyancy velocity $w=\beta (T-T_0)$, with $\beta ={\rho }_{0}g\alpha \Pi /\mu$. The leading heat equation retains both advection components and lateral thermal conduction: $u{T}_x+w{T}_z=\kappa T_{xx}$. Its convective heat input per span is conserved, $F_h=C_h\int w(T-T_0)dx$. Defining $Q=\beta F_h/C_h$ gives width proportional to $z^{2/3}$ and vertical speed/temperature excess proportional to $z^{-1/3}$.

For a slender porous thermal plume, Darcy law gives $w=\beta (T-T_0)$, with $\beta =Kg\alpha /\nu$. If $C_h={\rho }_0{c}_p$ and $F_h=C_h\int w(T-T_0)dx$ is the conserved excess convective heat flux per span, then the kinematic momentum flux is $M=\int w^{2}dx=\beta F_h/C_h$, and the buoyancy flux is $B=\int g\alpha (T-T_0)wdx=g\alpha F_h/C_h$. The hyperbolic-secant porous plume profile has $Q=(36\kappa Mz{)}^{1/3}$. Multiplication of $Q$ and $M$ by ${\rho }_0$ gives physical mass flux and the density-weighted Darcy second moment. Intrinsic pore-velocity momentum flux instead has an additional inverse-porosity factor. Momentum conservation here follows from the heat/velocity proportionality, not from an inertia balance: porous drag balances buoyancy locally.

The symmetric steady porous thermal plume has $w=W{\mathrm{sech}}^2(x/b)$, $T-T_0=w/\beta$ and $\psi =-Wb\tanh (x/b)$. Its exact normalizations are $W{b}^2=6\kappa z$ and $(4/3)W^{2}b=Q$, so $b=(48{\kappa }^2{z}^2/Q{)}^{1/3}$ and $W=(3Q^2/(32\kappa z){)}^{1/3}$. The streamfunction supplies lateral entrainment required by continuity.

Truncating the steady porous thermal plume profile below a head height and conserving supplied heat gives ${\dot{Z}}_h\simeq \int w^{2}dx/\int wdx=(2/3)W(Z_h)$. Thus $Z_h\propto t^{3/4}$. This is an integral bulk-head estimate, while fastest centreline parcels have characteristic speed $W$; the approximation does not resolve the transient nose.

A Darcy velocity $\mathbf{u}$ driven by pressure in a stationary porous skeleton dissipates mechanical power $\mu |\mathbf{u}{|}^2/\Pi$ per bulk volume. It follows by dotting Darcy law with the discharge: $-\mathbf{u}\cdot \nabla p=\mu |\mathbf{u}{|}^2/\Pi$. Its contribution to a heat equation can be neglected only when it is small relative to the relevant thermal-energy transfer.

Darcy velocity, also called Darcy flux or specific discharge, is fluid volume flux relative to the solid framework per unit total cross-sectional area. For a uniform axial discharge $Q$ through bulk area $A$,The area includes both solid and voids. Hence Darcy velocity is a superficial flux with units $\mathrm{m}{\mathrm{s}}^{-1}$, even though no fluid travels through the solid part. For a nonuniform discharge, $Q={\int }_A\mathbf{q}\cdot \mathbf{n}dA$, with $\mathbf{n}$ the unit normal; $Q/A$ is the area average of the normal Darcy flux.

If porosity is $ϵ$, the pore-volume-averaged liquid velocity is ${\mathbf{v}}_l$, and the local solid framework velocity is ${\mathbf{v}}_s$, the convention used here is $\mathbf{u}=ϵ({\mathbf{v}}_l-{\mathbf{v}}_s)$. In a stationary, fully saturated framework with all pore volume participating, ${\mathbf{v}}_l=\mathbf{q}/ϵ$. The notation $\varphi =ϵ$ for porosity is equally common. Darcy flow through a homogeneous slab gives $q=2.0\times 10^{-5}\mathrm{m}{\mathrm{s}}^{-1}$ but pore velocity $v=8.0\times 10^{-5}\mathrm{m}{\mathrm{s}}^{-1}$ at $\varphi =0.25$. Keeping these quantities separate lets Darcy law predict filter discharge while pore velocity predicts ideal advective travel time.^[ref]

Pore velocity is the fluid velocity averaged over pore volume relative to the solid framework. For uniform porosity $\varphi$, full saturation and all pores participating, $\mathbf{v}=\mathbf{q}/\varphi$, with Darcy velocity $\mathbf{q}$ defined per bulk area. The identity follows from volume averaging: extending the relative fluid velocity by zero in the solid gives a bulk average equal to the pore fraction times the pore-volume average. It does not assert that every local passage has open area $\varphi A$ or identical speed.

In a uniform advective model, travel over length $L$ takes $L/v=\varphi L/q$. If some pore volume is inaccessible to transport, the corresponding estimate uses effective porosity ${\varphi }_e$ and $v=q/{\varphi }_e$. A connected but nearly stagnant passage can exchange solute slowly, so even this mean does not describe all arrival times. Flux-weighted residence time in a porous layer explains one precise averaging convention.^[ref]

Darcy-Bénard convection is thermal convection in a porous layer heated from below, with velocity governed by Darcy's law. With fixed temperatures and impermeable boundaries, a normalized two-dimensional model uses streamfunction $\psi$ and temperature $\theta$:The control parameter is a porous-medium Rayleigh number. On $0<z<1$ with temperature values $0,-1$, the conductive state is $\psi =0$, $\theta =-z$.

A tall porous column with square width $a$, insulated impermeable sidewalls and a maintained upward temperature decrease $\Gamma$ admits nearly vertical, zero-net-flux overturning modes. Let $B={\rho }_{0}g\alpha K/\mu$ be buoyancy mobility and $\kappa =k_e/({\rho }_0{c}_f)$ the thermal transport coefficient accompanying Darcy flux. Neumann transverse modes have $k_{\perp }^2={\pi }^2(m^2+n^2)/a^2$. In the long-column limit the lowest nonconstant modes have growth rate proportional to $B\Gamma -\kappa {\pi }^2/a^2$. They are warm ascending and cold descending columns with return flow at distant ends. The two lowest modes, varying along either side direction, are degenerate; linear theory alone does not select their nonlinear superposition.

In a unit square with prescribed conductive boundary temperatures and impermeable boundaries, the stream function convention $\mathbf{u}=(-{\psi }_z,0,{\psi }_x)$ gives linear equations $\Delta \psi =R{\theta }_x$, ${\theta }_t={\psi }_x+\Delta \theta$, with $\psi =\theta =0$ on all sides. The velocity equation follows by taking the vertical-plane curl of Darcy law, and the temperature equation by linearizing the heat equation about $1-z$. Unlike a laterally unbounded Darcy layer, the conducting side boundaries produce a two-dimensional real marginal eigenspace at $R_c=8{\pi }^2$.

For $q=q_c+ϵq_1$, ${\theta }_t=ϵ\lambda \theta$, the first-order equations are $\Delta {\psi }_1-q_c^2{\theta }_{1x}=2q_c{q}_1{\theta }_{0x}$ and ${\psi }_{1x}+\Delta {\theta }_1=\lambda {\theta }_0$. Multiply them by ${\psi }_0,q_c^2{\theta }_0$ and integrate. Integration by parts and the leading equations cancel the homogeneous correction terms, giving $2q_c{q}_1⟨{\psi }_0{\theta }_{0x}⟩+q_c^2\lambda ⟨{\theta }_0^2⟩=0$. Since $⟨{\psi }_0{\theta }_{0x}⟩=-⟨|\nabla {\theta }_0{|}^2⟩$, the displayed title formula follows. The reciprocal norm ratio is incompatible with this solvability condition. Reflection parity diagonalizes the first-order splitting of the two marginal modes.

At zero growth rate put $P=\psi +iq\theta$, $q=\sqrt{R}$. The coupled Darcy equations become $\Delta P+iq{P}_x=0$. Substituting $P=e^{-iq(x-1/2)/2}F$ gives $\Delta F+(q^2/4)F=0$ with homogeneous Dirichlet boundary conditions. The square Dirichlet Laplacian eigenvalues give $q^2=4{\pi }^2(m^2+n^2)$, $m,n\ge 1$, whose minimum is $8{\pi }^2$. For $S=\sin \pi x\sin \pi z$, $a=q_c/2$, the complex multiples $F=S,iS$ give real pairs $(S\cos (a\xi ),-S\sin (a\xi )/q_c)$ and $(S\sin (a\xi ),S\cos (a\xi )/q_c)$, $\xi =x-1/2$. Their reflection parities are opposite, yielding two independent physical modes.

With permeability $K$, depth $d$, kinematic viscosity $\nu$ and the relevant effective thermal diffusivity ${\kappa }_T$, the Darcy Rayleigh number compares buoyancy-driven Darcy advection to thermal diffusion. It differs from the clear-fluid Rayleigh number by replacing $d^3$ with $Kd$. Effective heat capacity and porosity factors may be absorbed into the definition of ${\kappa }_T$ and the thermal-time scale, and must be kept consistent with the chosen Darcy law convention.

A perturbation proportional to $\sin (n\pi z)e^{st+ikx}$ has growth rateMinimizing the neutral curve $(k^2+n^2{\pi }^2{)}^2/k^2$ gives $n=1$, $k=\pi$ and $R_c=4{\pi }^2$. Thus the first convection rolls have horizontal wavelength twice the layer depth.