# 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.

## 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

and

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

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

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

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

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

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

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

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

with the inner and outer grid spacing related by

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

and

giving

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

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

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

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

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

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

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

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.