A Hybrid Model Reduction Method for Dual-Continuum Model with Random Inputs

1. Introduction

Numerous engineering disciplines—such as groundwater hydrology, hydrocarbon reservoir engineering, and composite materials science—study systems composed of multiple interacting continuous media. These systems exhibit complex behaviors arising from multi-component, multi-scale, and multi-physics coupling, which traditional single-continuum models, relying on homogeneity assumptions, often fail to capture. Standard dual-continuum models address this limitation by coupling mass transfers between two separate continua of distinct scales. However, due to incomplete knowledge about the physical properties and the presence of measurement noise, such models inherently contain uncertainties. In this work, we consider the numerical simulation methods for dual-continuum problems and discuss their related applications.

For the numerical simulation of dual-continuum model, there exist some efficient methods, including finite element method [11,15], finite volume method [9,18] and hybrid schemes that combine both finite element method and finite volume method [20,24]. In [21], the authors employ a finite element method with standard linear basis functions for spatial approximation, and combine implicit and explicit temporal schemes regarding time discretization. However, due to the complex heterogeneity of porous media—particularly the presence of multiple scales and high contrast—some form of model reduction is often necessary for practical flow simulations. The application of homogenization theory provides one pathway to deriving dual-continuum models, as illustrated by [1], who extended the earlier framework of [4]. Their work established a general double-porosity model for single-phase flow in fractured reservoirs, explicitly describing transport within and between each continuum. Further extensions have been proposed to enhance the representation of inter-continuum interactions: Wu and Pruess [33] employed multiple interacting continua to better capture transient inter-continuum flow, and Yan et al. [34] later generalized this concept into a versatile multi-continuum framework capable of handling an arbitrary number of continua.

While such models capture flow both inside and across continua, their accuracy relies on a strong separation of scales. In cases of weak scale separation, homogenized effective parameters can lose physical meaning, significantly compromising model fidelity. To overcome the limitations of classical homogenization and directly incorporate multiscale heterogeneity in a computationally efficient manner, the Multiscale Finite Element Method (MsFEM) [12,13] and its generalization, the Generalized Multiscale Finite Element Method (GMsFEM) [14], have been developed. These approaches integrate fine-scale heterogeneity while reducing computational cost. For simulating gas transport in the shale formation, the work of I. Akkutlu et al. in [2] develop a coupled multiscale and multi-continuum approach. Furthermore, active research efforts continue to develop novel model reduction techniques, such as the constraint energy minimizing generalized multiscale finite element method (CEM-GMsFEM) [6,7,17]. Related numerical methods for multi-continuum systems [5], including the non-local multi-continuum method (NLMC) [8,32], are also being advanced. These approaches are particularly effective in addressing the challenges posed by high-contrast properties and multiscale features inherent in multi-continuum media.

When employing a multiscale approach to address heterogeneity in model parameters, solving models with uncertain parameters remains a significant challenge. Typically, such uncertain parameters are characterized by a finite set of random variables. This representation allows the modeling framework to be described by Parameterized Partial Differential Equations (PPDEs) [25,29]. The work [23] presents some Reduced Basis (RB) methods for fluid infiltration problems through certain porous media modeled as dual-continuum with localized uncertainties. In this work, we focus on developing the reduced basis approximations by combining the uncoupled GMsFEM method and the Variable-separation (VS) method for dual-continuum models with localized uncertainties. The VS method aims to approximate a solution in the form:

$$p(x,t;\mu )\approx \sum _{i=1}^N{\zeta }_i\left(\mu \right)P_i(t,x),$$

where N is the number of separated terms, $\left\{P_i(t,x)\right\}$ are deterministic basis functions of space x variable and time variable t, and $\left\{{\zeta }_i\left(\mu \right)\right\}$ are stochastic basis functions of the random parameter $\mu$. This separated representation can also be achieved by other methods, such as the Proper Generalized Decomposition (PGD) [26,27,28] and the Empirical Interpolation Method (EIM) [3,19,30]. However, the VS method offers distinct advantages. Compared to PGD, which typically requires numerous iterations with an arbitrary initial guess at each enrichment step, the VS method is more computationally efficient. While the VS method shares similarities with EIM, it does not require the selection of optimal parameter values and interpolation nodes. Furthermore, for online computation, the VS method directly evaluates the solution using the separated representation. In contrast, EIM requires solving an algebraic system based on pre-computed interpolation nodes and parameter values. Consequently, the VS method is not only easier to implement but also enables faster online computation.

The organization of the paper is as follows. In Section 2, we introduce some preliminaries and notations. Section 3 is about problem formulation, where we show the existence and uniqueness of the weak solution, and provide fine-scale finite element discretization. In Section 4, an overview of the VS method is given. In Section 5, we describe the uncoupled GMsFEM and hybrid model reduction method. Numerical examples are presented in Section 6 to demonstrate the effectiveness of the proposed hybrid model reduction method. Finally, conclusions and closing remarks are provided.

2. Mathematical Model

In this section, we present some preliminaries and notations for the rest of paper. We first define the close set $[t_0,t_f]$ as the time interval. Let $(\Omega ,\mathcal{F},\mathcal{P})$ be a complete probability space, where $\Omega$ is the space of elementary events, $\mathcal{F}\subset 2^{\Omega }$ is the $\sigma$-algebra, and $\mathcal{P}$ is the probability measure on $\Omega$. For computation, the random field needs to be approximated by a prescribed number of random variables, $\mu \left(\omega \right)=\{{\mu }_1\left(\omega \right),{\mu }_2\left(\omega \right),\dots ,{\mu }_n\left(\omega \right)\}$, i.e., $\mu (·):\Omega \to \Gamma \subset {\mathbb{R}}^n(n\ge 1)$. We use the random vector $\mu \left(\omega \right)$ to characterize the stochastic property of the random field.

We denote the space of the random variables with second order moments by $L_{\mathcal{P}}^2(\Omega )$, which is defined by

$$L_{\mathcal{P}}^2(\Omega )=\{\mu :\omega \in \Omega \to \xi \left(\omega \right)\in \mathbb{R};{\int }_{\Omega }\mu {\left(\omega \right)}^2\mathcal{P}\left(d\mu \right)<\infty \}.$$

The inner product in $L_{\mathcal{P}}^2(\Omega )$ is defined by ${(\mu ,\xi )}_{L_{\mathcal{P}}^2(\Omega )}:={\int }_{\Omega }\mu \left(\omega \right)\xi \left(\omega \right)\mathcal{P}\left(d\mu \right)$, which induces the norm ${\parallel \mu \parallel }_{L^2}^2={\parallel \mu \parallel }_{L_{\mathcal{P}}^2}^2:={(\mu ,\mu )}_{L_{\mathcal{P}}^2}$.

Let $L^2\left(D\right)$ and $H^s\left(D\right)$ be the usual Lebesgue space and Sobolev space [16], respectively. With the notation defined above, we define the following tensor-product Hilbert space

$$L^2(t_s,t_f;L_{\mathcal{P}}^{\rho }(\Gamma ;H^s\left(D\right))):=\{\nu :[t_s,t_f]\times D\times \Gamma \to {\mathbb{R}|\parallel \nu \parallel }_{L^2(t_s,t_f;L_{\mathcal{P}}^{\rho }(\Gamma ;H^s(\Omega )))}<\infty \},$$

and the equipped norm is

$$\begin{array}{cc}\hfill {\parallel \nu \parallel }_{L^2(t_s,t_f;L_{\mathcal{P}}^{\rho }(\Gamma ;H^s\left(D\right)))}^2& ={\int }_{t_s}^{t_f}{\parallel \nu \parallel }_{L_{\mathcal{P}}^{\rho }(\Gamma ;H^s\left(D\right))}^{2}dt\hfill \\ \hfill & ={\int }_{t_s}^{t_f}({\int }_{\Gamma }{\parallel \nu (·,·;\gamma )\parallel }_{H^s\left(D\right)}^2\rho d_{\gamma })dt.\hfill \end{array}$$

To shorten the notation, we denote

$${\mathcal{H}}^s\left(D_t\right):=L^2(t_s,t_f;L_{\mathcal{P}}^{\rho }(\Gamma ;H^s\left(D\right))).$$

In particular,

$${\mathcal{H}}_0^1\left(D_t\right)=\{\nu \in {\mathcal{H}}^1\left(D_t\right):\nu {|}_{\partial D}=0\}.$$

When $s=0$, we employ the abbreviated notion ${\mathcal{L}}^2\left(D_t\right)$ to denote ${\mathcal{H}}^0\left(D_t\right)$ by convention.

For deterministic situation, the tensor spaces can be simplified by

$${\mathcal{L}}^2\left({\Omega }_t\right):=L^2(t_s,t_f;L^2\left(D\right)),{\mathcal{H}}^1\left({\Omega }_t\right):=L^2(t_s,t_f;H^1\left(D\right)).$$

3. Dual Continuum Model

In this paper, We consider the dual continuum model with localized uncertainties:

$$\begin{array}{c}\hfill \left\{\begin{array}{c}\begin{array}{cc}& {\displaystyle \frac{\partial p_1(x,t;\mu )}{\partial t}}-\mathrm{div}({\kappa }_1\nabla p_1(x,t;\mu ))+\sigma \left(x\right)(p_1(x,t;\mu )-p_2(x,t;\mu ))=f_1(x,t;\mu ),\hfill \\ & {\displaystyle \frac{\partial p_2(x,t;\mu )}{\partial t}}-\mathrm{div}({\kappa }_2\nabla p_2(x,t;\mu ))+\sigma \left(x\right)(p_2(x,t;\mu )-p_1(x,t;\mu ))=f_2(x,t;\mu ),\hfill \end{array}\hfill \end{array}\right.\end{array}$$

(1)

in a computational domain $D\subset {\mathbb{R}}^2$. For $i=1,2$, let $p_i$ denote the pressure, ${\kappa }_i$ the permeability, and $f_i$ the source term for the i-th continuum. The continua are coupled through mass exchange, with $\sigma$ representing the strength of the mass transfer between them. A notable application of the dual continuum model (1) is its ability to capture the global interactive effects between unresolved fractures and the surrounding matrix. In this work, we consider high-contrast, channelized media. The initial condition is prescribed as $p_i(x,t_0;\mu )=p_i^0(x;\mu )$ in D, and the boundary condition is given by $p_i(x,t;\mu )=g_i(t;\mu )$ on $\partial \mathcal{D}\times [t_0,t_f]$. Here, we assume the permeability fields ${\kappa }_i\left(x\right)(i=1,2)$ are uniformly bounded, i.e.,

$$\underline{\kappa }\le {\kappa }_i\left(x\right)\le \overline{\kappa },for\forall x\in \Omega .$$

and the transform term is also bounded, i.e,

$$\left|\sigma \left(x\right)\right|\le \overline{\sigma },for\forall x\in \Omega .$$

3.1. Variational Formulation of the Dual-Continuum Model

The model (1) leads to the variational problem as follows: for any test functions $v_i,i=1,2$,

$$\begin{array}{c}\hfill \left\{\begin{array}{c}\begin{array}{cc}& ({\displaystyle \frac{\partial p_1}{\partial t}},v_1)+a_1(p_1,v_1)+q(p_1-p_2,v_1)=s_1\left(v_1\right)\hfill \\ & ({\displaystyle \frac{\partial p_2}{\partial t}},v_2)+a_2(p_2,v_2)+q(p_2-p_1,v_2)=s_2\left(v_2\right),\hfill \end{array}\hfill \end{array}\right.\end{array}$$

(2)

where

$$\begin{array}{cc}\hfill & ({\displaystyle \frac{\partial p_1}{\partial t}},v_1)={\int }_D{\displaystyle \frac{\partial p_1}{\partial t}}{v}_{1}dx,({\displaystyle \frac{\partial p_2}{\partial t}},v_2)={\int }_D{\displaystyle \frac{\partial p_2}{\partial t}}{v}_{2}dx,\hfill \\ \hfill & a_1(p_1,v_1)={\int }_D{k}_1\nabla p_1\nabla v_{1}dx,a_2(p_2,v_2)={\int }_D{k}_2\nabla p_2\nabla v_{2}dx,\hfill \\ \hfill & q(p_m-p_f,v)={\int }_D\sigma (p_1-p_2)vdx,q(p_f-p_m,v)={\int }_D\sigma (p_2-p_1)vdx.\hfill \end{array}$$

and

$$s_1\left(v\right)={\int }_D{f}_1{v}_{1}dx-{\int }_D{\kappa }_1\nabla G_1\nabla v_{1}dx,s_2\left(v\right)={\int }_D{f}_1{v}_{2}dx-{\int }_D{\kappa }_2\nabla G_2\nabla v_{2}dx.$$

Here, $G_1$ and $G_2$ are continuous piecewise polynomials that match the values of $g_1$ and $g_2$, respectively, at the boundary nodes on $\partial D$.

Moreover, we assume that both $a_1(·,·;\mu )$ and $a_2(·,·;\mu )$ are affinely dependent on the parameter vector $\mu$, i.e., for some $Q_a,Q_b,Q_{u_0},Q_{\widehat{u}}\in \mathbb{N}$, the following affine decompositions hold:

$$\left\{\begin{array}{c}\begin{array}{cc}\hfill & a_1(w,v;\mu )=\sum _{i=1}^{Q_a}{Q}_1^i\left(\mu \right)a_1^i(w,v)\forall w,v\in H^1(\Omega ),\forall \mu \in \Gamma ,\hfill \\ \hfill & a_2(w,v;\mu )=\sum _{j=1}^{Q_b}{Q}_2^j\left(\mu \right)a_2^j(w,v)\forall w,v\in H^1(\Omega ),\forall \mu \in \Gamma ,\hfill \\ \hfill & s_1(v;\mu )=\sum _{n=1}^{N_a}{X}_n\left(\mu \right)s_{1,n}\left(v\right)\forall v\in H^1(\Omega ),\forall \mu \in \Gamma ,\hfill \\ \hfill & s_2(v;\mu )=\sum _{m=1}^{M_a}{Y}_m\left(\mu \right)s_{2,m}\left(v\right)\forall v\in H^1(\Omega ),\forall \mu \in \Gamma .\hfill \end{array}\hfill \end{array}\right.$$

(3)

If the bilinear forms $a_1(·,·;\mu )$ and $a_2(·,·;\mu )$ are not affine, an affine approximation can still be constructed using the VS method introduced in [10] or the Empirical Interpolation Method (EIM; see, e.g., [3,19,30]).

3.2. The Approximation Based on FEM

To obtain a reference (truth) solution, the finite element method (FEM) is used for the spatial discretization of the dual-continuum model, while the implicit Euler method is employed for time integration.

Let ${\mathcal{T}}_h$ be a uniform partition of the spatial domain D. On fine grid ${\mathcal{T}}_h$, the corresponding FE space is ${\mathbf{V}}^h\left(D\right)\subset H^1\left(D\right)\times H^1\left(D\right)$ and ${\mathbf{V}}_0^h(\Omega )\subset {\mathbf{V}}^h(\Omega )$ is the subspace with vanishing boundary values. Let $N_h$ be the number of vertices, and the dimension of FE space be $\mathcal{N}$. Typically, the dimension is very large. Furthermore, we suppose that $M^h(\Omega )$ is the finite dimensional subspace of $L^2(\Omega )$. Thus, we define

$${\mathcal{V}}^h\left({\Omega }_t\right):=L^2(t_s,t_f;{\mathbf{V}}^h(\Omega ))\subset {\mathcal{H}}^1\left({\Omega }_t\right),{\mathcal{M}}^h\left({\Omega }_t\right):=L^2(t_s,t_f;M^h(\Omega ))\subset {\mathcal{L}}^2\left({\Omega }_t\right).$$

With the notations defined above, the Galerkin formulation of the system (2) is to find $(p_{1,h},p_{2,h})\in {\mathcal{V}}^h\left({\Omega }_t\right)\times {\mathcal{V}}^h\left({\Omega }_t\right)$ such that

$$\begin{array}{c}\hfill \left\{\begin{array}{c}\begin{array}{cc}& ({\displaystyle \frac{\partial p_{1,h}}{\partial t}},v_{1,h})+a_1(p_{1,h},v_{1,h})+q(p_{1,h}-p_{2,h},v_{1,h})=s_{1,h}\left(v_{1,h}\right)\hfill \\ & ({\displaystyle \frac{\partial p_{2,h}}{\partial t}},v_{2,h})+a_2(p_{2,h},v_{2,h})+q(p_{2,h}-p_{1,h},v_{2,h})=s_{2,h}\left(v_{2,h}\right),\hfill \end{array}\hfill \end{array}\right.\end{array}$$

(4)

where both $v_{1,h}$ and $v_{2,h}$ are the respective test functions.

Applying the implicit Euler method to the system (4)yields the full discretized formulation of the dual-continuum model as follows:

$$\begin{array}{c}\hfill \left\{\begin{array}{c}\begin{array}{cc}& (p_{1,h}^{k+1},v_{1,h})+\Delta t{a}_1(p_{1,h}^{k+1},v_{1,h})+\Delta tq(p_{1,h}^{k+1}-p_{2,h}^{k+1},v_{1,h})=\Delta t{s}_{1,h}^{k+1}\left(v_{1,h}\right)+(p_{1,h}^k,v_{1,h}),\hfill \\ & (p_{2,h}^{k+1},v_{2,h})+\Delta t{a}_2(p_{2,h}^{k+1},v_{2,h})+\Delta tq(p_{2,h}^{k+1}-\Delta t{p}_{1,h}^{k+1},v_{2,h})=\Delta t{s}_{2,h}^{k+1}\left(v_{2,h}\right)+(p_{2,h}^k,v_{2,h}),\hfill \end{array}\hfill \end{array}\right.\end{array}$$

(5)

where $\Delta t:={\displaystyle \frac{t_f-t_s}{N_t}}$ is the time step, $t_k:=k\Delta t$, $0\le k\le N_t$. If the basis functions of space $V^h(\Omega )$ are denoted by ${\left\{{\varphi }_j\right\}}_{j=1}^{N_h}$, we can represent the variables in Eq. (5) with the linear combination of the corresponding basis functions as

$$p_{1,h}^k=\sum _{i=1}^{N_h}{p}_{1,h}^{k,i}\left(\mu \right){\varphi }_i,p_{2,h}^k=\sum _{i=1}^{N_h}{p}_{2,h}^{k,i}\left(\mu \right){\varphi }_i.$$

If affine assumption (3) holds, the full discretization (5) becomes

$$\begin{array}{c}\hfill \left\{\begin{array}{c}\begin{array}{cc}& \sum _{i=1}^{N_h}{P}_{1,h}^{k+1,i}({\varphi }_i,{\varphi }_l)+\Delta t\sum _{s=1}^{Q_a}\sum _{i=1}^{N_h}{Q}_1^s\left(\mu \right)P_{1,h}^{k+1,i}{a}_1^s({\varphi }_i,{\varphi }_l)+\Delta t\sum _{i=1}^{N_h}{P}_{1,h}^{k+1,i}q({\varphi }_i,{\varphi }_l)-\hfill \\ & \Delta t\sum _{j=1}^{N_h}{P}_{2,h}^{k+1,i}q({\varphi }_j,{\varphi }_l)=\Delta t\sum _{n=1}^{N_a}{X}_n\left(\mu \right)S_{1,n}^{k+1}\left({\varphi }_l\right)+\sum _{i=1}^{N_h}{P}_{1,h}^{k,i}({\varphi }_i,{\varphi }_l)\hfill \\ & \sum _{j=1}^{N_h}{P}_{2,h}^{k+1,j}({\varphi }_j,{\varphi }_m)+\Delta t\sum _{t=1}^{Q_b}\sum _{j=1}^{N_h}{Q}_2^t\left(\mu \right)P_{2,h}^{k+1,j}{a}_2^t({\varphi }_j,{\varphi }_m)+\Delta t\sum _{j=1}^{N_h}{P}_{2,h}^{k+1,j}q({\varphi }_j,{\varphi }_m)-\hfill \\ & \Delta t\sum _{i=1}^{N_h}{P}_{1,h}^{k+1,i}q({\varphi }_i,{\varphi }_m)=\Delta t\sum _{r=1}^{M_a}{Y}_r\left(\mu \right)S_{2,r}^{k+1}\left({\varphi }_m\right)+\sum _{j=1}^{N_h}{P}_{1,h}^{k,j}({\varphi }_j,{\varphi }_m),\hfill \end{array}\hfill \end{array}\right.\end{array}$$

(6)

where ${\varphi }_l$ and ${\varphi }_m$ are the test functions corresponding to the pressure variables.

Next, if we define

$$\begin{array}{c}\hfill \left\{\begin{array}{cc}\hfill & {\left(M_{1,h}\right)}_{i,l}=({\varphi }_i,{\varphi }_l),1\le i,l\le N_h,\hfill \\ \hfill & {\left(M_{2,h}\right)}_{i,l}=q({\varphi }_i,{\varphi }_l),1\le i,l\le N_h,\hfill \\ \hfill & {\left(K_{1,h}^s\right)}_{i,l}=a_1^s({\varphi }_i,{\varphi }_l),1\le i,l\le N_h\hfill \\ \hfill & {\left(K_{2,h}^t\right)}_{j,m}=a_2^t({\varphi }_j,{\varphi }_m),1\le j,m\le N_h\hfill \\ \hfill & {\left(B_{1,n}^{k+1}\right)}_l=S_{1,n}^{k+1}\left({\varphi }_l\right),1\le l\le N_h\hfill \\ \hfill & {\left(B_{2,r}^{k+1}\right)}_m=S_{2,r}^{k+1}\left({\varphi }_m\right),1\le m\le N_h\hfill \\ \hfill & K_{1,h}={\displaystyle \sum _{s=1}^{Q_a}}{Q}_1^s\left(\mu \right)K_{1,h}^s,K_{2,h}={\displaystyle \sum _{t=1}^{Q_b}}{Q}_2^t\left(\mu \right)K_{2,h}^t\hfill \\ \hfill & B_1^{k+1}={\displaystyle \sum _{n=1}^{N_a}}{X}_n\left(\mu \right)B_{1,n}^{k+1},B_2^{k+1}={\displaystyle \sum _{r=1}^{M_a}}{Y}_r\left(\mu \right)B_{2,r}^{k+1},\hfill \end{array}\right.\end{array}$$

(7)

then we can get the matrix form of full discrete system (6) as

$$\begin{array}{c}\hfill \left\{\begin{array}{c}\begin{array}{cc}& (M_{1,h}+\Delta t{K}_{1,h})P_{1,h}^{k+1}+\Delta t{M}_{2,h}(P_{1,h}^{k+1}-P_{2,h}^{k+1})=\Delta t{B}_1^{k+1}+M_{1,h}{P}_{1,h}^k\hfill \\ & (M_{1,h}+\Delta t{K}_{2,h})P_{2,h}^{k+1}+\Delta t{M}_{2,h}(P_{2,h}^{k+1}-P_{1,h}^{k+1})=\Delta t{B}_2^{k+1}+M_{2,h}{P}_{2,h}^k.\hfill \end{array}\hfill \end{array}\right.\end{array}$$

(8)

Combining these equations gives the matrix system

$$\left[\begin{array}{cc}{M}_{1,h}+\Delta t{K}_{1,h}+\Delta t{M}_{2,h}& -\Delta t{M}_{2,h}\\ -\Delta t{M}_{2,h}& M_{1,h}+\Delta t{K}_{2,h}+\Delta t{M}_{2,h}\end{array}\right]\left[\begin{array}{c}{P}_{1,h}^{k+1}\left(\mu \right)\\ P_{2,h}^{k+1}\left(\mu \right)\end{array}\right]=\left[\begin{array}{c}{F}_1^{k+1}\left(\mu \right)\\ F_2^{k+1}\left(\mu \right)\end{array}\right].$$

(9)

Here, $P_{1,h}^{k+1}\left(\mu \right)$ and $P_{2,h}^{k+1}\left(\mu \right)$ denote the discrete approximations of $p_1,p_2$ at the $(k+1)$-th time step, with the corresponding right-hand-side vectors given by

$$F_1^{k+1}\left(\mu \right):=\Delta t{B}_1^{k+1}+M_{1,h}{P}_{1,h}^k{F}_2^{k+1}\left(\mu \right):=\Delta t{B}_2^{k+1}+M_{1,h}{P}_{2,h}^k.$$

Under the affine assumption (3), the mass matrices $M_{1,h}$ and $M_{2,h}$, the stiffness matrices ${\left\{K_{1,h}^s\right\}}_{s=1}^{Q_a}$ and ${\left\{K_{2,h}^t\right\}}_{t=1}^{Q_b}$, and the right-hand-side vectors ${\left\{B_{1,n}^{k+1}\right\}}_{n=1}^{N_a}$ and ${\left\{B_{2,r}^{k+1}\right\}}_{r=1}^{M_a}$ can be precomputed once during the offline stage, despite the high computational cost involved. In the online stage, for each new parameter sample $\mu \in \Gamma$, all coefficients ${\left\{Q_1^s\left(\mu \right)\right\}}_{s=1}^{Q_a}$, ${\left\{Q_2^t\left(\mu \right)\right\}}_{t=1}^{Q_b}$, ${\left\{X_n\left(\mu \right)\right\}}_{n=1}^{N_a}$ and ${\left\{Y_r\left(\mu \right)\right\}}_{r=1}^{M_a}$ need to be updated.

4. A Variable-Separation Method for the Stochastic Dual-Continuum Model

In this section, we present a variable-separation method applicable to the stochastic dual-continuum model. The corresponding stochastic solutions are expressed in the following form.

$$\left\{\begin{array}{c}\begin{array}{cc}\hfill & p_1(x,t;\mu )\approx p_{1,N}(x,t;\mu ):=\sum _{i=1}^N{\xi }_i\left(\mu \right)P_{1,i}(t,x)\hfill \\ \hfill & p_2(x,t;\mu )\approx p_{2,N}(x,t;\mu ):=\sum _{i=1}^N{\lambda }_i\left(\mu \right)P_{2,i}(t,x).\hfill \end{array}\hfill \end{array}\right.$$

(10)

Here, N denotes the number of separation terms. The stochastic basis functions ${\left\{{\xi }_i\left(\mu \right)\right\}}_{i=1}^N$ and ${\left\{{\lambda }_j\left(\mu \right)\right\}}_{j=1}^N$ belong to the space $L_{\mathcal{P}}^2(\Gamma )$, while ${\left\{P_{1,i}(t,x)\right\}}_{i=1}^N$ and ${\left\{P_{2,j}(t,x)\right\}}_{j=1}^N$ are corresponding deterministic basis functions in ${\mathcal{V}}^h\left({\Omega }_t\right)$.

Let $\Xi \in \Gamma$ be a set of training samples. We define

$$\begin{array}{c}\hfill \left\{\begin{array}{c}{e}_{p_1}(x,t;\mu ):=p_1(x,t;\mu )-p_{1,k-1}(x,t;\mu ),\hfill \\ e_{p_2}(x,t;\mu ):=p_2(x,t;\mu )-p_{2,k-1}(x,t;\mu ).\hfill \end{array}\right.\end{array}$$

For $\forall v_1,v_2\in V^h(\Omega )$, we then have

$$\begin{array}{c}\hfill \left\{\begin{array}{cc}\hfill & \left({\displaystyle \frac{\partial e_{p_1}}{\partial t}},v_1;\mu \right)+a_1(e_{p_1},v_1;\mu )+q(e_{p_1}-e_{p_2},v_1;\mu )\hfill \\ \hfill & =s_1(v_1;\mu )-\left({\displaystyle \frac{\partial p_{1,k-1}}{\partial t}},v_1;\mu \right)-a_1(p_{1,k-1},v_1;\mu )-q(p_{1,k-1}-p_{2,k-1},v_1;\mu )\hfill \\ \hfill & \left({\displaystyle \frac{\partial e_{p_2}}{\partial t}},v_2;\mu \right)+a_2(e_{p_2},v_2;\mu )+q(e_{p_2}-e_{p_1},v_2;\mu )\hfill \\ \hfill & =s_2(v_2;\mu )-\left({\displaystyle \frac{\partial p_{2,k-1}}{\partial t}},v_2;\mu \right)-a_2(p_{2,k-1},v_2;\mu )-q(p_{2,k-1}-p_{1,k-1},v_2;\mu )\hfill \end{array}\right.\end{array}$$

Let $r_{p_1}^k(v_1;\mu ),r_{p_2}^k(v_2;\mu )\in {\mathcal{V}}^{h,*}\left({\Omega }_t\right)$. We have

$$\begin{array}{cc}\hfill & r_{p_1}^k(v_1;\mu )=\left\{\begin{array}{ccc}\hfill & s_1(v_1;\mu ),\hfill & \hfill k=1,\\ \hfill & s_1(v_1;\mu )-\left({\displaystyle \frac{\partial p_{1,k-1}}{\partial t}},v_1;\mu \right)-a_1(p_{1,k-1},v_1;\mu )-q(p_{1,k-1}-p_{2,k-1},v_1;\mu ),\hfill & \hfill k\ge 2,\end{array}\right.\hfill \end{array}$$

$$\begin{array}{cc}\hfill & r_{p_2}^k(v_2;\mu )=\left\{\begin{array}{ccc}\hfill & s_2(v_2;\mu ),\hfill & \hfill k=1,\\ \hfill & s_2(v_2;\mu )-\left({\displaystyle \frac{\partial p_{2,k-1}}{\partial t}},v_2;\mu \right)-a_2(p_{2,k-1},v_2;\mu )-q(p_{2,k-1}-p_{1,k-1},v_2;\mu ),\hfill & \hfill k\ge 2,\end{array}\right.\hfill \end{array}$$

where ${\mathcal{V}}^{h,*}\left({\Omega }_t\right)$ is the dual space of ${\mathcal{V}}^h\left({\Omega }_t\right)$, and ${\mathcal{M}}^{h,*}\left({\Omega }_t\right)$ is the dual space of ${\mathcal{M}}^h\left({\Omega }_t\right)$. Then we have

$$\begin{array}{c}\hfill \left\{\begin{array}{cc}\hfill & \left({\displaystyle \frac{\partial e_{p_1}}{\partial t}},v_1;\mu \right)+a_1(e_{p_1},v_1;\mu )+q(e_{p_1}-e_{p_2},v_1;\mu )=r_{p_1}^k(v_1;\mu )\hfill \\ \hfill & \left({\displaystyle \frac{\partial e_{p_2}}{\partial t}},v_2;\mu \right)+a_2(e_{p_2},v_2;\mu )+q(e_{p_2}-e_{p_1},v_2;\mu )=r_{p_2}^k(v_2;\mu )\hfill \end{array}\right.\end{array}$$

(12)

By Riesz representation theory, there exist functions ${\widehat{e}}_{p_1}^k\left(\mu \right),{\widehat{e}}_{p_2}^k\left(\mu \right)\in {\mathcal{V}}^h\left({\Omega }_t\right)$, such that

$$\begin{array}{c}\hfill \left\{\begin{array}{c}{({\widehat{e}}_{p_1}^k\left(\mu \right),v_1;\mu )}_{{\mathcal{H}}^1}=r_{p_1}^k(v_1;\mu ),\hfill \\ {({\widehat{e}}_{p_2}^k\left(\mu \right),v_2;\mu )}_{{\mathcal{H}}^1}=r_{p_2}^k(v_2;\mu ).\hfill \end{array}\right.\end{array}$$

(13)

Then we can rewrite Eq. (12) as

$$\begin{array}{c}\hfill \left\{\begin{array}{c}\left({\displaystyle \frac{\partial e_{p_1}}{\partial t}},v_1;\mu \right)+a_1(e_{p_1},v_1;\mu )+q(e_{p_1}-e_{p_2},v_1;\mu )={({\widehat{e}}_{p_1}^k\left(\mu \right),v_1;\mu )}_{{\mathcal{H}}^1},\hfill \\ \left({\displaystyle \frac{\partial e_{p_2}}{\partial t}},v_2;\mu \right)+a_2(e_{p_2},v_2;\mu )+q(e_{p_2}-e_{p_1},v_2;\mu )={({\widehat{e}}_{p_2}^k\left(\mu \right),v_2;\mu )}_{{\mathcal{H}}^1}.\hfill \end{array}\right.\end{array}$$

(14)

Consequently, the dual norm of the residuals $\left\{r_u^k(v_1;\mu )\right\}$, $\left\{r_{\lambda }^k(v_2;\mu )\right\}$, $\left\{r_f^k(\tilde{f};\mu )\right\}$ can be evaluated through the Riesz representation,

$$\begin{array}{c}\hfill \left\{\begin{array}{c}\parallel r_{p_1}^k(v_1;\mu ){\parallel }_{{\left({\mathcal{H}}^1\right)}^{*}}:=\underset{v_1\in {\mathcal{H}}^1}{sup}{\displaystyle \frac{r_{p_1}^k(v_1;\mu )}{\parallel v_1{\parallel }_{{\mathcal{H}}^1}}}={\parallel {\widehat{e}}_{p_1}\left(\mu \right)\parallel }_{{\mathcal{H}}^1},\hfill \\ \parallel r_{p_2}^k(v_2;\mu ){\parallel }_{{\left({\mathcal{H}}^1\right)}^{*}}:=\underset{v_2\in {\mathcal{H}}^1}{sup}{\displaystyle \frac{r_{p_2}^k(v_2;\mu )}{\parallel v_2{\parallel }_{{\mathcal{H}}^1}}}={\parallel {\widehat{e}}_{p_2}\left(\mu \right)\parallel }_{{\mathcal{H}}^1}.\hfill \end{array}\right.\end{array}$$

Next, we define an error estimator for the dual-continuum model (4) by

$${\Delta }_k(x,t;\mu ):=\sqrt{\parallel {\widehat{e}}_{p_1}^k{\left(\mu \right)\parallel }_{{\mathcal{H}}^1}^2+{\parallel {\widehat{e}}_{p_2}^k\left(\mu \right)\parallel }_{{\mathcal{H}}^1}^2}.$$

The computation of the residuals is crucial to the VS procedure. To efficiently compute $\parallel {\widehat{e}}_{p_1}^k{\left(\mu \right)\parallel }_{{\mathcal{H}}^1}$, $\parallel {\widehat{e}}_{p_2}^k{\left(\mu \right)\parallel }_{{\mathcal{H}}^1}$, we apply an offline-online procedure presented in [22,31]. The detailed procedure is presented in Appendix A.

At step k, we choose

$${\mu }_k:=\left\{\begin{array}{cc}\mathrm{choose}\mathrm{randomly}\mathrm{in}\Xi ,\hfill & \hfill k=1,\\ arg\underset{\mu \in \Xi }{max}{\Delta }_k\left(\mu \right)\hfill & \hfill k\ge 2.\end{array}\right.$$

Let $e_{p_1}(x,t;{\mu }_k)$ and $e_{p_2}(x,t;{\mu }_k)$be the solutions of system (12) with $\mu ={\mu }_k$. We take

$$p_{1,k}(t,x)=e_{p_1}(x,t;{\mu }_k),p_{2,k}(t,x)=e_{p_2}(x,t;{\mu }_k).$$

We define the error expansions as

$$e_{p_1}(x,t;\mu )={\xi }_k\left(\mu \right)p_{1,k}(t,x),e_{p_2}(x,t;\mu )={\lambda }_k\left(\mu \right)p_{2,k}(t,x).$$

Next, in the first equation of system (12) we set $v_1=p_{1,k}(x,t_J)$, and in the second equation we set $v_2=p_{2,k}(x,t_J)$, where

$$J:=arg\underset{j=1,\dots ,N_t}{max}\left(r_{p_1}^k(p_{1,k}(x,t_j);{\mu }_k)+r_{p_2}^k(p_{2,k}(x,t_j);{\mu }_k)\right).$$

By (3), (10), (12) and (A.24), we can get the linear system with regard to ${\xi }_k\left(\mu \right)$ and ${\lambda }_k\left(\mu \right)$ as follows

$$\begin{array}{c}\hfill \left\{\begin{array}{c}{L}^{11}\left(\mu \right){\xi }_k\left(\mu \right)+L_s^{12}\left(\mu \right){\lambda }_k\left(\mu \right)=r_{p_1}^k(p_{1,k}(x,t_J);\mu ),\hfill \\ L^{21}\left(\mu \right){\lambda }_k\left(\mu \right)+L^{22}\left(\mu \right){\xi }_k\left(\mu \right)=r_{p_2}^k(p_{2,k}(x,t_J);\mu ).\hfill \end{array}\right.\end{array}$$

(15)

The coefficients are defined as follows:

$$\begin{array}{c}\hfill \left\{\begin{array}{cc}\hfill & L_{11}\left(\mu \right)=({\displaystyle \frac{\partial p_{1,k}(t,x)}{\partial t}},p_{1,k}(x,t_J))+{\displaystyle \sum _{i=1}^{Q_a}}{Q}_1^i\left(\mu \right)a_1^i(p_{1,k}(t,x),p_{1,k}(x,t_J))+q(p_{1,k}(t,x),p_{1,k}(x,t_J))\hfill \\ \hfill & L_{12}\left(\mu \right)=-q(p_{2,k}(t,x),p_{1,k}(x,t_J))\hfill \\ \hfill & L_{21}\left(\mu \right)=({\displaystyle \frac{\partial p_{2,k}(t,x)}{\partial t}},p_{p,k}(x,t_J))+{\displaystyle \sum _{j=1}^{Q_b}}{Q}_2^j\left(\mu \right)a_2^i(p_{2,k}(t,x),p_{21,k}(x,t_J))+q(p_{2,k}(t,x),p_{2,k}(x,t_J))\hfill \\ \hfill & L_{22}\left(\mu \right)=-q(p_{1,k}(t,x),p_{2,k}(x,t_J))\hfill \end{array}\right.\end{array}$$

and

$$\begin{array}{c}\hfill L_g^{11}\left(\gamma \right)=\beta (f^k(t,x),f_J^k\left(x\right)),L_g^{12}\left(\gamma \right)=-({\lambda }^{2k-1}(t,x),f_J^k\left(x\right)),L_g^{13}\left(\gamma \right)=-({\lambda }^{2k}(t,x),f_J^k\left(x\right)).\end{array}$$

The right terms have the following affine expressions:

$$\begin{array}{cc}\hfill r_{p_1}^k(p_{1,k}(x,t_J);\mu )& =\sum _{n=1}^{N_a}{X}_n\left(\mu \right)s_{1,n}\left(p_{1,k}(x,t_J)\right)-\sum _{i=1}^{k-1}{\xi }_i\left(\mu \right)({\displaystyle \frac{\partial p_{1,i}(t,x)}{\partial t}},p_{1,k}(x,t_J))\hfill \\ \hfill & -\sum _{s=1}^{Q_a}\sum _{i=1}^{k-1}{Q}_1^s\left(\mu \right){\xi }_i\left(\mu \right)a_1^s(p_{1,i}(t,x),p_{1,k}(x,t_J))\hfill \\ \hfill & -\sum _{i=1}^{k-1}{\xi }_i\left(\mu \right)q(p_{1,i}(t,x),p_{1,k}(x,t_J))+\sum _{j=1}^{k-1}{\lambda }_j\left(\mu \right)q(p_{2,j}(t,x),p_{1,k}(x,t_J)),\hfill \end{array}$$

$$\begin{array}{cc}\hfill r_{p_2}^k(p_{2,k}(x,t_J);\mu )& =\sum _{m=1}^{M_a}{Y}_m\left(\mu \right)s_{2,m}\left(p_{2,k}(x,t_J)\right)-\sum _{j=1}^{k-1}{\lambda }_j\left(\mu \right)({\displaystyle \frac{\partial p_{2,i}(t,x)}{\partial t}},p_{2,k}(x,t_J))\hfill \\ \hfill & -\sum _{t=1}^{Q_b}\sum _{j=1}^{k-1}{Q}_2^t\left(\mu \right){\lambda }_j\left(\mu \right)a_2^t(p_{2,i}(t,x),p_{2,k}(x,t_J))\hfill \\ \hfill & -\sum _{j=1}^{k-1}{\lambda }_j\left(\mu \right)q(p_{2,j}(t,x),p_{2,k}(x,t_J))+\sum _{i=1}^{k-1}{\xi }_i\left(\mu \right)q(p_{1,i}(t,x),p_{2,k}(x,t_J)).\hfill \end{array}$$

Therefore, the stochastic basis functions ${\xi }_k\left(\mu \right)$, ${\lambda }_k\left(\mu \right)$ can be obtained by system (15). We outline the VS method for stochastic dual-continuum problems in Algorithm 1.

Algorithm 1 The VS method for stochastic dual-continuum problems.

Input: The stochastic dual-continuum problem (3), a set of training samples $\Xi \in \Gamma$, and the error tolerance $\epsilon$. Output: The variable-separation representations $p_{1,N}(x,t;\mu ):={\displaystyle \sum _{i=1}^N}{\xi }_i\left(\mu \right)P_{1,i}(t,x),p_{2,N}(x,t;\mu ):={\displaystyle \sum _{i=1}^N}{\lambda }_i\left(\mu \right)P_{2,i}(t,x)$.   1:   Initialization: Residuals $r_{p_1}^1(\tilde{u};\mu ):=s_1(v_1;\mu )$, $r_{p_2}^1(\tilde{u};\mu ):=s_2(v_2;\mu )$, a random sample ${\mu }_1\in \Xi$;   2:   Calculate the deterministic physical basis functions $P_{1,k}(t,x)=e_{p_1}(x,t;{\mu }_k)$, $P_{2,k}(t,x)=e_{p_2}(x,t;{\mu }_k)$ by solving (12) with $\mu ={\mu }_k$;   3:   Compute the stochastic basis functions ${\xi }^k\left(\mu \right)$, ${\lambda }^k\left(\mu \right)$ by (15);   4:   Update $\Xi$ with $\Xi =\Xi \setminus {\mu }_k$, and take the approximation $p_{1,k}(x,t;\mu ):={\displaystyle \sum _{i=1}^k}{\xi }_i\left(\mu \right)P_{1,i}(t,x)$ and $p_{2,k}(x,t;\mu ):={\displaystyle \sum _{i=1}^k}{\lambda }_i\left(\mu \right)P_{2,i}(t,x)$   5:   k=k+1;   6:   Update the residuals $r_{p_1}^k(v_1;\mu )$ and $r_{p_2}^k(v_2;\mu )$, and compute ${\Delta }_k(x,t;\mu )$;   7:   Take ${\mu }_k$ as ${\mu }_k:=arg\underset{\mu \in \Xi }{max}{\Delta }_k(x,t;\mu )$;   8:   Return Step 2 if ${\mu }_k\ge \epsilon$, otherwise terminate.

5. The Hybrid Model Reduction Method Based on Variable-Separation

Derived from the governing equations of the dual-continuum model, the stochastic and deterministic basis functions used in the VS framework (Section 4) are computed via a traditional finite element method on a fine grid, employing a backward Euler temporal discretization. If the dual-continuum governing equations have strong multiscale features, then we have to use a very fine mesh to resolve the features in all scales. This computation of computing snapshots may be very expensive. To overcome the difficulty and improve the offline computation, we can use the uncoupled Generalized Multiscale Finite Element Method (uncoupled GMsFEM) to compute the basis functions. In this section, we briefly present the uncoupled GMsFEM and the hybrid model reduction method.

5.1. Uncoupled GMsFEM

For the parameterized dual-continuum model, when constructing the multiscale reduced basis space using the uncoupled Generalized Multiscale Finite Element Method (GMsFEM), the multiscale basis functions on each continuum are constructed independently, taking into account only the permeability ${\kappa }_i(x;\mu )$ while neglecting the effect of the interaction term $\sigma$.

Before presenting the framework of uncoupled GMsFEM approximation for the dual-continuum model, we first establish the numerical discretization setting, as illustrated in Figure 1. Let H be the characteristic size of a generic coarse cell. We consider that the spatial domain $\Omega$ to be uniformly partitioned by a coarse mesh ${\mathcal{T}}_H$. The partition ${\mathcal{T}}_H$ is referred to as a coarse grid and a generic element K within it is called a coarse element. Moreover, A corresponding fine-grid partition, denoted by ${\mathcal{T}}_h$, is obtained by refining the coarse mesh ${\mathcal{T}}_H$, where $h>0$ represents the fine mesh size. Let $N_h$ denote the number of elements in ${\mathcal{T}}_H$. The vertices of the coarse grid are denoted by ${\left\{x_j\right\}}_{j=1}^{N_c}$, where $N_c$ is the total number of coarse nodes. We define the neighborhood ${\omega }_j$ as the union of all coarse elements containing $x_j$ in their closure, i.e.,

$${\omega }_j=\cup \{K_i\in {\mathcal{T}}_H|x_j\in \overline{K_i}\}.$$

More specially, on each coarse neighborhood ${\omega }_j$, for each i-th continuum, we first solve the following local eigenvalue problem:

$$\begin{array}{c}\hfill \left\{\begin{array}{c}\begin{array}{cc}\hfill -\mathrm{div}({\kappa }_i\left(x\right)\nabla {\psi }_{\ell ,i}^{\left(j\right)})& ={\lambda }_{i,\ell }\widetilde{\kappa \left(x\right)}{\psi }_{\ell ,i}^{\left(j\right)}\mathrm{in}{\omega }_j,\hfill \\ \hfill {\kappa }_i\left(x\right)·\nabla {\psi }_{\ell ,i}^{\left(j\right)}·n& =0\mathrm{on}\partial {\omega }_j,\hfill \end{array}\hfill \end{array}\right.\end{array}$$

(16)

where $\widetilde{{\kappa }_i\left(x\right)}={\kappa }_i\left(x\right){\displaystyle \sum _{i=1}^{N_c}}{H}^2{|\nabla {\chi }_{j,i}|}^2$. ${\left\{{\chi }_{j,i}\right\}}_{i=1}^{N_c}$ represents the partition of unity functions defined in the local support ${\omega }_j$, which is defined as follows:

$$\left\{\begin{array}{c}\begin{array}{cc}\hfill -\mathrm{div}({\kappa }_i\left(x\right)\nabla {\chi }_{j,i})& =0\forall K\mathrm{in}{\omega }_j\hfill \\ \hfill {\chi }_{j,i}& ={\chi }_{j,i}^0\mathrm{on}\partial K,\hfill \end{array}\hfill \end{array}\right.$$

To construct local multiscale basis functions on ${\omega }_j$ corresponding to the i-th continuum ($i=1,2$), we now solve local spectral problems: find the eigenfunctions ${\psi }_{\ell ,i}^{\left(j\right)})\in V_h\left({\omega }_j\right)$ and ${\lambda }_{i,\ell }\in \mathbb{R}$ such that

$$a_i^{\left(j\right)}\left({\psi }_{\ell ,i}^{\left(j\right)}\right),v_i)={\lambda }_{i,\ell }{s}_i^{\left(j\right)}\left({\psi }_{\ell ,i}^{\left(j\right)}\right),v_i),$$

for all $v_i\in V_h\left({\omega }_j\right)$, where $s_i^{\left(j\right)}\left({\psi }_{\ell ,i}^{\left(j\right)}\right),v_i)$ is defined as follows

$$s_i^{\left(j\right)}(u_i,v_i)={\int }_{{\omega }_j}\widetilde{{\kappa }_i}{u}_i{v}_{i}dx.$$

We denote the basis functions defined on the fine grid as ${\left\{{\varphi }_{j,m}\right\}}_{m=1}^{N_h}$, then the matrix form of eigenvalue problem (16) can be written as

$$\begin{array}{c}\hfill A{\psi }_{\ell ,i}^{\left(j\right)}={\lambda }_{\ell ,i}^{\left(j\right)}S{\psi }_{\ell ,i}^{\left(j\right)}.\end{array}$$

(17)

Here, the definition of matrices A and S are

$${\left(A\right)}_{m,n}={({\kappa }_i\left(x\right)\nabla {\varphi }_{j,m},\nabla {\varphi }_{j,n})}_{L^2\left({\omega }_j\right)},{\left(S\right)}_{m,n}={(\widetilde{{\kappa }_i\left(x\right)}{\varphi }_{j,m},{\varphi }_{j,n})}_{L^2\left({\omega }_j\right)},$$

respectively.

By solving the eigenvalue problem (17), we can get a series of eigenvalues ${\left\{{\lambda }_{\ell ,i}^{\left(j\right)}\right\}}_{i=1}^{N_h}$ and eigenfunctions ${\left\{{\psi }_{\ell ,i}^{\left(j\right)}\right\}}_{i=1}^{N_h}$. Then we choose the $M_i$ lowermost eigenvalues and the corresponding eigenvectors of the eigenvalue problem. Then we can get the local snapshot space

$$V_{\mathrm{ms}}^i\left({\omega }_j\right):=\mathrm{span}\{{\psi }_{\ell ,i}^{\left(j\right)},1\le \ell \le M_j\}.$$

Further, We use the partition of unity functions to paste the snapshot functions and get the multiscale basis function space of i-th continuum

$$V_{\mathrm{ms}}^i(\Omega ):=\sum _{j=1}^{N_c}{V}_{\mathrm{ms}}^i\left({\omega }_j\right)=\mathrm{span}\{{\psi }_{\ell ,i}^{\left(j\right),\mathrm{ms}}:{\psi }_{\ell ,i}^{\left(j\right),\mathrm{ms}}={\chi }_{j,i}{\psi }_{\ell ,i}^{\left(j\right)},1\le j\le N_c,1\le \ell \le M_i\}.$$

Finally, the global multiscale corresponding to the dual continuum model as

$${\mathit{V}}_{\mathrm{ms}}(\Omega )=V_{\mathrm{ms}}^1(\Omega )\times V_{\mathrm{ms}}^2(\Omega ).$$

The multiscale basis functions can be repeatedly used online. We can use a matrix $R_{re}^i$ to store these basis functions, i.e.,

$$R_{re}^i=[{\psi }_{1,i}^{\left(1\right),\mathrm{ms}},\dots ,{\psi }_{M_i,i}^{\left(1\right),\mathrm{ms}},\dots ,{\psi }_{1,i}^{\left(N_c\right),\mathrm{ms}},\dots ,{\psi }_{M_i,i}^{\left(N_c\right),\mathrm{ms}}].$$

Then define the global matrix in the dual continuum model as

$$R_{re}=[R_{re}^1,R_{re}^2].$$

At the online stage, we use the precomputed multiscale basis functions to compute dual-continuum problem on the coarse grid.

5.2. The Hybrid Model Reduction Method

In this subsection, we will construct the hybrid model reduction method by combing the variable-separation method and uncoupled GMsFEM, which gives the solution in the following form

$$\left\{\begin{array}{c}\begin{array}{cc}\hfill & p_1^{vs-GMsFEM}(x,t;\mu ):=\sum _{i=1}^N{\xi }_i\left(\mu \right)P_{1,i}^{gms}(t,x)\hfill \\ \hfill & p_2^{vs-GMsFEM}(x,t;\mu ):=\sum _{i=1}^N{\lambda }_i\left(\mu \right)P_{2,i}^{gms}(t,x).\hfill \end{array}\hfill \end{array}\right.$$

(18)

In the separation form (18), the physical basis functions ${\left\{P_{1,i}^{gms}(t,x)\right\}}_{i=1}^N$ and ${\left\{P_{2,i}^{gms}(t,x)\right\}}_{i=1}^N$ are learned using the uncoupled GMsFEM.

$$\begin{array}{c}\hfill \left\{\begin{array}{c}\begin{array}{cc}& {\left(R_{re}^1\right)}^T(M_{1,h}+\Delta t{K}_{1,h})R_{re}^1{P}_{1,h}^{gms,k+1}+\Delta t[{\left(R_{re}^1\right)}^T{M}_{2,h}{R}_{re}^1{P}_{1,h}^{gms,k+1}-\hfill \\ & {\left(R_{re}^1\right)}^T{M}_{2,h}{R}_{re}^2{P}_{2,h}^{gms,k+1}]=\Delta t{\left(R_{re}^1\right)}^T{B}_1^{k+1}+{\left(R_{re}^1\right)}^T{M}_{1,h}{R}_{re}^1{P}_{1,h}^{gms,k}\hfill \\ & {\left(R_{re}^2\right)}^T(M_{1,h}+\Delta t{K}_{2,h})R_{re}^2{P}_{2,h}^{gms,k+1}+\Delta t[{\left(R_{re}^2\right)}^T{M}_{2,h}{R}_{re}^2{P}_{2,h}^{gms,k+1}-\hfill \\ & {\left(R_{re}^2\right)}^T{M}_{2,h}{R}_{re}^1{P}_{1,h}^{gms,k+1}]=\Delta t{\left(R_{re}^2\right)}^T{B}_2^{k+1}+{\left(R_{re}^2\right)}^T{M}_{2,h}{R}_{re}^2{P}_{2,h}^{gms,k}.\hfill \end{array}\hfill \end{array}\right.\end{array}$$

(19)

To simplify the symbols, we denote the matrices by

$$\left\{\begin{array}{c}\begin{array}{cc}\hfill A_1& ={\left(R_{re}^1\right)}^T(M_{1,h}+\Delta t{K}_{1,h})R_{re}^1,\hfill \\ \hfill A_2& ={\left(R_{re}^2\right)}^T(M_{1,h}+\Delta t{K}_{2,h})R_{re}^2,\hfill \\ \hfill C_{11}& ={\left(R_{re}^1\right)}^T{M}_{2,h}{R}_{re}^1,\hfill \\ \hfill C_{12}& ={\left(R_{re}^1\right)}^T{M}_{2,h}{R}_{re}^2,\hfill \\ \hfill C_{22}& ={\left(R_{re}^2\right)}^T{M}_{2,h}{R}_{re}^2,\hfill \\ \hfill C_{21}& ={\left(R_{re}^1\right)}^T{M}_{2,h}{R}_{re}^1,\hfill \\ \hfill F_1^{gms,k+1}\left(\mu \right)& =\Delta t{\left(R_{re}^1\right)}^T{B}_1^{k+1}+{\left(R_{re}^1\right)}^T{M}_{1,h}{R}_{re}^1{P}_{1,h}^{gms,k},\hfill \\ \hfill F_1^{gms,k+1}\left(\mu \right)& =\Delta t{\left(R_{re}^2\right)}^T{B}_2^{k+1}+{\left(R_{re}^2\right)}^T{M}_{2,h}{R}_{re}^2{P}_{2,h}^{gms,k}.\hfill \end{array}\hfill \end{array}\right.$$

Combining these equations gives the matrix system

$$\left[\begin{array}{cc}{A}_1+\Delta t{C}_{11}& \Delta t{C}_{12}\\ \Delta t{C}_{22}& A_2+\Delta t{C}_{22}\end{array}\right]\left[\begin{array}{c}{P}_{1,h}^{gms,k+1}\left(\mu \right)\\ P_{2,h}^{gms,k+1}\left(\mu \right)\end{array}\right]=\left[\begin{array}{c}{F}_1^{gms,k+1}\left(\mu \right)\\ F_2^{gms,k+1}\left(\mu \right)\end{array}\right].$$

(20)

Here, the terms $P_{1,h}^{gms,k+1}\left(\mu \right)$ and $P_{2,h}^{gms,k+1}\left(\mu \right)$ represent the discrete numerical approximations for $p_1,p_2$ at the time step $(k+1)$, obtained using the uncoupled GMsFEM.

In Algorithm 2, we outline the main steps for the hybrid model reduction method based on the VS method.

Algorithm 2 Hybrid model reduction method for the stochastic dual continuum model

Input: The stochastic dual continuum model (4) Output: The hybrid reduced model solutions $p_1^{vs-GMsFEM}(x,t;\mu ):={\displaystyle \sum _{i=1}^N}{\xi }_i\left(\mu \right)P_{1,i}^{gms}(t,x),p_2^{vs-GMsFEM}(x,t;\mu ):={\displaystyle \sum _{i=1}^N}{\lambda }_i\left(\mu \right)P_{2,i}^{gms}(t,x)$.   1:   Initialization: Residuals $r_{p_1}^1(\tilde{u};\mu ):=s_1(v_1;\mu )$, $r_{p_2}^1(\tilde{u};\mu ):=s_2(v_2;\mu )$, a random sample ${\mu }_1\in \Xi$;   2:   Run the algorithm corresponding to system (20) similar to Algorithm 1 to obtain the physical basis functions $\left\{P_{1,k}^{gms}(t,x)\right\}$ and $\left\{P_{2,k}^{gms}(t,x)\right\}$ with $\mu ={\mu }_k$;   3:   Compute the stochastic basis functions ${\xi }^k\left(\mu \right)$, ${\lambda }^k\left(\mu \right)$ similar to (15);   4:   Update $\Xi$ with $\Xi =\Xi \setminus {\mu }_k$, and take the approximation $p_1^{hb}(x,t;\mu ):={\displaystyle \sum _{i=1}^k}{\xi }_i\left(\mu \right)P_{1,i}^{gms}(t,x)$ and $p_2^{hb}(x,t;\mu ):={\displaystyle \sum _{i=1}^k}{\lambda }_i\left(\mu \right)P_{2,i}^{gms}(t,x)$   5:   k=k+1;   6:   Update the residuals $r_{p_1}^k(v_1;\mu )$ and $r_{p_2}^k(v_2;\mu )$, and compute ${\Delta }_k(x,t;\mu )$;   7:   Take ${\mu }_k$ as ${\mu }_k:=arg\underset{\mu \in \Xi }{max}{\Delta }_k(x,t;\mu )$;   8:   Return Step 2 if ${\mu }_k\ge \epsilon$, otherwise terminate.

6. Numerical Results

In this section, we present two numerical examples to illustrate the applicability of the proposed reduced basis methods for solving the stochastic dual continuum models. In Section 6.1, we consider an example to illustrate performance of the hybrid model reduction method for solving the dual continuum model. In Section 6.2, we use the proposed model reduction method to a dual continuum model with high-dimensional random inputs. Before presenting the individual examples, we describe the computational domain, that is, spatial domain $D={[0,1]}^2$. The time interval is from 0 to 1.

Let $p_1^{MR}$ and $p_2^{MR}$ be the solutions of the reduced basis methods, and $(p_{1,h},p_{2,h})$ be the reference solutions, which are solved by FEM on a fine grid. Then the relative mean errors in the weighted $L^2$ norm are defined by

$$\begin{array}{c}\hfill \begin{array}{c}\hfill e_j={\displaystyle \frac{1}{N}}\sum _{i=1}^N{\displaystyle \frac{\parallel p_j^{MR}(x,{\mu }_i)-p_j^h(x,{\mu }_i){\parallel }_{L^2}}{\parallel p_j^h(x,{\mu }_i){\parallel }_{L^2}}},\end{array}\end{array}$$

(21)

where $j=1,2$, ${\parallel v\parallel }_{L^2}^2={\int }_D{v}^{2}dx$, and N is the number of random samples.

6.1. Experiment 1

To demonstrate the efficacy of the hybrid model reduction method, we first consider the dual-continuum with random inputs as follows.

$$\begin{array}{c}\hfill \left\{\begin{array}{c}\begin{array}{cc}& {\displaystyle \frac{\partial p_1}{\partial t}}-\mathrm{div}({\kappa }_1(x,\mu )\nabla p_1)+\sigma (x,\mu )(p_1-p_2)=f_1(t,x,\mu ),\hfill \\ & {\displaystyle \frac{\partial p_2}{\partial t}}-\mathrm{div}({\kappa }_2(x,\mu )\nabla p_2)+\sigma (x,\mu )(p_2-p_1)=f_2(t,x,\mu ),\hfill \end{array}\hfill \end{array}\right.\end{array}$$

(22)

where the diffusion coefficients, the transfer coefficient, and the source terms are defined by

$$\left\{\begin{array}{c}\begin{array}{cc}\hfill & {\kappa }_1(x,\mu )=(1+{\mu }^2)(1+x_1)(1+x_2)+(1+\cos \left(\mu \right)+{\mu }^4){\kappa }_{12}\left(x\right),\hfill \\ \hfill & {\kappa }_2(x,\mu )=(1+2{\mu }^2)(1+x_1^2)(1+x_2^2)+(1+\sin \left(\mu \right)+{\mu }^2){\kappa }_{22}\left(x\right),\hfill \\ \hfill & \sigma (x,\mu )=(1+\mu )(1+x_1)+(1+{\mu }^2)x_1^2{x}_2^2,\hfill \\ \hfill & f_1(t,x,\mu )=2{\mu }^2(1-x_1-x_2-t)+2{\mu }^4{t}^2{x}_1^2{x}_2^2(1-x_1-x_2-t)\hfill \\ \hfill & -2{\mu }^2(1+{\mu }^2)\left((1+2x_1^3)(1+x_2^4)t+(1+x_1)(1+2x_2+t)\right)\hfill \\ \hfill & f_2(t,x,\mu )=0.\hfill \end{array}\hfill \end{array}\right.$$

Here $x=(x_1,x_2)\in \Omega$ and the random variable $\mathit{\xi }$ obeys the beta distribution, i.e., $\mathit{\xi }\sim \mathrm{Beta}(1,1)$. ${\kappa }_{12}\left(x\right)$ and ${\kappa }_{22}\left(x\right)$ are both high-contrast functions independent of $\mu$, whose maps are depicted in Figure 2.

The boundary conditions and initial condition are set as

$$\left\{\begin{array}{c}\begin{array}{cc}\hfill & g_1\left(x\right)=x_1^2+x_2^2,\hfill \\ \hfill & g_2\left(x\right)={(x-1)}_1^2+{(x-1)}_2^2,\hfill \\ \hfill & p_1(x,t_0)=x_1,\hfill \\ \hfill & p_2(x,t_0)=x_2,\hfill \end{array}\hfill \end{array}\right.$$

respectively. The partition of spatial domain is defined on $120\times 120$ uniform grid to compute the reference solution. The model reduction solutions are computed on $10\times 10$ coarse mesh. The time step is $△t=0.02$ up to the final simulation time $T=1$. For offline computation in this example, we randomly choose 200 random samples for the training set $\Xi$ and select 5 optimal samples.

Figure 3 shows the mean profiles of the solutions $p_i(x,t;\mu ),i=1,2$ for the reference, GMsFEM and hybrid model reduction methods. The results demonstrate that the mean profiles from both reduction methods are nearly identical to the reference.

To evaluate the performance of the hybrid model reduction method, we compute the average relative $L^2$ error over an ensemble of 1000 randomly selected parameter samples.

The numerical setup is configured as follows: the number of local basis functions per coarse element is fixed at $L=5$, and the number of variable-separation terms is also set to five. To systematically investigate the impact of the coarse mesh resolution on the approximation accuracy, we conduct a series of experiments with varying coarse mesh sizes, specifically $H=\{1/2,1/4,1/8,1/10,1/12\}$.

Table 1. The relative $L^2$ errors with different coarse mesh size H for $p_1(x,\mathit{\xi })$ and $p_2(x,\mathit{\xi })$.

$Coarsemeshsize$	$H=1/2$	$H=1/4$	$H=1/8$	$H=1/10$	$H=1/12$	$H=1/15$
$e_{p_1}^{fg}$	5.0389E-01	1.2556E-01	1.3176E-01	3.2681E-03	3.4457E-03	1.2881E-03
$e_{p_1}^{fv}$	5.0337E-01	1.2563E-01	1.3191E-01	4.6316E-03	4.1590E-02	5.9059E-03
$e_{p_2}^{fg}$	2.2118E-01	1.3833E-01	8.7244E-02	2.5336E-03	1.8137E-03	1.4234E-03
$e_{p_2}^{fv}$	2.2208E-01	1.4099E-01	8.6095E-02	7.0355E-03	8.9305E-03	5.4730E-03

Table 1. The relative $L^2$ errors with different coarse mesh size H for $p_1(x,\mathit{\xi })$ and $p_2(x,\mathit{\xi })$.

$Coarsemeshsize$	$H=1/2$	$H=1/4$	$H=1/8$	$H=1/10$	$H=1/12$	$H=1/15$
$e_{p_1}^{fg}$	5.0389E-01	1.2556E-01	1.3176E-01	3.2681E-03	3.4457E-03	1.2881E-03
$e_{p_1}^{fv}$	5.0337E-01	1.2563E-01	1.3191E-01	4.6316E-03	4.1590E-02	5.9059E-03
$e_{p_2}^{fg}$	2.2118E-01	1.3833E-01	8.7244E-02	2.5336E-03	1.8137E-03	1.4234E-03
$e_{p_2}^{fv}$	2.2208E-01	1.4099E-01	8.6095E-02	7.0355E-03	8.9305E-03	5.4730E-03

Figure 4. The average relative errors corresponding to the time t.

Figure 4. The average relative errors corresponding to the time t.

6.2. Experiment 2

In this subsection, we consider the dual-continuum with random inputs as follows.

$$\begin{array}{c}\hfill \left\{\begin{array}{c}\begin{array}{cc}& {\displaystyle \frac{\partial p_1}{\partial t}}-\mathrm{div}({\kappa }_1(x,\mu )\nabla p_1)+\sigma (x,\mu )(p_1-p_2)=f_1(t,x,\mu ),\hfill \\ & {\displaystyle \frac{\partial p_2}{\partial t}}-\mathrm{div}({\kappa }_2(x,\mu )\nabla p_2)+\sigma (x,\mu )(p_2-p_1)=f_2(t,x,\mu ),\hfill \end{array}\hfill \end{array}\right.\end{array}$$

(23)

where the diffusion coefficients, the transfer coefficient, and the source terms are defined by

$$\left\{\begin{array}{c}\begin{array}{cc}\hfill & {\kappa }_1(x,\mu )=(1+2{\mu }^2)(1+\exp (x_1+x_2))+(2+\sin \left(\mu \right)){\kappa }_{12}\left(x\right),\hfill \\ \hfill & {\kappa }_2(x,\mu )=(2+\exp \left(\mu \right))(2+\sin \left(2x_1{x}_2\right))+(1+\sin \left(2\mu \right)+{\mu }^2){\kappa }_{22}\left(x\right),\hfill \\ \hfill & \sigma (x,\mu )=(1+2\mu )(1+x_1)+(3+{\mu }^2)x_1^2{x}_2^2,\hfill \\ \hfill & f_1(t,x,\mu )=(\mu +{\mu }^2)(1+x_1+x_2+t)+({\mu }^4-{\mu }^2)t^2{x}_1^2{x}_2^2(1+x_1+x_2+t),\hfill \\ \hfill & f_2(t,x,\mu )=(\mu -{\mu }^3)(x_1+x_2+2)+({\mu }^2+{\mu }^4)2x_1{x}_2(x_1+x_2+2).\hfill \end{array}\hfill \end{array}\right.$$

Here $x=(x_1,x_2)\in \Omega$ and the random variable $\mathit{\mu }$ obeys the beta distribution, i.e., $\mathit{\mu }\sim \mathrm{Beta}(2,2)$. ${\kappa }_{12}\left(x\right)$ and ${\kappa }_{22}\left(x\right)$ are both high-contrast functions independent of $\mu$, whose maps are depicted in Figure 5.

The conductivity values in the background both are fixed to be 1, while the conductivity values in the channels are 10000. In this experiment, we randomly take 2000 random samples to test the proposed iterative method. The boundary conditions and initial condition are set as

$$\left\{\begin{array}{c}\begin{array}{cc}\hfill & g_1\left(x\right)=2-x_1^2-x_2^2,\hfill \\ \hfill & g_2\left(x\right)=2+x_1+x_2,\hfill \\ \hfill & p_1(x,t_0)=2-x_1^2-x_2^2,\hfill \\ \hfill & p_2(x,t_0)=2+x_1+x_2,\hfill \end{array}\hfill \end{array}\right.$$

respectively. In this experiment, the solutions of the multi-continuum model for the three methods: FEM on $120\times 120$ fine grid (reference solutions), GMsFEM on $12\times 12$ coarse mesh size and the proposed model reduction method.

To evaluate the approximation accuracy for the hybrid model reduction method, we randomly selected 2000 parameter samples and compute the average relative $L^2$-errors. Using a fixed number of separated terms ($N=10$) and five local basis functions in GMsFEM, the mean and the variance of variables $p_1(t,x;\mu )$ and $p_2(t,x;\mu )$ are illustrated at a random time layer in Figure 6 and Figure 7, respectively.

The figure demonstrates that the solutions for both $p_1(t,x,\xi )$ and $p_2(t,x,\xi )$ closely match their respective reference solutions.

To systematically assess the impact of the coarse mesh resolution on accuracy, we conduct experiments using progressively refined meshes with sizes $H=\{1/2,1/4,1/8,1/12,1/24\}$. We compare the relative $L^2$ errors of all methods in Table 2.

These results, presented in the table, demonstrate that the approximation accuracy of both the hybrid method and GMsFEM improves—that is, their respective errors decrease—as the coarse mesh is refined.

Figure 7. The variance of $p_1(x,\mu )$ (1st row) and $p_2(x,\mu )$ (2nd row) by different methods.

Figure 7. The variance of $p_1(x,\mu )$ (1st row) and $p_2(x,\mu )$ (2nd row) by different methods.

Figure 8. The relative $L^2$-error of numerical solutions: $p_1(t,x;\mu )$(left) and $p_2(t,x;\mu )$(right) with different number of separated terms.

Figure 8. The relative $L^2$-error of numerical solutions: $p_1(t,x;\mu )$(left) and $p_2(t,x;\mu )$(right) with different number of separated terms.

7. Conclusions

In this paper, we present a hybrid model reduction method for solving stochastic dual-continuum models. In this work, In this work, the possible uncertainties considered originate in the PDE coefficient, the transfer term, and the source term. Directly solving these equations on a fine grid presents a significant computational challenge, particularly in many-query scenarios. To dramatically improve efficiency, our method integrates the uncoupled Generalized Multiscale Finite Element Method (GMsFEM) with a low-rank variable-separation technique to construct compact representations of the solutions. This hybrid approach enables the computation of a large number of stochastic realizations at minimal cost. In the numerical experiments section, we apply the proposed model reduction method to dual-continuum models with random inputs. The results confirm that our hybrid method achieves an effective trade-off between computational efficiency and accuracy.