Multiscale Hyperbolic-Parabolic Models for Nonlinear Reactive Transport in Heterogeneously Fractured Porous Media

Abstract: We study nonlinear reactive transport in a layered porous medium separated by an $\varepsilon$-thin, highly heterogeneous fracture whose aperture and obstacle pattern vary periodically. Species transport in the bulk is governed by parabolic reaction--diffusion equations, coupled to a convection-diffusion-reaction problem in the fracture with nonlinear wall and obstacle reactions and Peclet number of order $O(\varepsilon<sup>{-1})$. Via multiscale analysis as $\varepsilon \to 0$, when the fracture collapses to a flat interface, we derive a new type of homogenized model consisting of bulk diffusion--reaction equations coupled through nonlinear interface conditions and a first-order semilinear hyperbolic system on the interface. We prove well-posedness and regularity of the limit system, construct a multiscale approximation with boundary-layer correctors, and derive quantitative error estimates in suitable energy norms.

• The paper derives a homogenized model coupling nonlinear parabolic bulk transport with a semilinear hyperbolic fracture interface, including averaged wall and obstacle reactions in the critical scaling regime.
• The analysis proves well-posedness and classical regularity of the effective system, while boundary-layer constructions capture oscillatory walls, perforations, and outflow effects omitted by standard homogenization.
• Quantitative estimates show bulk errors of order O(ε^{1/2}) and fracture errors of order O(ε), providing a rigorous basis for reduced models of convection-dominated reactive flow in heterogeneous fractures.

Setting and microscopic model

The paper studies nonlinear reactive transport of $\mathcal{M}$ chemical species in a two-dimensional porous medium consisting of two bulk domains separated by a thin fracture whose aperture is of order $\mathcal{O}(\varepsilon)$ and whose walls oscillate $\varepsilon$-periodically; the fracture additionally contains a periodic pattern of internal obstacles (2602.16439). The bulk transport is governed by parabolic reaction-diffusion systems with constant diagonal diffusion tensors $\mathbb{D}^\pm$, while transport inside the perforated fracture is convection-dominated: diffusion is scaled by $\varepsilon$, so that the Péclet number is of order $\varepsilon^{-1}$. The velocity field in the fracture is conservative, generated by a potential solving a cell problem on the periodicity cell $Y_0$ with prescribed inflow profile $\mathrm{v}_0$; its longitudinal average $\hat{\mathrm{v}} = \langle v_1\rangle_{Y_0}$ is strictly positive, which fixes the direction of convective transport along the interface.

The coupling between subdomains occurs through nonlinear boundary conditions on three distinct parts of the fracture boundary: side-dependent flux conditions $\boldsymbol{\Upsilon}^\pm$ on the oscillating upper/lower walls, exchange fluxes $\mathcal{O}(\varepsilon)$0 scaled by $\mathcal{O}(\varepsilon)$1 on those same walls as seen from the fracture, and reactions $\mathcal{O}(\varepsilon)$2 scaled by $\mathcal{O}(\varepsilon)$3 on the obstacle boundaries $\mathcal{O}(\varepsilon)$4. The paper focuses on the critical regime $\mathcal{O}(\varepsilon)$5, in which the interfacial reactions act at the same asymptotic order as bulk and fracture transport and therefore enter the leading terms of the homogenized model; other scalings are explicitly deferred to future work.

The data assumptions are strong but standard for quantitative asymptotics: all nonlinearities are $\mathcal{O}(\varepsilon)$6-smooth, uniformly bounded with derivatives, Lipschitz continuous, vanish at zero, have compact support in space and time, and the interfacial chemistry is inactive on an initial interval $\mathcal{O}(\varepsilon)$7 — a delay used crucially later to exclude corner singularities when proving regularity.

Formal asymptotics and structure of the limit model

In the fracture, the ansatz combines a macroscopic profile $\mathcal{O}(\varepsilon)$8 with periodic correctors $\mathcal{O}(\varepsilon)$9 and $\varepsilon$0 solving convection-diffusion cell problems on $\varepsilon$1. Solvability of these problems follows from the Fredholm alternative applied to a model cell problem whose homogeneous adjoint admits only constants as solutions — a consequence of the incompressibility and impermeability of the velocity field. The resulting solvability condition produces the averaged reaction terms

$\varepsilon$2

which drive a first-order semilinear hyperbolic system on the collapsed interface $\varepsilon$3:

$\varepsilon$4

This hyperbolic limit subsystem is the qualitative novelty relative to prior thin-layer homogenizations such as Gahn–Neuss-Radu–Knabner or Pop–Bogers–Kumar, which obtain parabolic or purely reactive effective models; among related works only Mikelić–van Duijn and Mel'nyk–Rohde produce hyperbolic dispersion limits, none in this fractured-porous-medium geometry with nonlinear side- and obstacle-dependent reactions.

In the bulk domains, outer expansions are matched with inner boundary-layer expansions near the oscillating walls via harmonic correctors $\varepsilon$5 and $\varepsilon$6 posed on semi-infinite periodic half-strips $\varepsilon$7. Because $\varepsilon$8 are assumed even, $\varepsilon$9 is odd about $\mathbb{D}^\pm$0, so it contributes no first-order correction to the effective flux, and $\mathbb{D}^\pm$1 satisfies a normalization $\mathbb{D}^\pm$2 ensuring linear growth at infinity. Matching yields the homogenized transmission law

$\mathbb{D}^\pm$3

with $\mathbb{D}^\pm$4 the tangential average over the wall profile. A remark notes explicitly that if $\mathbb{D}^\pm$5 were not even, an additional term involving $\mathbb{D}^\pm$6 multiplying $\mathbb{D}^\pm$7 appears and the justification requires further corner layers — a clear limitation of the current framework tied to geometric symmetry.

Well-posedness and regularity of the homogenized system

Existence and uniqueness of a weak solution to the coupled parabolic–hyperbolic system are proved by an iterative scheme alternating linear parabolic solves (via Fourier sine–cosine representations) with explicit characteristic integration of the hyperbolic subsystem. Convergence relies on two abstract lemmas on recursive inequalities — one involving a contraction factor $\mathbb{D}^\pm$8, one with summable constants — applied to uniform norms of trace differences and to energy norms. The contraction is achieved on small time intervals ($\mathbb{D}^\pm$9), and the solution is extended to $\varepsilon$0 by successive restarts. Uniqueness follows by the same estimates on differences of two solutions.

A substantial regularity result establishes that the weak solution is classical: $\varepsilon$1 up to the closure of the space-time cylinders and $\varepsilon$2. The proof proceeds through a bootstrapping argument combining parabolic Schauder theory with Hölder estimates for the semilinear hyperbolic problem obtained directly along characteristics. Mixed Dirichlet–Neumann corners at $\varepsilon$3 and $\varepsilon$4 are handled by odd reflection across the vertical edges, using the compact-support and delayed-onset assumptions on the data to preserve compatibility. This regularity is not merely technical: it supplies exactly the smoothness required for the residual estimates underlying the error analysis.

Boundary layers near the outflow boundary

Because the regular part of the expansion does not satisfy the Dirichlet data on the vertical segment $\varepsilon$5, the paper constructs boundary-layer correctors $\varepsilon$6, $\varepsilon$7 on a periodically perforated semi-infinite strip $\varepsilon$8. The key analytical ingredient is Lemma 4.5: the stationary convection-diffusion problem on $\varepsilon$9 has a unique $\varepsilon^{-1}$0 solution decaying exponentially together with its first and second derivatives as $\varepsilon^{-1}$1. The proof combines truncated-domain approximations, energy estimates exploiting the divergence structure of $\varepsilon^{-1}$2 (the convective term integrates to a vanishing boundary flux), the maximum principle with the Hopf–Oleinik lemma, and a geometric decay argument showing that slice maxima contract by a factor $\varepsilon^{-1}$3 per period, giving decay rate $\varepsilon^{-1}$4. Derivative decay then follows from scale-invariant interior and boundary Schauder estimates on period-sized covers. This exponential stabilization justifies cutting off the boundary layer away from the outflow edge.

Approximation and error estimates

The global approximation $\varepsilon^{-1}$5 assembles the homogenized solution, bulk boundary-layer correctors (cut off by $\varepsilon^{-1}$6 away from the walls), the fracture corrector $\varepsilon^{-1}$7, and the outflow boundary layer $\varepsilon^{-1}$8 cut off at scale $\varepsilon^{-1}$9 with $Y_0$0 — an exponent chosen to balance the residuals left by truncation against the size of the layer. Small $Y_0$1 mismatches of the Dirichlet traces at the fracture ends are neutralized by additional cutoff functions contributing only $Y_0$2 in $Y_0$3.

The main theorem gives quantitative justification:

$Y_0$4

$Y_0$5

The proof couples three energy estimates across the interface: the bulk estimate contains a surface term measuring the fracture error, the fracture estimate contains weighted surface integrals of the bulk errors, and a trace inequality specific to thin perforated domains closes the loop. Gronwall arguments in time then yield the stated rates. As corollaries, the raw convergence rates follow: $Y_0$6 for bulk concentrations against $Y_0$7, and $Y_0$8 for fracture concentrations against $Y_0$9; in the cylindrical (non-perforated) special case, the cross-sectionally averaged fracture concentration converges at rate $\mathrm{v}_0$0. Notably, the approximation reveals rapidly oscillatory behavior of both the solution and its gradient along the interfaces in the bulk domains — information invisible at the homogenized level alone, illustrating why corrector-level results carry more physical content than convergence statements.

Limitations and open questions

Several restrictions bound the scope of the results. The regime $\mathrm{v}_0$1 is treated exclusively; formal calculations indicate that for $\mathrm{v}_0$2 the homogenized problem decouples into independent bulk problems with a trivial fracture state at leading order, but detecting fracture effects at higher asymptotic orders remains open. Evenness of the aperture functions $\mathrm{v}_0$3 is essential to the present justification of the bulk asymptotics near vertical boundaries. The rectangle geometry is assumed, though the authors note extension to smooth planar domains meeting a compatibility condition should be routine. All results are two-dimensional; three-dimensional configurations with connected perforated plates are proposed but untreated. Finally, no discretization of the mixed hyperbolic–parabolic homogenized system exists within standard finite element or finite volume frameworks, and the paper identifies the design of dedicated numerical schemes — enabled by the proven regularity — as an open task.

Conclusion

The paper delivers a complete rigorous program for reactive transport through a heterogeneously fractured porous medium under convection-dominated fracture dynamics: derivation of a new homogenized model coupling nonlinear diffusion-reaction systems in the bulks to a semilinear hyperbolic interface system through averaged nonlinear transmission laws, proof of well-posedness and $\mathrm{v}_0$4 regularity of the limit problem, construction of a multiscale approximation with wall, obstacle, and outflow boundary-layer correctors, and energy-norm error estimates of order $\mathrm{v}_0$5 in the fracture and $\mathrm{v}_0$6 in the bulk. Together these results provide one of the few quantitatively justified dimensionally-reduced models in which microscale geometry, nonlinear wall chemistry, and strong advection jointly shape the macroscopic transport.