Two-Layer Darcy–Forchheimer Model

• The model establishes a coupled framework by integrating distinct porous regions with nonlinear inertial resistance and precise interface transmission conditions.
• It utilizes varying formulations—such as Brinkman–Forchheimer/Darcy and matrix–fracture reductions—to capture diverse flow physics in heterogeneous media.
• Robust numerical techniques, including mixed finite element and multiscale methods, validate convergence and accuracy across applications like thermal-fluid and fracture flow modeling.

A two-layer Darcy–Forchheimer model is a coupled porous-media formulation in which adjacent regions are assigned distinct constitutive laws, dimensional reductions, or transport closures, while remaining linked by explicit transmission conditions. In the canonical porous-flow setting, a more permeable layer or fracture carries nonlinear inertial resistance of Darcy–Forchheimer or Brinkman–Forchheimer type, whereas a less permeable surrounding layer follows Darcy flow; in reduced-order thermo-fluid settings, the term can also denote a thickness-averaged Darcy–Brinkman–Forchheimer fluid layer coupled to a conductive base-solid layer through interfacial heat exchange. Across these variants, the defining ingredients are domain decomposition, layer-specific momentum laws, and interface conditions that enforce mass transfer together with pressure, traction, or heat-flux transmission (Caucao et al., 2023, Knabner et al., 2013, Saito et al., 30 Sep 2025).

1. Canonical configurations

The literature uses the expression two-layer Darcy–Forchheimer model in several closely related but not identical senses. One class of models treats two full-dimensional porous regions with different permeability scales. In the formulation developed for a coupled Brinkman–Forchheimer/Darcy system, a bounded Lipschitz domain $\Omega \subset \mathbb{R}^n$ is decomposed into two non-overlapping subdomains: a more permeable region $\Omega_{BF}$, where Brinkman–Forchheimer applies, and a less permeable region $\Omega_D$, where Darcy applies. The interface is $$\Gamma := \partial \Omega_{BF} \cap \partial \Omega_D$$.

A second class is mixed-dimensional. Knabner and Roberts analyze a matrix–fracture configuration in which $\Omega = \Omega_1 \cup \gamma \cup \Omega_2$, with $\gamma$ a $(d-1)$-dimensional planar fracture separating two matrix subdomains. The matrix obeys Darcy flow, while the fracture supports tangential Darcy–Forchheimer flow and exchanges mass with the matrix through leak-off terms.

A third class arises in reduced-order heat-transfer design. There, the two layers are not two porous-flow subdomains of the same type, but rather an upper thermal-fluid layer governed by a thickness-averaged Darcy–Brinkman–Forchheimer momentum balance and a lower base-solid layer governed by conduction, coupled by an analytically derived interfacial heat-transfer coefficient.

Configuration	Layer laws	Principal coupling
$\Omega_{BF}$–$\Omega_D$ decomposition	Brinkman–Forchheimer / Darcy	Normal-velocity continuity and momentum continuity
Matrix–fracture $\Omega_1 \cup \gamma \cup \Omega_2$	Darcy / Darcy–Forchheimer	Leak-off mass balance and Robin-type pressure jump
Mixed-dimensional fracture reduction	Darcy in bulk / Forchheimer on $\Omega_{BF}$0	Flux exchange and reduced fracture balance
Thermal-fluid/base-solid reduction	Darcy–Brinkman–Forchheimer / heat conduction	Interfacial heat exchange

This multiplicity of usages is important. A common misconception is that two-layer Darcy–Forchheimer necessarily means that both layers satisfy the same nonlinear Darcy–Forchheimer law. In the cited literature, hybrid Darcy/Forchheimer and Brinkman/Forchheimer/Darcy combinations are at least as central as the fully nonlinear two-layer case (Knabner et al., 2013, Arrarás et al., 2018).

2. Governing equations and transmission mechanisms

In the coupled Brinkman–Forchheimer/Darcy setting, the more permeable layer satisfies

$\Omega_{BF}$1

$\Omega_{BF}$2

while the less permeable layer satisfies

$\Omega_{BF}$3

$\Omega_{BF}$4

The interface conditions are continuity of normal velocity,

$\Omega_{BF}$5

and continuity of momentum in the normal direction,

$\Omega_{BF}$6

A Lagrange multiplier $\Omega_{BF}$7 enforces the normal-velocity constraint and represents the interface pressure trace.

In the matrix–fracture model of Knabner and Roberts, each matrix subdomain obeys Darcy’s law

$\Omega_{BF}$8

whereas the fracture supports tangential Darcy–Forchheimer flow,

$\Omega_{BF}$9

$\Omega_D$0

Here $\Omega_D$1 is the leak-off contribution from the matrix. The matrix and fracture are connected by the Robin-type pressure jump law

$\Omega_D$2

with $\Omega_D$3 and $\Omega_D$4.

In the mixed-dimensional reduction used for fractured porous media, the fracture is collapsed to a one-dimensional manifold $\Omega_D$5, and the constitutive law becomes

$\Omega_D$6

with balance

$\Omega_D$7

and interface exchange

$\Omega_D$8

These formulations show that interface transmission is not uniquely prescribed by pressure continuity. Depending on the physical setting, the relevant coupling may be traction continuity, a Robin-type pressure jump, or a reduced leak-off balance (Caucao et al., 2023, Knabner et al., 2013, Arrarás et al., 2018).

3. Variational structure, monotonicity, and limiting relations

The coupled Brinkman–Forchheimer/Darcy model admits a mixed/dual-mixed weak formulation. In $\Omega_D$9, the formulation is standard mixed, with $$\Gamma := \partial \Omega_{BF} \cap \partial \Omega_D$$0 for velocity; in $$\Gamma := \partial \Omega_{BF} \cap \partial \Omega_D$$1, it is dual-mixed, with $$\Gamma := \partial \Omega_{BF} \cap \partial \Omega_D$$2 so that local mass conservation and normal-flux transmission are built into the space. The pressure/multiplier space is

$$\Gamma := \partial \Omega_{BF} \cap \partial \Omega_D$$3

with the multiplier $$\Gamma := \partial \Omega_{BF} \cap \partial \Omega_D$$4 entering through interface terms

$$\Gamma := \partial \Omega_{BF} \cap \partial \Omega_D$$5

and the constraint

$$\Gamma := \partial \Omega_{BF} \cap \partial \Omega_D$$6

The analysis rests on strong monotonicity of the nonlinear operator on the kernel of the constraint and on continuous and discrete inf–sup conditions for the coupled spaces.

In the matrix–fracture setting, Knabner and Roberts formulate the hybrid problem in Banach spaces

$$\Gamma := \partial \Omega_{BF} \cap \partial \Omega_D$$7

and

$$\Gamma := \partial \Omega_{BF} \cap \partial \Omega_D$$8

The fracture nonlinearity enters through

$$\Gamma := \partial \Omega_{BF} \cap \partial \Omega_D$$9

and the interface transmissibility through

$\Omega = \Omega_1 \cup \gamma \cup \Omega_2$0

A recurrent technical inequality is

$\Omega = \Omega_1 \cup \gamma \cup \Omega_2$1

which yields strict monotonicity of the Forchheimer operator and underpins uniqueness.

The hybrid Darcy-matrix/Forchheimer-fracture model is also obtained as a limit of a fully Forchheimer system in which matrix coefficients $\Omega = \Omega_1 \cup \gamma \cup \Omega_2$2 tend to zero. The weak convergences stated in the paper include

$\Omega = \Omega_1 \cup \gamma \cup \Omega_2$3

together with

$\Omega = \Omega_1 \cup \gamma \cup \Omega_2$4

This limit interpretation clarifies that Darcy flow in one layer can be viewed as a vanishing-inertia regime of a more general two-layer Forchheimer model.

A related generalization is available on compact Riemannian manifolds of dimension $\Omega = \Omega_1 \cup \gamma \cup \Omega_2$5 or $\Omega = \Omega_1 \cup \gamma \cup \Omega_2$6. There the transmission problem is stated in $\Omega = \Omega_1 \cup \gamma \cup \Omega_2$7-based Sobolev spaces, with interface conditions

$\Omega = \Omega_1 \cup \gamma \cup \Omega_2$8

$\Omega = \Omega_1 \cup \gamma \cup \Omega_2$9

where $\gamma$0 is strongly positive on tangential fields. Existence and uniqueness for the nonlinear transmission problem are established by a layer-potential representation combined with a fixed-point theorem for sufficiently small data (Caucao et al., 2023, Knabner et al., 2013, Kohr et al., 2016).

4. Discretization strategies and nonlinear solvers

For the coupled Brinkman–Forchheimer/Darcy problem, the finite element discretization uses Bernardi–Raugel elements for the Brinkman–Forchheimer velocity, Raviart–Thomas $\gamma$1 elements for the Darcy velocity, piecewise constants for the pressures, and continuous piecewise linear elements on a coarser trace mesh for the interface multiplier. The discrete spaces are designed so that normal components on $\gamma$2 are evaluated consistently and the discrete lifting operator yields a robust discrete inf–sup condition. Under regularity assumptions

$\gamma$3

the paper derives a Céa/Strang-type estimate and the first-order bound

$\gamma$4

The nonlinear Forchheimer term is treated by Newton-type iterations or Picard linearization. In the “tombstone” test with $\gamma$5, $\gamma$6, $\gamma$7, $\gamma$8, and $\gamma$9, the scheme achieved first-order convergence in all variables; Newton iteration counts were 4 for $(d-1)$0 and up to 9 for $(d-1)$1. In the rectangular heterogeneous example with $(d-1)$2, $(d-1)$3, $(d-1)$4, and $(d-1)$5, Newton counts ranged from 1 for $(d-1)$6 to 8 for $(d-1)$7.

For mixed-dimensional fractured media, a staggered finite-volume discretization produces a nonlinear saddle-point system whose only nonlinear block is the fracture block $(d-1)$8. This structure supports a monolithic full approximation scheme multigrid method with mixed-dimensional restriction, prolongation, and smoothing. The matrix uses a classical five-point Vanka relaxation, while the fracture uses a nonlinear three-point Vanka relaxation. Because the nonlinearity is localized in the fracture, no coupling between matrix and fracture unknowns is needed in the smoother. With W(2,2) cycles and damping $(d-1)$9, numerical experiments reported around 8 iterations to reduce the residual by $\Omega_{BF}$0 for $\Omega_{BF}$1 and $\Omega_{BF}$2, and about 8–11 iterations to reduce the residual by $\Omega_{BF}$3 across fracture permeabilities from $\Omega_{BF}$4 to $\Omega_{BF}$5 and Forchheimer coefficients from $\Omega_{BF}$6 to $\Omega_{BF}$7.

A third numerical line couples generalized multiscale finite element methods with a multipoint flux mixed finite element method based on lowest-order $\Omega_{BF}$8 spaces and symmetric trapezoidal quadrature. The quadrature makes the velocity block locally eliminable, yielding a cell-centered symmetric positive definite pressure system,

$\Omega_{BF}$9

Offline snapshot spaces, spectral reduction, and residual-driven online enrichment are then used to resolve high-contrast heterogeneous Darcy–Forchheimer media. In the reported examples, Newton iteration counts were 7, 9, 11, 12, and 14 for $\Omega_D$0, whereas Picard required 48, 153, 491, 1376, and 3057 iterations in the same sequence; a second example showed Newton 7–14 versus Picard 56–4257 (Caucao et al., 2023, Arrarás et al., 2018, He et al., 2020).

5. Reduced-order thermal and local-thermal-nonequilibrium extensions

In porous heat-sink optimization, the term two-layer Darcy–Forchheimer model denotes a thickness-averaged thermal-fluid layer of thickness $\Omega_D$1 coupled to a base-solid layer of thickness $\Omega_D$2. The fluid momentum equation is

$\Omega_D$3

with incompressibility

$\Omega_D$4

The temperature fields satisfy

$\Omega_D$5

$\Omega_D$6

The interfacial heat-transfer coefficient is derived from the two-layer theory as

$\Omega_D$7

The porous/void heterogeneity is represented by two design variables: $\Omega_D$8, which selects void versus porous phase, and $\Omega_D$9, which controls graded lattice density. RVE-calibrated maps provide $\Omega_1 \cup \gamma \cup \Omega_2$0, $\Omega_1 \cup \gamma \cup \Omega_2$1, and $\Omega_1 \cup \gamma \cup \Omega_2$2. Inlet pressure drops of $\Omega_1 \cup \gamma \cup \Omega_2$3, $\Omega_1 \cup \gamma \cup \Omega_2$4, and $\Omega_1 \cup \gamma \cup \Omega_2$5 Pa were used in the optimization cases, and full-scale validation reported approximately 20–30 percent higher maximum Nusselt numbers for optimized voided lattices than for conventional plate-fin and uniform lattice heat sinks while maintaining lower pressure losses.

A different reduced-order extension is the one-dimensional Darcy–Forchheimer model under local thermal nonequilibrium for transpiration cooling. There the momentum and continuity system is

$\Omega_1 \cup \gamma \cup \Omega_2$6

$\Omega_1 \cup \gamma \cup \Omega_2$7

The full LTNE energy model is

$\Omega_1 \cup \gamma \cup \Omega_2$8

$\Omega_1 \cup \gamma \cup \Omega_2$9

and the simplified system neglecting fluid conduction is

$\Omega_{BF}$00

The paper proves uniqueness of the temperature system and of the coupled mass–momentum problem under an interface-density condition, and numerical comparisons justify neglecting fluid conduction in the assembled-1D model. In the reported two-domain hot-gas/porous coupling, six iterations suffice (Saito et al., 30 Sep 2025, Müller et al., 2022).

6. Applications, interpretation, and modeling caveats

The principal application domain is flow through heterogeneous porous media with strong permeability contrast. In the matrix–fracture setting, Darcy flow in the matrix and Darcy–Forchheimer flow in the fracture is appropriate when matrix velocities remain moderate but fracture velocities are high enough that inertial corrections matter. The fracture coefficient $\Omega_{BF}$01 then increases effective resistance with velocity magnitude, producing larger pressure drops along the fracture and limiting throughput at high flow rates.

In the Brinkman–Forchheimer/Darcy setting, the Brinkman term $\Omega_{BF}$02 is not merely a numerical convenience. It regularizes velocity in high-permeability regions, supports $\Omega_{BF}$03-conforming velocities, and gives a natural traction expression through the symmetric gradient. When $\Omega_{BF}$04 is negligible, the momentum balance reduces formally to Darcy–Forchheimer,

$\Omega_{BF}$05

so the coupled framework can be interpreted as an extension of a pure two-layer Darcy–Forchheimer law. Conversely, in the Knabner–Roberts analysis, the hybrid Darcy/Forchheimer formulation is the weak limit of a fully Forchheimer system as matrix inertial coefficients tend to zero.

Several interface misconceptions recur in this area. First, pressure continuity is not universal: some formulations use traction continuity, others use a Robin-type pressure jump involving normal transmissibility $\Omega_{BF}$06, and manifold formulations may add tangential slip through a strongly positive operator $\Omega_{BF}$07. Second, two-layer need not imply two full-dimensional subdomains; a codimension-one fracture model can still constitute a two-layer Darcy–Forchheimer system in the mixed-dimensional sense. Third, nonlinearity is often localized. In the mixed-dimensional multigrid setting, only the fracture block is nonlinear, which is why uncoupled smoothers remain effective.

The assumptions behind the principal theories are also specific. The matrix–fracture weak-limit identification uses $\Omega_{BF}$08 and Sobolev embedding $\Omega_{BF}$09; the manifold transmission theory requires compact boundaryless manifolds of dimension $\Omega_{BF}$10 or $\Omega_{BF}$11, absence of non-trivial Killing fields, and sufficiently small data; the heat-sink reduction assumes steady incompressible laminar flow, fully developed velocity and temperature profiles across thickness, isotropic homogenized porous properties, and low-to-moderate Reynolds numbers. In the topology-optimization study, the reduced-order model overestimates velocities and Nusselt numbers by about 7–23 percent relative to full-scale finite-element validation, which is acceptable there as a low-fidelity surrogate during optimization.

Taken together, these results show that the two-layer Darcy–Forchheimer model is best understood as a framework rather than a single PDE system. Its core structure is the same across variants—layerwise constitutive laws, explicit transmission operators, and monotone nonlinear resistance—while the precise form of the layers, interfaces, and solution theory changes with the physics being represented (Knabner et al., 2013, Caucao et al., 2023, Kohr et al., 2016).