Multiple Continuum Model

Two forms of multiple continuum models are available in PFLOTRAN, one based on thermal conduction and the other incorporating multicomponent reactive transport. Both models only account for diffusion in the secondary continua (matrix) modeled as a disconnected one-dimensional domains referred to as the DCDM (Dual Continuum Disconnected Matrix) model (Lichtner, 2000), but include advection in the primary continuum.

In contrast to the approach used by TOUGH’s MINC method in which primary and secondary continua are treated as a single system of equations solved simultaneously, in the PFLOTRAN implementation the primary and secondary continua are solved separately treating the secondary continua as a 1D system of equations. Parallelization is highly efficient because the secondary continuum equations are embarrassingly parallel. The results should be identical between the two implementations. Details of the implementation of the DCDM model in PFLOTRAN and its scalability on multiple processor cores are presented in Lichtner and Karra (2014).

Thermal Conduction

A thermal conduction model employing a multiple continuum model has been added to modes MPHASE and TH. The formulation is based on Pruess and Narasimhan (1985) using a integrated finite volume approach to develop equations for fracture (primary continuum) and matrix (secondary continua) temperatures \(T_{\alpha}\) and \(T_{\beta}\), respectively, with fracture volume fraction denoted by \(\epsilon_{{\alpha}}\). In what follows the matrix porosity is assumed to be zero.

In terms of partial differential equations the heat conservation equations may be written as

(1)\[\begin{split}\frac{{{\partial}}}{{{\partial}}t} \epsilon_{{\alpha}}\Big[\varphi_{{\alpha}}\rho_{{\alpha}}U_{{\alpha}}+ (1-\varphi_{{\alpha}}) \rho_r C_r T_{{\alpha}}\Big] &+ {\boldsymbol{\nabla}}\cdot \Big({\boldsymbol{q}}_{{\alpha}}\rho_{{\alpha}}H_{{\alpha}}- \kappa_a{\boldsymbol{\nabla}}T_{{\alpha}}\Big) \nonumber\\ &= -\sum_\beta A_{{{\alpha}}{{\beta}}} \kappa_{{{\alpha}}{{\beta}}}\frac{{{\partial}}T_{{\beta}}}{{{\partial}}n},\end{split}\]

and

(2)\[\frac{{{\partial}}}{{{\partial}}t} \rho_r C_r T_{{\beta}}+ \frac{{{\partial}}}{{{\partial}}\xi} \Big(-\kappa_{{\beta}}\frac{{{\partial}}T_{{\beta}}}{{{\partial}}\xi}\Big) = 0,\]

for fracture and matrix temperatures \(T_{{\alpha}}\) and \(T_{{\beta}}\), respectively, where \(\xi\) represents the matrix coordinate assumed to be an effective 1D domain. The boundary condition

(3)\[T_{{\beta}}(\xi_M,\,t\,|\,{\boldsymbol{r}}) = T_{{\alpha}}({\boldsymbol{r}},\,t),\]

is imposed at the fracture-matrix interface, where \(\xi_M\) denotes the outer boundary of the matrix.

Using the control volume configuration shown in Figure 10, with fracture nodes designated by the subscript \(n\) and matrix nodes by \(m\), the integrated finite volume form of the heat transport equation for the \(n\)th fracture control volume is given by

(4)\[\begin{split}&\Big[\varphi_{{\alpha}}\Big(\big(\rho_{{\alpha}}U_{{\alpha}}\big)_n^{k+1} - \big(\rho_{{\alpha}}U_{{\alpha}}\big)_n^{k}\Big) + (1-\varphi_{{\alpha}})\Big(\big(\rho_r C_r T_{{\alpha}}\big)_n^{k+1} - \big(\rho_r C_r T_{{\alpha}}\big)_n^{k}\Big)\Big] \frac{V_n^{{\alpha}}}{\Delta t} \nonumber\\ &\qquad + \sum_{n'} \Big[\big(q_{{\alpha}}\rho_{{\alpha}}H_{{\alpha}}\big)_{nn'} + \frac{\kappa_{nn'}^{{\alpha}}}{d_n+d_{n'}}\big(T_{{{\alpha}}n} - T_{{{\alpha}}n'}\big) \Big] A_{nn'}^{{\alpha}}\nonumber\\ &\qquad + \sum_{l = 1}^{N_{{\beta}}}\frac{\kappa_{nM}^{{{\alpha}}{{\beta}}_l}}{d_n+d_{M}}\big(T_{{{\alpha}}n}-T_{{{\beta}}_l M}\big) A_{nM}^{{{\beta}}_l} = 0,\end{split}\]

where \(V_n^{{\alpha}}\) denotes the fracture volume, and

(5)\[\begin{split}\Big(\big(\rho_r C_r T_{{\beta}}\big)_m^{k+1} - \big(\rho_r C_r T_{{\beta}}\big)_m^{k}\Big) \frac{V_m^{{\beta}}}{\Delta t} &+ \sum_{m'} \frac{\kappa_{mm'}^{{\beta}}}{d_m+d_{m'}}\big(T_{{{\beta}}m} - T_{{{\beta}}m'}\big) A_{mm'}^{{\beta}}\nonumber\\ &+ \delta_{mM}^{}\frac{\kappa_{nM}^{{{\alpha}}{{\beta}}}}{d_n+d_{M}}\big(T_{{{\alpha}}n} - T_{{{\beta}}M}\big) A_{nM}^{{\beta}}= 0,\end{split}\]

for the \(m\)th matrix node with volume \(V_m^{{\beta}}\). The matrix node designated by \(M\) refers to the outer most node in contact with the fracture (see Figure [fnodestruct]). More than one type of matrix geometry is included in the above equations as indicated by the sum over \(l\) in Eqn. (4), where \(N_{{\beta}}\) denotes the number of different secondary continua. However, it should be noted that the current implementation in PFLOTRAN only allows coupling to a single secondary continuum \((N_{{\beta}}=1)\). The fracture volume \(V_n^{{\alpha}}\) is related to the REV volume \(V_n\) by

(6)\[\epsilon_{{\alpha}}= \frac{V_n^{{\alpha}}}{V_n}.\]

Thermal conductivity at the interface between two control volumes is calculated using the harmonic average

(7)\[\kappa_{ll'} = \frac{\kappa_l \kappa_{l'}(d_l+d_{l'})}{d_l \kappa_{l'}+d_{l'}\kappa_l}.\]

(8)\[\qquad \bigg|\quad\mathop{\bullet}_{\ \ \ \, \displaystyle 1 \ {{\beta}}}\quad\bigg| \qquad \cdots \qquad \bigg|\quad\mathop{\bullet}_{\ \ \ \, \displaystyle l \ {{\beta}}}\quad\bigg| \qquad \cdots \qquad \bigg|\quad\mathop{\bullet}_{\ \ \ \, \displaystyle M \ {{\beta}}}\quad \bigg|\quad\mathop{\bullet}_{\ \ \ \, \displaystyle n \ {{\alpha}}}\quad\bigg|\nonumber\]

For better convergence uniform logarithmic spacing is used for the matrix nodes

(9)\[\Delta \xi_m = \rho \,\Delta \xi_{m-1},\]

specifying \(\Delta\xi_M\) and \(l_M\) for the outer most matrix node and matrix block size, respectively. The factor \(\rho\) is determined from the constraint

(10)\[l_M = 2\sum_{m=1}^{M} \Delta \xi_m,\]

which requires that \(\rho\) satisfy the equation

(11)\[\frac{l_M}{2\Delta \xi_1} = \frac{\rho^M-1}{\rho-1},\]

with the inner and outer grid spacing related by

(12)\[\Delta\xi_M = \rho^{M-1} \Delta \xi_1.\]

Control volumes in DCDM multiple continuum model with fracture aperture \(2\delta\) and matrix block size \(d\).

According to the geometry in Figure 10 assuming a 3D orthogonal set of fractures,

(13)\[V_n = (d+2\delta)^3,\]

and

(14)\[V_n^{{\alpha}}= (d+2\delta)^3 - d^3,\]

giving

(15)\[\begin{split}\epsilon_{\alpha} &= 1-\frac{d^3}{(d+2\delta)^3} = 1-\Bigg(\dfrac{1}{1+\dfrac{2\delta}{d}}\Bigg)^3,\\ & ~\simeq~ \frac{6\delta}{d}.\end{split}\]

The fracture aperture \(2\delta\) is found to be in terms of \(\epsilon_{{\alpha}}\) and \(d\)

(16)\[2\delta = d \Big(\frac{1}{(1-\epsilon_{{\alpha}})^{1/3}} -1\Big).\]

A list of different sub-continua geometries and parameters implemented in PFLOTRAN is given in Table [tdcdmgeom]. Different independent and dependent parameters for the nested cube geometry are listed in Table [tnestedcube]. The interfacial area \(A_{nn'}^{{\alpha}}\) between fracture control volumes is equal to \(\Delta y \Delta z\), \(\Delta z \Delta x\), \(\Delta x \Delta y\) for \(x\), \(y\), and \(z\) directions, respectively.

In the case of nested cubes there are four possible parameters \((\epsilon_{{\alpha}}, \, 2\delta, \, l_m,\, l_f)\), where \(l_m\) denotes the matrix block size and \(l_f\) refers to the fracture spacing, two of which are independent.

The fracture-matrix interfacial area \(A_{nM}\) per unit volume is equal to

(17)\[A_{nM}^{{\beta}}= \frac{{{\mathcal N}}_{{\beta}}}{V} A_{{\beta}}^0,\]

where the number density \({{\mathcal N}}_{{\beta}}/V\) of secondary continua of type \({{\beta}}\) is equal to

(18)\[\frac{{{\mathcal N}}_{{\beta}}}{V} = \frac{1}{V} \frac{V_{{\beta}}}{V_{{\beta}}^0} = \frac{\epsilon_{{\beta}}}{V_{{\beta}}^0},\]

and \(A_{{\beta}}^0\) and \(V_{{\beta}}^0\) refer to the area and volume of each geometric type as listed in Table [tdcdmgeom].

Geometry	Area \(A_{{\beta}}^0\)	Volume \(V_{{\beta}}^0\)
Slab	\(A\)	\(A l\)
Nested Cubes	\(6d^2\)	\(d^3\)
Nested Spheres	\(4 \pi R^2\)	\(4/3 \pi R^3\)

Table: DCDM geometric parameters.

The primary-secondary coupling term can then be written in the form

(19)\[\sum_{{\beta}}\frac{\kappa_{nM}^{{{\alpha}}{{\beta}}}}{d_n+d_{M}}\big(T_n^{{\alpha}}-T_{M}^{{\beta}}\big) A_{nM}^{{\beta}}= V_n \sum_{{\beta}}\frac{\epsilon_{{\beta}}\kappa_{nM}^{{{\alpha}}{{\beta}}}}{d_n+d_{M}}\big(T_n^{{\alpha}}-T_{M}^{{\beta}}\big) \frac{A_{{\beta}}^0}{V_{{\beta}}^0}.\]

Independent		Dependent
\(\epsilon_{{\alpha}}\)	\(l_f\)	\(2\delta = l_f - l_m\)	\(l_m = l_f(1-\epsilon_{{\alpha}})^{1/3}\)
\(\epsilon_{{\alpha}}\)	\(l_m\)	\(2\delta = l_f - l_m\)	\(l_f = l_m(1-\epsilon_{{\alpha}})^{-1/3}\)
\(2\delta\)	\(l_f\)	\(\epsilon_{{\alpha}}= 1-(l_m/l_f)^3\)	\(l_m = l_f - 2\delta\)
\(2\delta\)	\(l_m\)	\(\epsilon_{{\alpha}}= 1-(l_m/_f)^3\)	\(l_f = l_m + 2\delta\)
\(2\delta\)	\({\epsilon}_{\alpha}\)	\(l_m = 2\delta \Big(\dfrac{1}{(1-\epsilon_{{\alpha}})^{1/3}}-1\Big)^{-1}\)	\(l_m = l-2\delta\)

Table: Independent and dependent nested cube parameters.

Reactive Transport Dual Continuum Model

The implementation of a dual continuum model for reactive transport is based on the DCDM model. The primary continuum equations have the form

(20)\[\frac{\partial}{\partial t} \Big(\epsilon_\alpha \varphi_\alpha \sum_p s_p^\alpha \Psi_{jp}^\alpha\Big) + \nabla\cdot\sum_p \epsilon_\alpha \Omega_{jp}^\alpha = -\sum_{p\beta} A_{\alpha\beta} \Omega_{jp}^{\alpha\beta} - \epsilon_\alpha \sum_m \nu_{jm}^{} I_{mp}^\alpha - \epsilon_\alpha \frac{\partial S_{jp}^\alpha}{\partial t},\]

where now an additional term appears on the right-hand side representing mass transfer between primary and secondary continua with

(21)\[\Omega_{jp}^{\alpha\beta}(r,\,t) = \Omega_{jp}^\beta (\xi_{\alpha\beta},\,t|r).\]

The secondary continuum mass conservation equations have a similar form but without the factor \(\epsilon_\alpha\) and the coupling term. Imposition of symmetry at the boundary of the secondary continuum leads to the equation

(22)\[\frac{\partial}{\partial t} \Big(\varphi_\beta \sum_p s_p^\beta\Psi_{jp}^\beta\Big) + \nabla_\xi\cdot\sum_p \Omega_{jp}^\beta = - \sum_m \nu_{jm}^{} I_{mp}^\beta - \frac{\partial S_{jp}^\beta}{\partial t},\]

where the gradient operator \(\nabla_\xi\) refers to the effective one-dimensional secondary continuum geometry. Similar considerations apply to mass and heat flow for primary and secondary continuum conservation equations.