Higher-Order Multiscale Computational Method for Multi-Continuum Problems in Highly Heterogeneous Media

Abstract: This paper presents a high-accuracy higher-order multiscale method for solving multi-continuum problems in in highly heterogeneous media. First, microscopic unit cell functions are defined, leading to the derivation of macroscopic homogenized equations and formulas for calculating effective parameters, which yield a higher-order multi-scale (HOMS) asymptotic solution. Subsequently, the pointwise approximation properties of this solution to the original equations are analyzed, and its convergence rate in the integral norm is rigorously established under certain assumptions. Furthermore, a multiscale numerical algorithm is developed by integrating the finite element method (FEM), finite difference method, and interpolation technique. Finally, numerical experiments demonstrate the high accuracy, efficiency, and stability of the proposed HOMS numerical algorithm.

• The paper develops a second-order multiscale expansion for two-continuum parabolic flow with strong exchange, introducing cell problems and cross-gradient coupling terms absent from first-order models.
• The analysis shows that first-order approximations retain an O(1) residual, while the higher-order method achieves an O(ε) error bound in uniform-in-time L² and H¹ norms under symmetry and H⁴-regularity assumptions.
• The numerical method reproduces fine-scale pressure oscillations with roughly 1–7% L² errors and 2.7–3.5× speedups over fine-mesh FEM, although gradient errors reach about 17% in a challenging 3D channel case.

Problem setting and motivation

The paper addresses transient flow in media where two overlapping continua—typically fractures and porous matrix—coexist at the same spatial points and exchange mass through a strong inter-continuum coupling term. The governing system is a pair of parabolic equations for the pressures $u_1^\varepsilon$ and $u_2^\varepsilon$, with multiscale porosities $$c_l^\varepsilon(\boldsymbol{x}) = c_l(\boldsymbol{x}/\varepsilon)$$, permeabilities $$\kappa_{l,ij}^\varepsilon(\boldsymbol{x}) = \kappa_{l,ij}(\boldsymbol{x}/\varepsilon)$$, and an exchange term scaled as $$\frac{1}{\varepsilon}Q_l^\varepsilon(\boldsymbol{x})(u_{3-l}^\varepsilon - u_l^\varepsilon)$$. This $O(1/\varepsilon)$ scaling of the transfer coefficient is the defining feature of the multi-continuum regime: it enforces near-equilibrium between continua at leading order while retaining a nontrivial $O(1)$ effective exchange in the homogenized limit. The coefficients are assumed 1-periodic in the microscale variable, uniformly elliptic and bounded (assumptions $(A_1)$–$(A_3)$).

The authors position this work against two lines of literature: single-continuum homogenization of Stokes/Darcy systems (Allaire; Hillairet), and multi-continuum modeling from Barenblatt's dual-porosity theory through MINC, triple-continuum models, and recent multicontinuum homogenization frameworks (Park and Hoang; Xie et al.; Ammosov et al.). Their stated gap is that prior multiscale treatments of multi-continuum problems deliver only first-order accuracy, which they show to be insufficient when $\varepsilon$ is a fixed engineering parameter.

The higher-order multiscale (HOMS) construction

The solution is expanded as $u_2^\varepsilon$0, and matched orders of $u_2^\varepsilon$1 after applying the chain rule. The leading term is independent of the fast variable. The first corrector takes the form

$u_2^\varepsilon$2

where $u_2^\varepsilon$3 solves the standard gradient cell problem and $u_2^\varepsilon$4 solves a new cell problem driven by the exchange coefficient $u_2^\varepsilon$5. A notable structural consequence appears in the homogenized system: beyond the classical effective storage $u_2^\varepsilon$6, permeability $u_2^\varepsilon$7, and exchange coefficient $u_2^\varepsilon$8, cross-coupling terms $u_2^\varepsilon$9 and $$c_l^\varepsilon(\boldsymbol{x}) = c_l(\boldsymbol{x}/\varepsilon)$$0 appear, so each continuum's homogenized equation contains gradient terms of the other continuum's pressure. These arise because the first corrector couples pressure differences across continua through $$c_l^\varepsilon(\boldsymbol{x}) = c_l(\boldsymbol{x}/\varepsilon)$$1, and their formulas involve both $$c_l^\varepsilon(\boldsymbol{x}) = c_l(\boldsymbol{x}/\varepsilon)$$2 and $$c_l^\varepsilon(\boldsymbol{x}) = c_l(\boldsymbol{x}/\varepsilon)$$3 cell functions.

The second-order corrector involves five families of auxiliary cell functions per continuum—$$c_l^\varepsilon(\boldsymbol{x}) = c_l(\boldsymbol{x}/\varepsilon)$$4, $$c_l^\varepsilon(\boldsymbol{x}) = c_l(\boldsymbol{x}/\varepsilon)$$5, $$c_l^\varepsilon(\boldsymbol{x}) = c_l(\boldsymbol{x}/\varepsilon)$$6, $$c_l^\varepsilon(\boldsymbol{x}) = c_l(\boldsymbol{x}/\varepsilon)$$7, and $$c_l^\varepsilon(\boldsymbol{x}) = c_l(\boldsymbol{x}/\varepsilon)$$8—each defined by a Dirichlet cell problem on $$c_l^\varepsilon(\boldsymbol{x}) = c_l(\boldsymbol{x}/\varepsilon)$$9. All cell problems use homogeneous Dirichlet rather than periodic boundary conditions; the authors justify this by citing prior work showing Dirichlet conditions are admissible replacements and are computationally more convenient. This is a deliberate deviation from classical periodic-cell theory whose error consequences are absorbed into the estimates below.

Error analysis

Two results structure the approximation theory. First, pointwise analysis of the residual problems shows that the FOMS solution leaves an $$\kappa_{l,ij}^\varepsilon(\boldsymbol{x}) = \kappa_{l,ij}(\boldsymbol{x}/\varepsilon)$$0 residual (the source $$\kappa_{l,ij}^\varepsilon(\boldsymbol{x}) = \kappa_{l,ij}(\boldsymbol{x}/\varepsilon)$$1 does not vanish as $$\kappa_{l,ij}^\varepsilon(\boldsymbol{x}) = \kappa_{l,ij}(\boldsymbol{x}/\varepsilon)$$2), whereas the HOMS solution reduces the residual to $$\kappa_{l,ij}^\varepsilon(\boldsymbol{x}) = \kappa_{l,ij}(\boldsymbol{x}/\varepsilon)$$3. The authors draw the practical conclusion directly: since $$\kappa_{l,ij}^\varepsilon(\boldsymbol{x}) = \kappa_{l,ij}(\boldsymbol{x}/\varepsilon)$$4 is fixed in applications, only the second-order expansion delivers acceptable accuracy—a claim that motivates the entire construction and is confirmed numerically.

Second, under additional symmetry assumptions $$\kappa_{l,ij}^\varepsilon(\boldsymbol{x}) = \kappa_{l,ij}(\boldsymbol{x}/\varepsilon)$$5–$$\kappa_{l,ij}^\varepsilon(\boldsymbol{x}) = \kappa_{l,ij}(\boldsymbol{x}/\varepsilon)$$6 (coefficients symmetric about mid-planes of the unit cell, ensuring continuity of conormal derivatives of all cell functions on $$\kappa_{l,ij}^\varepsilon(\boldsymbol{x}) = \kappa_{l,ij}(\boldsymbol{x}/\varepsilon)$$7), and regularity $$\kappa_{l,ij}^\varepsilon(\boldsymbol{x}) = \kappa_{l,ij}(\boldsymbol{x}/\varepsilon)$$8, the main theorem establishes

$$\kappa_{l,ij}^\varepsilon(\boldsymbol{x}) = \kappa_{l,ij}(\boldsymbol{x}/\varepsilon)$$9

The proof proceeds via the energy method: the boundary flux terms over interior cell interfaces cancel exactly thanks to Lemma on conormal continuity, Poincaré–Friedrichs and Young inequalities bound the energy terms, and Gronwall's inequality closes the estimate. The convergence rate is first order in $$\frac{1}{\varepsilon}Q_l^\varepsilon(\boldsymbol{x})(u_{3-l}^\varepsilon - u_l^\varepsilon)$$0 in both $$\frac{1}{\varepsilon}Q_l^\varepsilon(\boldsymbol{x})(u_{3-l}^\varepsilon - u_l^\varepsilon)$$1 and $$\frac{1}{\varepsilon}Q_l^\varepsilon(\boldsymbol{x})(u_{3-l}^\varepsilon - u_l^\varepsilon)$$2 norms uniformly in time—an improvement over the pointwise $$\frac{1}{\varepsilon}Q_l^\varepsilon(\boldsymbol{x})(u_{3-l}^\varepsilon - u_l^\varepsilon)$$3 behavior of the FOMS ansatz. It should be noted that the rate depends on the plane-symmetry assumption $$\frac{1}{\varepsilon}Q_l^\varepsilon(\boldsymbol{x})(u_{3-l}^\varepsilon - u_l^\varepsilon)$$4, which restricts the class of unit cells for which the theorem holds.

Numerical algorithm and experiments

The algorithm combines FEM on the unit cell (for all auxiliary functions), a Crank–Nicolson-type coupled finite element/finite difference scheme in time for the homogenized system, element-averaged gradients for reconstructing $$\frac{1}{\varepsilon}Q_l^\varepsilon(\boldsymbol{x})(u_{3-l}^\varepsilon - u_l^\varepsilon)$$5 and second derivatives, and interpolation to evaluate the assembled HOMS solution at arbitrary points. Four examples validate the method against fine-mesh FEM reference solutions:

Example	Domain	$$\frac{1}{\varepsilon}Q_l^\varepsilon(\boldsymbol{x})(u_{3-l}^\varepsilon - u_l^\varepsilon)$$6	FEM time	HOMS time	Rel. err. $$\frac{1}{\varepsilon}Q_l^\varepsilon(\boldsymbol{x})(u_{3-l}^\varepsilon - u_l^\varepsilon)$$7 ($$\frac{1}{\varepsilon}Q_l^\varepsilon(\boldsymbol{x})(u_{3-l}^\varepsilon - u_l^\varepsilon)$$8)	Rel. err. $$\frac{1}{\varepsilon}Q_l^\varepsilon(\boldsymbol{x})(u_{3-l}^\varepsilon - u_l^\varepsilon)$$9 ($O(1/\varepsilon)$0)
2D porous	$O(1/\varepsilon)$1	1/8	168 s	61 s	~6%, 6%	~8%, 9%
3D porous	$O(1/\varepsilon)$2	1/5	8970 s	2629 s	~6%, 1%	~12%, 10%
2D channel	$O(1/\varepsilon)$3	1/8	232 s	77 s	~7%, 5%	~10%, 8%
3D channel	$O(1/\varepsilon)$4	1/4	14893 s	4218 s	~2%, 2%	~17%, 15%

The speedups range from roughly 2.7× (2D) to 3.5× (3D), achieved with far coarser meshes—for instance, 565,341 elements for the multiscale reference versus 93,750 + 83,536 for the homogenized and cell problems in Example 2. Across all cases, the HOMS solution substantially outperforms both the homogenized and FOMS solutions, and only the HOMS reconstruction reproduces local oscillations of the fine-scale pressure field. Error curves remain stable over $O(1/\varepsilon)$5, supporting applicability to time-dependent problems. The largest observed errors occur in the $O(1/\varepsilon)$6 seminorm for the 3D channel case (~17% and 15%), indicating that gradient accuracy degrades most in high-contrast channel geometries even though $O(1/\varepsilon)$7 accuracy there is excellent.

Limitations and open questions

The paper's scope carries several explicit restrictions. The model covers exactly two interacting continua with linear exchange; extension to $O(1/\varepsilon)$8 continua, nonlinear (non-Darcy, poromechanical, multiphase) coupling, random (non-periodic) microstructures, and three-scale (micro–meso–macro) systems are all identified by the authors as unresolved directions. The convergence theorem requires the mid-plane symmetry assumption $O(1/\varepsilon)$9 and $O(1)$0 regularity of the homogenized solution, neither of which is verified for general geometries or rough data. All validation uses FEM reference solutions rather than analytical solutions, since closed-form solutions of the multiscale problem are unavailable; consequently the reported error percentages measure agreement with a very fine discretization, not exact error. Finally, the reported speedups are moderate (under 4×); whether the offline cell-problem cost amortizes better in large-scale or many-query settings is not examined.

Conclusion

The paper extends second-order two-scale asymptotic analysis to dual-continuum parabolic systems with $O(1)$1 inter-continuum exchange, deriving a homogenized system with cross-gradient coupling terms, proving an $O(1)$2 uniform-in-time error estimate in $O(1)$3 and $O(1)$4 norms under symmetry assumptions, and implementing a hybrid FEM/finite difference algorithm. Numerical results confirm that the second-order correction is necessary—first-order expansions leave $O(1)$5 pointwise error—and demonstrate relative errors of 1–7% in $O(1)$6 with computational savings of roughly 3× relative to fine-grid FEM. The framework is limited to linear, periodic, two-continuum settings, and its extension to nonlinear, stochastic, and multi-scale-hierarchical regimes remains open.