Hybrid Transport Modelling Approach: UMT3D

Transport of mass within the porous matrix is also simulated with a continuum approach, i.e., with the standardised three-dimensional multi-species transport model MT3DMS (Zheng & Wang 1999). MT3DMS was selected due to its compatibility to the MODFLOW model. In contrast, mass transport within the network of tubes is modelled discretely with a one-dimensional advective transport model similar to that presented by Birk (2001) in CAVE. Coupling between the continuum transport model and the transport in the conduits is achieved using the two-step method SNIA (Sequential Non-Iterative Approach, compare Chapter 2.4.3 and Walter et al. 1994; Steefel & MacQuarrie 1996). The transport of solutes in the continuum is solved first, followed by a second step during which the solute concentrations in the conduit system are updated and the solute mass is transported. The resulting coupled hybrid transport model was named UMT3D (three-dimensional Underground Mine Transport model). For the clarity of the difference between the UMT3D and CAVE model, CAVE, primarily as the hybrid flow model for the conduits and continuum, can also calculate mass transport in the conduits but not in the continuum. However, the option of conduit transport calculation in CAVE is not utilised when flow data are calculated for the UMT3D model.

The general transport equation (Zheng & Wang 1999) was extended to include a further sink/source term, i.e., ξ to couple MT3DMS to the conduit transport model in UMT3D. Thus, the three-dimensional transport equation in the continuum is:

( ) ∑

Cm (M L^-3): aqueous concentration of component k in the porous matrix Dm,ij (L² T^-1): hydrodynamic dispersion coefficient tensor in the porous matrix vm,i (L T^-1): linear pore water velocity in the porous matrix

qm,s (T^-1): volumetric flux of water per unit volume of porous matrix representing sources (positive) and sinks (negative)

θm (-): porosity of the porous matrix

k s

Cm_, (M L^-3): concentration of component k of the sources or sinks to the porous matrix

ξ ^{k (M L-3 T-1}): volumetric mass flux rate of solute transfer between the porous matrix and the conduit system per unit volume of component k

RXNm,k (M L^-3 T^-1): chemical reaction term of the porous matrix with respect to component k.

For the calculation of solute transport from the continuum cell to a conduit node i, the mass transfer rate is determined by multiplying the cell concentration of the pore water in the porous matrix of component k, _Cm^k_,i, at node i, with the respective exchange flow rate Γi and dividing the mass flux by the volume of the respective cell, Vm,i , i.e.,

Alternatively, if solute mass is transported from conduit node i into the matrix cell, then _Ci^k_,m is replaced by the nodal concentration, C_i^k. Fig. 3 illustrates this exchange for both cases. Such mass transfer rates are then treated as mass sink/source terms in MT3DMS since conduit nodes act similarly to other sinks/sources (e.g., wells) in MT3DMS. The difference between wells as implemented in MT3DMS and conduit nodes is that the mass removed by means of a well is not returned to the porous matrix, while the mass removed with an entry conduit node may be returned completely through the exit conduit nodes.

Amount, location and required time for these returns mainly depend on (i) the transport velocity in the different conduits, (ii) the magnitude of the exchange coefficients between the exit conduit nodes and the porous matrix, (iii) the magnitude of conduit sink terms, and (iv) the length of the different conduit tubes.

Fig. 3: Illustration of porous material cell-conduit exchange.

A one-dimensional transport equation that solely considers advection is applied to each tube, i.e.,

Cj(M L^-3): aqueous concentration of component k in conduit tube j qj (L T^-1): flux of water in tube j

t (T): time

z_j (L): distance along a respective Cartesian coordinate axis in the respective tube j.

Note there are no sink/source terms in this transport equation. The mass exchange rates from the porous matrix, from the six potentially connecting tubes and from conduit sink/source terms such as direct recharge, fixed head and fixed concentration as applied to the different conduit nodes are considered in terms of initial or boundary concentration values to the transport equation 16. The resulting concentration values are obtained by a weighted arithmetic mean of the single flow and transport components for each transport time step. Such an approach is common in mixing cell models (Bajracharya & Barry 1993).

ΓiC_i,m V_i

ΓiC_i V_i C_i,m= concentration of the cell

C_i = concentration of the conduit node i

V_i = volume of the cell

Mass exchange rate from the cell to the conduit node i

Mass exchange rate from the conduit node ito the cell

Γi = exchange flow rate ΓiC_i,m

V_i

ΓiC_i V_i C_i,m= concentration of the cell

C_i = concentration of the conduit node i

V_i = volume of the cell

Mass exchange rate from the cell to the conduit node i

Mass exchange rate from the conduit node ito the cell

Γi = exchange flow rate

Mathematically, a weighted arithmetic mean of the concentration value of component k at conduit node i can be expressed as follows:

6

subscript l: first or the maximum number of tube subsections or segments in the different tubes depending on the flow direction.

The tubes or conduits can further be divided into a (user-defined) number of sub-sections or segments as necessary to minimise numerical dispersion / improve numerical stability. Tubes within a network can be considerably longer than cell widths, lengths or thickness.

3.2.1 Transport Solvers for the Conduit System

Two numerical methods were implemented in UMT3D to solve equation 16 which are (i) the explicit finite difference (FD) method and (ii) a mass-conservative semi-Lagrangian scheme referred to as EMCNOT. The explicit finite difference (FD) method was already incorporated in the transport module for the conduit system in CAVE (Birk 2001). The EMCNOT scheme was developed by Liu et al.

(2001) for modelling advection-dominated contaminant transport problems and is an explicit mass conservative scheme without a time step limit. By means of 1-D and 2-D benchmark problems, Liu et al. (2001) demonstrated that by using EMCNOT as the solver for the transport equations, there is no stringent stability constraint on the transport time step. Moreover, it was shown that this scheme was essentially free of spurious oscillation and numerical dispersion. The EMCNOT was implemented in UMT3D with the intention to reduce numerical dispersion and increase computational efficiency of transport calculations (see Chapter 5.3).

3.2.2 Time Criterion for the Transport Solvers for the Conduit System The time criterion, i.e., the transport time step size for the FD method is the minimum residence time value within a tube segment of the conduit system multiplied by a user defined Courant number, Cr. With the EMCNOT method, a transport time step size of up to the minimum residence time value of a respective pipe, (trj)min in the conduit system multiplied by a user defined Courant number, Cr can theoretically be used for both models. The residence time of groundwater, tr in a specific pipe transported by advection under steady state flow conditions is determined by dividing the length of a pipe, Lj by the flow velocity in a respective pipe. In mathematical form, the maximum transport time step size, ∆^tmax can be expressed as:

min

Such a time criterion may vary with each flow time step since the flow rate in each tube may change with each flow time step. In Chapter 5.3 is demonstrated how time efficient the EMCNOT scheme can be to solve advective transport in the conduit system over the standard finite difference (FD) method. On top of this efficiency, the numerical dispersion can significantly be reduced with the EMCNOT method.

3.2.3 Mass Balance in the Conduit System

Mass balance in the conduit system is determined in a similar way as in the matrix. To check performance of both transport models, mass balance calculations are also carried out independently for the different conduit nodes in MT3DMS as sink/source terms. In contrast to the conduit transport model by Birk (2001) the above-described model uses a global approach to calculate mass balance and thus improves computational efficiency. Moreover, the global variable arrays and subroutines in the modified conduit transport model are fully compatible with those used in MT3DMS. All solution options in the original MT3DMS code are retained with the conduit transport model implemented as a new package named “CON”.