Papers
Topics
Authors
Recent
Search
2000 character limit reached

Inviscid Incompressible Porous Medium Equations

Updated 7 July 2026
  • IPM equations are active-scalar systems where a transported scalar determines an incompressible velocity via a nonlocal Darcy law.
  • They are analyzed through Eulerian and Lagrangian formulations that reveal stratified steady states, asymptotic stability, and diverse normalization schemes.
  • The study encompasses weak solution constructions via convex integration, critical regularity challenges, and finite-time singularity mechanisms.

The inviscid incompressible porous medium equations (IPM) are a class of active-scalar systems in which a transported scalar—typically denoted ρ\rho or θ\theta—determines an incompressible Darcy velocity through a nonlocal elliptic law. In two dimensions, a standard normalization is

tρ+uρ=0,u=p(0,ρ),u=0,\partial_t \rho + u\cdot \nabla \rho = 0,\qquad u = -\nabla p -(0,\rho),\qquad \nabla\cdot u=0,

while equivalent sign conventions and rotated formulations, such as u+p=(0,ρ)u+\nabla p=(0,\rho) or u+p=(ρ,0)u+\nabla p=(-\rho,0), also occur in the literature. The subject includes Eulerian and Lagrangian formulations, stratified steady states, asymptotic stability and instability mechanisms, convex-integration weak solutions, and several singularity scenarios ranging from reduced one-dimensional models to finite-time blow-up in wedge geometries (Constantin et al., 2014, Elgindi, 2014, Jr, 2011, Bianchini et al., 2024, Dembski, 3 Nov 2025).

1. Governing equations and equivalent formulations

In the whole-plane formulation used for the 2D IPM equation, the scalar satisfies a pure transport law and the velocity is determined by Darcy’s law plus incompressibility: tρ+(u)ρ=0,u=π(0,ρ),u=0.\partial_t \rho + (u\cdot\nabla)\rho = 0,\qquad u = -\nabla \pi -(0,\rho),\qquad \nabla\cdot u = 0. Taking divergence gives

Δπ=x2ρ,-\Delta \pi = \partial_{x_2}\rho,

so the velocity becomes a zero-order singular integral of ρ\rho. One convenient representation is

u=(Δ)1R1ρ,u=\nabla^\perp(-\Delta)^{-1}R_1\rho,

with =(x2,x1)\nabla^\perp=(-\partial_{x_2},\partial_{x_1}) and θ\theta0 the first Riesz transform (Xie et al., 2024).

The same equation is often written with a transported scalar θ\theta1 and constitutive relation

θ\theta2

posed on θ\theta3. In that normalization, the vorticity satisfies

θ\theta4

which makes IPM structurally comparable to 2D Euler and SQG, but with a scalar-to-velocity law tied to Darcy forcing rather than a stream-function law (Constantin et al., 2014).

The same active-scalar structure is studied on θ\theta5, on the strip θ\theta6 with no-penetration, and on wedge domains. In the wedge setting one frequently uses the rotated normalization

θ\theta7

which is equivalent up to orientation of gravity. The stream-function form then becomes

θ\theta8

This suggests that much of the IPM literature is invariant under simple sign and rotation conventions, while the analytic questions depend more strongly on geometry, regularity class, and the structure of the scalar field than on a particular normalization (Dembski, 3 Nov 2025).

2. Lagrangian formulation and analytic trajectories

For classical IPM solutions, the Lagrangian flow map θ\theta9 is defined by

tρ+uρ=0,u=p(0,ρ),u=0,\partial_t \rho + u\cdot \nabla \rho = 0,\qquad u = -\nabla p -(0,\rho),\qquad \nabla\cdot u=0,0

Because tρ+uρ=0,u=p(0,ρ),u=0,\partial_t \rho + u\cdot \nabla \rho = 0,\qquad u = -\nabla p -(0,\rho),\qquad \nabla\cdot u=0,1 is transported, one has

tρ+uρ=0,u=p(0,ρ),u=0,\partial_t \rho + u\cdot \nabla \rho = 0,\qquad u = -\nabla p -(0,\rho),\qquad \nabla\cdot u=0,2

Using the vorticity identity tρ+uρ=0,u=p(0,ρ),u=0,\partial_t \rho + u\cdot \nabla \rho = 0,\qquad u = -\nabla p -(0,\rho),\qquad \nabla\cdot u=0,3, the 2D Biot–Savart law, and incompressibility of the flow map, IPM admits a closed Lagrangian singular-integral system: tρ+uρ=0,u=p(0,ρ),u=0,\partial_t \rho + u\cdot \nabla \rho = 0,\qquad u = -\nabla p -(0,\rho),\qquad \nabla\cdot u=0,4 where tρ+uρ=0,u=p(0,ρ),u=0,\partial_t \rho + u\cdot \nabla \rho = 0,\qquad u = -\nabla p -(0,\rho),\qquad \nabla\cdot u=0,5 is the Poisson bracket. In this representation the source term is a conserved Lagrangian quantity built from the initial scalar and the second component of the flow map (Constantin et al., 2014).

The same work derives an evolution equation for tρ+uρ=0,u=p(0,ρ),u=0,\partial_t \rho + u\cdot \nabla \rho = 0,\qquad u = -\nabla p -(0,\rho),\qquad \nabla\cdot u=0,6 involving a Calderón–Zygmund matrix kernel

tρ+uρ=0,u=p(0,ρ),u=0,\partial_t \rho + u\cdot \nabla \rho = 0,\qquad u = -\nabla p -(0,\rho),\qquad \nabla\cdot u=0,7

analytic away from the origin and of mean zero on spheres. Those structural properties are the key input in the time-analyticity argument.

Under the local well-posedness hypothesis recalled there—specifically tρ+uρ=0,u=p(0,ρ),u=0,\partial_t \rho + u\cdot \nabla \rho = 0,\qquad u = -\nabla p -(0,\rho),\qquad \nabla\cdot u=0,8 with tρ+uρ=0,u=p(0,ρ),u=0,\partial_t \rho + u\cdot \nabla \rho = 0,\qquad u = -\nabla p -(0,\rho),\qquad \nabla\cdot u=0,9, so that u+p=(0,ρ)u+\nabla p=(0,\rho)0—the Lagrangian trajectories are real analytic in time as long as the chord-arc parameter

u+p=(0,ρ)u+\nabla p=(0,\rho)1

remains finite. Thus, for 2D IPM in the classical regime, Eulerian regularity at the level of Hölder-continuous velocity gradient implies real-analytic temporal regularity of particle paths (Constantin et al., 2014).

3. Stratified equilibria and asymptotic stability

A basic structural fact is that sufficiently regular stationary IPM solutions in a bounded simply connected domain with u+p=(0,ρ)u+\nabla p=(0,\rho)2 are necessarily stratified: u+p=(0,ρ)u+\nabla p=(0,\rho)3 The proof uses the stationary transport equation u+p=(0,ρ)u+\nabla p=(0,\rho)4, incompressibility, and the Darcy law to show first that

u+p=(0,ρ)u+\nabla p=(0,\rho)5

and then that the same quantity equals u+p=(0,ρ)u+\nabla p=(0,\rho)6, forcing u+p=(0,ρ)u+\nabla p=(0,\rho)7. This identifies vertical stratification as the canonical steady-state class for IPM (Elgindi, 2014).

Around the linear stratification u+p=(0,ρ)u+\nabla p=(0,\rho)8, the perturbation u+p=(0,ρ)u+\nabla p=(0,\rho)9 satisfies a partially damped dynamics. On u+p=(ρ,0)u+\nabla p=(-\rho,0)0, one obtains the relaxation identity

u+p=(ρ,0)u+\nabla p=(-\rho,0)1

which shows that the deviation from the background profile is a Lyapunov functional and that the velocity is dissipated by stratification rather than by viscosity. The same paper proves global asymptotic stability for sufficiently small perturbations of u+p=(ρ,0)u+\nabla p=(-\rho,0)2: on u+p=(ρ,0)u+\nabla p=(-\rho,0)3, small data in u+p=(ρ,0)u+\nabla p=(-\rho,0)4 with u+p=(ρ,0)u+\nabla p=(-\rho,0)5 generate global solutions converging back to u+p=(ρ,0)u+\nabla p=(-\rho,0)6, while on u+p=(ρ,0)u+\nabla p=(-\rho,0)7 small data in u+p=(ρ,0)u+\nabla p=(-\rho,0)8, u+p=(ρ,0)u+\nabla p=(-\rho,0)9, converge to a stratified stationary limit, because purely vertical modes are undamped and persist (Elgindi, 2014).

A quantitative refinement is available for the three fundamental domains tρ+(u)ρ=0,u=π(0,ρ),u=0.\partial_t \rho + (u\cdot\nabla)\rho = 0,\qquad u = -\nabla \pi -(0,\rho),\qquad \nabla\cdot u = 0.0, tρ+(u)ρ=0,u=π(0,ρ),u=0.\partial_t \rho + (u\cdot\nabla)\rho = 0,\qquad u = -\nabla \pi -(0,\rho),\qquad \nabla\cdot u = 0.1, and tρ+(u)ρ=0,u=π(0,ρ),u=0.\partial_t \rho + (u\cdot\nabla)\rho = 0,\qquad u = -\nabla \pi -(0,\rho),\qquad \nabla\cdot u = 0.2. For the linear stratification tρ+(u)ρ=0,u=π(0,ρ),u=0.\partial_t \rho + (u\cdot\nabla)\rho = 0,\qquad u = -\nabla \pi -(0,\rho),\qquad \nabla\cdot u = 0.3, and on tρ+(u)ρ=0,u=π(0,ρ),u=0.\partial_t \rho + (u\cdot\nabla)\rho = 0,\qquad u = -\nabla \pi -(0,\rho),\qquad \nabla\cdot u = 0.4 also for quasi-linearly stratified backgrounds

tρ+(u)ρ=0,u=π(0,ρ),u=0.\partial_t \rho + (u\cdot\nabla)\rho = 0,\qquad u = -\nabla \pi -(0,\rho),\qquad \nabla\cdot u = 0.5

global asymptotic stability holds provided the buoyancy frequency tρ+(u)ρ=0,u=π(0,ρ),u=0.\partial_t \rho + (u\cdot\nabla)\rho = 0,\qquad u = -\nabla \pi -(0,\rho),\qquad \nabla\cdot u = 0.6 dominates the perturbation size and, in the quasi-linear case, the weighted tρ+(u)ρ=0,u=π(0,ρ),u=0.\partial_t \rho + (u\cdot\nabla)\rho = 0,\qquad u = -\nabla \pi -(0,\rho),\qquad \nabla\cdot u = 0.7 size of tρ+(u)ρ=0,u=π(0,ρ),u=0.\partial_t \rho + (u\cdot\nabla)\rho = 0,\qquad u = -\nabla \pi -(0,\rho),\qquad \nabla\cdot u = 0.8. For tρ+(u)ρ=0,u=π(0,ρ),u=0.\partial_t \rho + (u\cdot\nabla)\rho = 0,\qquad u = -\nabla \pi -(0,\rho),\qquad \nabla\cdot u = 0.9, the whole-space result gives

Δπ=x2ρ,-\Delta \pi = \partial_{x_2}\rho,0

and Δπ=x2ρ,-\Delta \pi = \partial_{x_2}\rho,1 as Δπ=x2ρ,-\Delta \pi = \partial_{x_2}\rho,2. On Δπ=x2ρ,-\Delta \pi = \partial_{x_2}\rho,3 and Δπ=x2ρ,-\Delta \pi = \partial_{x_2}\rho,4, the solution converges not to Δπ=x2ρ,-\Delta \pi = \partial_{x_2}\rho,5 itself but to an asymptotic vertical profile Δπ=x2ρ,-\Delta \pi = \partial_{x_2}\rho,6, with sharp decay rates for the decaying Δπ=x2ρ,-\Delta \pi = \partial_{x_2}\rho,7-dependent modes. The argument uses an anisotropic commutator estimate that lowers the regularity threshold to any real Δπ=x2ρ,-\Delta \pi = \partial_{x_2}\rho,8 in the purely linear case (Jo et al., 2022).

4. Weak solutions, relaxation, and Muskat mixing

For weak-solution theory, IPM is often rewritten as a differential inclusion. In the constant-viscosity case one introduces

Δπ=x2ρ,-\Delta \pi = \partial_{x_2}\rho,9

and a flux variable ρ\rho0, so that the system becomes a linear conservation law

ρ\rho1

together with the pointwise constitutive constraint

ρ\rho2

The corresponding wave cone is

ρ\rho3

and the exact relaxation of the constitutive set is

ρ\rho4

This explicit relaxed set replaces earlier ρ\rho5-configuration arguments and is the basis of a convex-integration construction of infinitely many weak solutions, including nontrivial weak solutions with compact support in time and mixing solutions for the unstable Muskat problem with flat interface. In that setting, the coarse-grained density selected by Otto’s relaxed formulation corresponds to the maximally mixing subsolution (Jr, 2011).

The same h-principle extends to viscosity jump ρ\rho6. In that case weak solutions in ρ\rho7 are again recovered by convex integration once a subsolution is available; nontrivial weak solutions with compact support in time and mixing solutions to the unstable Muskat problem with initial flat interface persist. The viscosity contrast changes the geometry of the relaxation set: the paper identifies a pinch singularity preventing the two fluids from mixing wherever there is neither Rayleigh–Taylor nor vorticity at the interface. It also verifies that the connection between subsolutions and Otto’s Lagrangian relaxed solution, previously established for ρ\rho8, remains valid for ρ\rho9 (Mengual, 2020).

This weak-solution theory is conceptually distinct from homogenization of ideal flow through a perforated solid region. For 2D Euler in a perforated domain, the macroscopic limit is either the full-plane Euler flow or Euler in the exterior of an effective impermeable obstacle, depending on the scaling of hole size and distance; no Darcy-type or classical IPM law appears in the limit. That distinction separates active-scalar IPM from inviscid homogenization through porous geometries (Lacave et al., 2014).

5. Critical regularity, instability, and small-scale formation

A different branch of the literature studies how much of the stratified stability picture survives at critical or near-critical regularity. One mechanism is long-time small-scale formation. For smooth IPM solutions on u=(Δ)1R1ρ,u=\nabla^\perp(-\Delta)^{-1}R_1\rho,0, u=(Δ)1R1ρ,u=\nabla^\perp(-\Delta)^{-1}R_1\rho,1, and the strip, the potential energy

u=(Δ)1R1ρ,u=\nabla^\perp(-\Delta)^{-1}R_1\rho,2

satisfies

u=(Δ)1R1ρ,u=\nabla^\perp(-\Delta)^{-1}R_1\rho,3

This monotonicity implies u=(Δ)1R1ρ,u=\nabla^\perp(-\Delta)^{-1}R_1\rho,4. For several geometric classes of initial data—odd-in-u=(Δ)1R1ρ,u=\nabla^\perp(-\Delta)^{-1}R_1\rho,5 data in u=(Δ)1R1ρ,u=\nabla^\perp(-\Delta)^{-1}R_1\rho,6, symmetric torus data, “bubble” data, and “layered” rearrangements of stratified states—the same paper shows that any global smooth solution must develop unbounded Sobolev norms as u=(Δ)1R1ρ,u=\nabla^\perp(-\Delta)^{-1}R_1\rho,7. These constructions yield nonlinear instability of a large class of stratified steady states in the sense of infinite-time derivative growth, even when global smoothness is assumed (Kiselev et al., 2021).

At the critical Lipschitz level, IPM is mildly ill-posed in u=(Δ)1R1ρ,u=\nabla^\perp(-\Delta)^{-1}R_1\rho,8 near arbitrary vertically stratified backgrounds u=(Δ)1R1ρ,u=\nabla^\perp(-\Delta)^{-1}R_1\rho,9. Writing =(x2,x1)\nabla^\perp=(-\partial_{x_2},\partial_{x_1})0, the perturbation equation contains the variable-coefficient term =(x2,x1)\nabla^\perp=(-\partial_{x_2},\partial_{x_1})1. For every sufficiently small =(x2,x1)\nabla^\perp=(-\partial_{x_2},\partial_{x_1})2, and for any =(x2,x1)\nabla^\perp=(-\partial_{x_2},\partial_{x_1})3 with =(x2,x1)\nabla^\perp=(-\partial_{x_2},\partial_{x_1})4 and =(x2,x1)\nabla^\perp=(-\partial_{x_2},\partial_{x_1})5, there exists initial data =(x2,x1)\nabla^\perp=(-\partial_{x_2},\partial_{x_1})6 with

=(x2,x1)\nabla^\perp=(-\partial_{x_2},\partial_{x_1})7

such that the local solution satisfies

=(x2,x1)\nabla^\perp=(-\partial_{x_2},\partial_{x_1})8

for universal =(x2,x1)\nabla^\perp=(-\partial_{x_2},\partial_{x_1})9. The construction uses smooth compactly supported perturbations θ\theta00 for which the Riesz term θ\theta01 is large while the θ\theta02 norm is small, producing norm inflation in arbitrarily short time. Notably, this applies even when θ\theta03, a regime usually regarded as physically stable (Xie et al., 2024).

At the θ\theta04-critical Sobolev threshold θ\theta05, the stable equation near θ\theta06 is strongly ill-posed. For arbitrarily small θ\theta07, there exist perturbations θ\theta08 with

θ\theta09

such that the corresponding classical solution of IPM with initial data θ\theta10 satisfies

θ\theta11

while remaining bounded in θ\theta12 for every θ\theta13. The mechanism combines a small θ\theta14 perturbation that cancels the background profile near the origin with a multiscale construction generating strong hyperbolic deformation and explosive growth of second derivatives (Bianchini et al., 2024).

Several recent results isolate singularity mechanisms in symmetry classes or reduced models. One infinite-energy class of Castro–Córdoba–Gancedo–Orive type solutions uses the ansatz

θ\theta15

which reduces 2D IPM to a one-dimensional nonlocal equation for θ\theta16: θ\theta17 This reduced equation has the explicit self-similar blow-up solution

θ\theta18

blowing up at θ\theta19. The stability theory identifies a sharp regularity threshold: small θ\theta20 perturbations preserve finite-time blow-up with profile θ\theta21, whereas arbitrarily small θ\theta22 perturbations can destroy asymptotic self-similarity. Via a nonlinear change of variables, this becomes asymptotic stability of the steady family θ\theta23 for the Proudman–Johnson equation (Collot et al., 23 Jul 2025).

Finite-time singularity formation has also been established for the full 2D IPM in wedge domains, without boundary mass. Using 1-homogeneous solutions

θ\theta24

the equation reduces to a one-dimensional angular system

θ\theta25

with θ\theta26. A self-similar ansatz θ\theta27, θ\theta28 yields finite-time gradient blow-up at θ\theta29. The resulting 2D solutions are Lipschitz initially, can be chosen compactly supported, vanish on the boundary, are smooth away from the origin, and satisfy

θ\theta30

The vanishing of density on the boundary is a defining feature: the blow-up mechanism overcomes the full regularizing effect of transport rather than relying on boundary-pinned mass (Dembski, 3 Nov 2025).

A different reduced scenario starts from the 2D periodic half-plane and extracts a boundary-layer model for the boundary trace. The resulting one-dimensional periodic equation is

θ\theta31

where

θ\theta32

This model is locally well-posed for smooth periodic data and blows up in finite time for smooth, bounded, nonnegative, even data with θ\theta33 and θ\theta34 on θ\theta35. The proof uses a Beale–Kato–Majda-type criterion together with a weighted integral inequality for θ\theta36, and positions the model as a boundary-layer analogue of the Córdoba–Córdoba–Fontelos equation (Kiselev et al., 2024).

Broader porous-medium terminology should be handled carefully. A geometric variational theory for an incompressible fluid moving through an elastic porous medium derives coupled fluid–solid equations and recovers Biot’s wave equations in suitable parameter regimes; that framework is a poromechanical generalization rather than the active-scalar IPM equation. Conversely, homogenization of ideal Euler flow through perforated media yields transparency or impermeability, not an effective Darcy/IPM law. This suggests that “porous medium” in PDE usage covers several mathematically distinct limits, of which active-scalar IPM is only one (Farkhutdinov et al., 2020, Lacave et al., 2014).

Topic to Video (Beta)

No one has generated a video about this topic yet.

Whiteboard

No one has generated a whiteboard explanation for this topic yet.

Follow Topic

Get notified by email when new papers are published related to Inviscid Incompressible Porous Medium Equations (IPM).