MOD_PROCESS_MRMT.f90

Description

mod_process_MRMT.f90 is a Fortran 90 module with a numerical methodology that permits including of general multi-rate mass transfer (MRMT) processes into existing numerical codes with relative ease and efficiency.

What can be modeled with mod_process_MRMT.f90

Simultaneous effects of sorption and diffusion in different types of immobile zones such as sorption sites, clay layers and aggregates of particles.
Late-time solute transport behavior.
Transport behavior of conservative solute in continuous time random walks.
Saturated flow in heterogeneous media.

Technical support

Orlando Silva: This email address is being protected from spambots. You need JavaScript enabled to view it.

Modeling multi-rate mass transfer processes

Numerical representation

Main subroutines implemented in mod_process_MRMT.f90

Linking mod_process_MRMT.f90 to other codes

References

Software (Download)

Modeling multi-rate mass transfer processes

One of the more thorough and generalized approaches of the mass transfer between mobile and immobile domains in groundwater is the multi-rate model (Villermaux, 1987; Brusseau et al., 1989; Valocchi, 1990; Sardin et al., 1991; Haggerty and Gorelick, 1995). The Multi-Rate Mass Transfer (MRMT) approach describes simultaneous mass transfer processes (see Figure 1), governed by multiple first-order equations and a set of mass transfer rate coefficients following a probability distribution (Haggerty and Gorelick, 1995). Adding a mathematical series as a sink-source term to the flow and solute transport equations is the common manner to include MRMT in numerical simulations (e.g., Haggerty and Gorelick, 1995; Carrera et al., 1998; Salamon et al., 2006; Silva et al., 2009).

Although several of the MRMT formulations are simple, most of them has been implemented to solve a particular numerical application or are integrated in an existing code and cannot be exported to other conventional simulators. Fortran module mod_process_MRMT.f90 accounts for multi-rate mass transfer processes by adding a mathematical series as a sink-source term to the flow and solute transport equations. The module can be quite easily embedded into any existing numerical code for flow and solute transport. This numerical implementation of multi-rate mass transfer is equivalent to other existing formulations, and also is able to describe other phenomena distinct from solute transport (Silva et al., 2009).

Figure 1. Schematic representation of multiple mass transfer processes that occur in a porous medium (taken from Haggerty and Gorelick, 1995). (a) Illustration of a porous medium with heterogeneous mass transfer processes indicated by arrows. (b) Illustration of how a multirate mass transfer model would simultaneously describe the various mass transfer processes.

Numerical representation

Accounting for MRMT can be achieved in two ways: (a) using an appropriate mesh with nodes representing the immobile zones (e.g., Neumann et al., 1982), or (b) by eliminating the unknown in the immobile region as an explicit state variable (e.g., Haggerty and Gorelick, 1995; Carrera et al., 1998). Here, we have adopted the later approach because: first, it maintains the number of unknowns unchanged and, second, it is actually simpler to implement.

Figure 2 displays a schematic representation of a hypothetical numerical mesh that includes both the mobile and immobile domains. We assume that each node m of the mobile zone is connected to all adjacent nodes of the mesh and to all the immobile blocks. Node im,j of the immobile region is only connected to node m. Geometrically, node im,j overlaps with node m. Numerically, a state variable u at node im,j can be solved explicitly as a function of um, the state variable at the mobile node m. Therefore, node im,j need not be an “uncertain” node, but can be considered as a zero-D node.

Figure 2. Hypothetical numerical discretization of the mobile and immobile domains.

Main subroutines implemented in mod_process_MRMT.f90

Subroutine	Task description
Create_	It is a constructor function called at the beginning of the program. It nullifies all the pointer parameters and sets all the other parameters to zero if they are integer and real.
Destroy_	Destructor function is called in the end as it deallocates all the pointers allocated.
Read_XML_	It reads immobile attributes from *.xml files. Checks if there is any error in opening the XML file. The subroutine read_xml_loc_ is called twice. First, without optional attribute, so it reads the parameters of the immobile zones from the XML file. Later, when it is called with the optional attribute this, the parameters are allocated with space and values.
Initialize_	Allocates the variables and takes the initial condition of the state variable in the mobile region to set an initial condition in the immobile domain. This subroutine is called after reading values from the input XML file.
ContriToMatrices_	Computes the contribution to the storage matrix D of the numerical system. This is one of the main computational subroutines and should be called whenever there is a change in time step.
ContriToSink_	Computes the contribution to the source/sink term b of numerical system. This subroutine should be called at every time step in conjunction with ContriToMatrices_ in order to solve governing equations with MRMT.
UpdateConc_	Updates state variable of immobile region. The subroutine takes as arguments the current and previous time step values of the state variable in mobile region. UpdateConc_ should be called only after calling the contributions subroutines as they require previous time step state variables.
TotSolMass_	Computes the total amount of a state variable u (e.g., total mass of solute) in all the immobile zones for a particular mobile node, the total amount and the average value of u in the immobile region.
WriteConc_	Writes the state variable (e.g., concentration) in the immobile region in a DAT file.
WriteSolMass_	Writes the values of solute mass in a separate DAT output file.
GetSolMass_	Get the values of total mass of solute and the average concentration in the immobile zones

Linking mod_process_MRMT.f90 to other codes

The structure and arguments of mod_process_MRMT.f90 has been designed such as state variables and parameters characterizing the immobile domain only are accessible from the module. This guarantees a minimal information exchange with the program units of any host code to which the user would hope to link mod_process_MRMT.f90. That feature helps to get a straightforward implementation of the present MRMT method on a general advection-dispersion transport simulator programmed in Fortran, as shown in Figure 3. As illustrated in that figure, only a small number of modifications are needed.

One starts by declaring the use of mod_process_MRMT module

use mod_process_MRMT

Also, declare as many variables of type t_immobile as required by the user's application. For instance, to create a variable ImmbReg of that type

type(t_immobile) :: ImmbReg

Immediately after declaration in the main flow and transport code, include the following sentence

call Create_(ImmbReg)

This calling nullifies the attributes of Immbreg. Initialization also requires reading input files in the corresponding section of the standard code

call Read_XML_(ImmbReg,’name.xml’)

where name.xml is the input file containing information of immobile domain. After initialization in the host code one has to initialize all variables of type t_immobile, i.e.

call Initialize_(ImmbReg,mobnodes,conc0)

where mobnodes is an integer variable with the number of nodes in the mobile domain and conc0 represents the variable containing the initial values of the state variable in the mobile region. The next action concerns with the simulation process itself. After calling the host code subroutines, modules or program units that construct the system matrices A, D and b, one has to modify D and b to include MRMT. This action is shown in Figure 3 for the case of including MRMT into the solute transport problem. The subroutine ContriToMatrices_ calculates de contribution of MRMT to numerical matrix D. Include the following calling statement

call ContriToMatrices_(ImmbReg,DT,Theta,Dcorr)

ContriToMatrices_ takes the time step DT and time integration factor Theta as input arguments, returning in Dcorr a correction to the storage matrix of the numerical system due to MRMT. Analogously, a correction is needed to the sink/source term b, which can be done including the following statement

call ContriToSink_(ImmbReg,DT,Theta,Uold,Bcorr)

where Uold is a real array storing the old values of the state variable for the mobile zone. The corresponding correction to the b term is returned in Bcorr. Note that subroutines ContriToMatrices_ and ContriToSink_ must be called at every time step. They can be called within the subroutines that form matrices D and b or after them, following the structure of the host code and user’s preferences, but before execution of the unit program that solve the numerical system. After solving a time step of the transport problem (Figure 3), the state variable u of the immobile domain must be updated. To do that, include the statement

call UpdateConc_(ImmbReg,DT,Theta,Uold,U)

where the input arguments Uold and U are real arrays to store the old and current values of the state variable in the mobile domain. If required, the total amount and the average value of u in the immobile domains can be computed calling the subroutine TotSolMass_. Note that those quantities are saved within the internal structure of the variable ImmbReg. However, they can be obtained by calling the subroutine GetSolMass_. There are also available writing subroutines to follow the evolution of the state variable u at immobile regions (WriteConc_), as well as total amount and the average value of u in the immobile domains (WriteSolMass_).

Finally, it is necessary and advisable to deallocate all the attributes of all variables of type t_immobile that were used in the simulation. Therefore, just before the end of the standard transport code, include the following statement

call Destroy_(ImmbReg)

At the present, the module mod_process_MRMT.f90 was linked to TRACONF (Carrera et al., 1989), a Fortran program for the simulation of water flow and solute transport in confined aquifers. We provide the input files for three examples through which the present formulation is applied to simulate an hypothetical radially convergent tracer test, and to solve two problems of radial flow to a pumping well that are described in the work of Haggerty and Gorelick (1995).

Figure 3. Linking mod_process_MRMT.f90 to a conventional numerical code for flow and transport.

References

Brusseau, M.L., Jessup, R.E., Rao, P.S.C., 1989. Modeling the transport of solutes influenced by multiprocess nonequilibrium. Water Resources Research 25(9), 1971-1988.

Carrera, J., Galarza, G., Medina, A., 1989. TRACONF, Programa de elementos finitos para la solución de las ecuaciones de flujo y transporte en acuíferos confinados. E.T.S.I. Caminos, Canales y Puertos, Universitat Politècnica de Catalunya.

Carrera, J., Sánchez-Vila, X., Benet, I., Medina, A., Galarza, G., Guimerà, J., 1998. On matrix diffusion: formulations, solution methods and qualitative effects. Hydrogeology Journal 6, 178-190.

Haggerty, R., Gorelick, S.M., 1995. Multiple-rate mass transfer for modeling diffusion and surface reactions in media with pore-scale heterogeneity. Water Resources Research 31(10), 2383-2400.

Neumann, S.P., Preller, C., Narasimhan, T.N., 1982. Adaptive explicit-implicit quasi-three-dimensional finite element model of flow and subsidence in multiaquifer systems. Water Resources Research 18(5), 1551-1562.

Salamon, P., Fernàndez-Garcia, D., Gómez-Hernández, J.J., 2006. Modeling mass transfer processes using random walk particle tracking. Water Resources Research 42, W11417, doi:10.1029/2006WR004927.

Sardin, M., Schweinch, D., Leij, F.J., van Genuchten, M.T., 1991. Modeling the nonequilibrium transport of linearly interacting solutes in porous media. Water Resources Research 27(9), 2287-2307.

Silva, O., Carrera, J., Kumar, S., Dentz, M., Willmann, M., Alcolea, A., 2009. A simple numerical implementation to simulate multi-rate mass transfer processes (in preparation).

Valocchi, A.J., 1990. Use of temporal moment analysis to study reactive solute transport in aggregated porous media. Geoderma 46(1-3), 233-247.

Villermaux, J., 1987. Chemical engineering approach to dynamic modeling of linear chromatography. A flexible method for representing complex phenomena from simple concepts. Journal of Chromatography 406, 11-26.