MF-USG Example Problems for Matrix Diffusion Transport Package
Example MD1: Single Discrete Fracture
Panday, S., 2024; USG-Transport Version 2.4.0: Transport and Other Enhancements to MODFLOW-USG, GSI Environmental, July 2024 http://www.gsi-net.com/en/software/free-software/USG-Transport.html
Several benchmark and verification simulations have been conducted with the MDT Package modules to test accuracy and performance. The code has been tested in 1-, 2-, and 3-dimensions, against analytical solutions as well as against other numerical codes; specifically,MT3D (Zheng and Wang, 1999). The following example problems are provided to demonstrate application of the MDT Process. It is recommended that users familiarize themselves with the different simulation options, code accuracy under various conditions, and input/output structures of the MDT Process via these test problems.
[1]:
import os
import matplotlib.pyplot as plt
import numpy as np
import flopy
import flopy.utils.binaryfile as bf
from flopy.mfusg import (
MfUsg,
MfUsgBas,
MfUsgBct,
MfUsgDisU,
MfUsgLpf,
MfUsgMdt,
MfUsgOc,
MfUsgPcb,
MfUsgSms,
MfUsgWel,
)
from flopy.modflow import ModflowChd, ModflowDis
from flopy.plot import PlotCrossSection, PlotMapView
from flopy.utils import HeadUFile
from flopy.utils.gridgen import Gridgen
Matrix diffusion from a single fracture (Figure Ex 10) is used to demonstrate the MODFLOWUSG MDT package’s ability to duplicate the semi-analytical solution method detailed in Falta and Wang (2017). This combination of an extremely small fracture in a much larger unfractured surface area for diffusion represents an extreme case for matrix diffusion. In this example, tritium flows through a 100 μm fracture with a Darcy velocity of 0.1 m/day at a constant concentration for 30 years (Falta and Wang, 2017; Sudicky and Frind, 1982).
The system was modeled as a one-dimension (1-D) grid along the fracture in the direction of groundwater flow with unit thickness (perpendicular to groundwater flow) representing the thickness of the fracture. To replicate Falta and Wang (2017), 61 1-m wide cells were used parallel to the groundwater flow with source concentrations assigned to the first upgradient cell. Constant head boundaries were assigned to the first and last grid cells. Model input parameters are provided in Table 3. A fracture-matrix interfacial area (Amd) of 2 m2 was used in each model cell to account for matrix diffusion from both the top and bottom of the fracture. The model was run for 50 years using a time-step of 0.1 years. Similar to the Falta and Wang (2017) setup, longitudinal dispersion was not included in the model.
A comparison of the MODFLOW-USG MDT package with the semi-analytical and analytical solutions of Falta and Wang (2017) is shown in Figure Ex 11. While the source is on, the plume continues to expand with concentrations decreasing with increasing distance from the source (year 5 and year 25 curves in Figure Ex 11). However, once the source is removed in year 30, back diffusion of tritium from the unfractured media into the fracture is what feeds the plumes observed for years 31, 33, and 50. As stated by Falta and Wang (2017):
The semi-analytical method seems to perform quite well for this demanding case, almost perfectly matching the concentrations during the loading period. At 31 years, one year after the source is removed, the semi-analytical method is able to capture the main features of the concentration reversal near the source, and it continues to perform well as the back diffusion proceeds past 50 years.
[2]:
model_ws = "Ex7_MD"
mf = MfUsg(
version="mfusg",
structured=False,
model_ws=model_ws,
modelname="Ex7_MD",
exe_name="mfusg_gsi",
)
[3]:
# one-dimension grid (61 1-m wide cells) along the fracture in the direction of groundwater flow with unit thickness
ms = flopy.modflow.Modflow()
nrow = 1
ncol = 61
delc = 1.0
delr = 1.0
nlay = 1
top = 0.0001
botm = 0.0
dis = flopy.modflow.ModflowDis(
ms, nlay, nrow, ncol, delr=delr, delc=delc, laycbd=0, top=top, botm=botm
)
[4]:
gridgen_ws = os.path.join(model_ws, "gridgen")
if not os.path.exists(gridgen_ws):
os.mkdir(gridgen_ws)
g = Gridgen(ms.modelgrid, model_ws=gridgen_ws)
g.build()
gridprops = g.get_gridprops_disu5()
anglex = g.get_anglex()
[5]:
nper = 2
perlen = [30.0, 20.0]
disu = MfUsgDisU(mf, **gridprops, itmuni=5, lenuni=1, nper=nper, perlen=perlen)
[6]:
bas = MfUsgBas(mf, strt=10.0)
[7]:
ipakcb = 41
hk = 2190.0
vka = 2190.0
lpf = MfUsgLpf(mf, ipakcb=ipakcb, constantcv=1, novfc=1, laytyp=0, hk=hk, vka=vka)
[8]:
sms = MfUsgSms(
mf,
hclose=1.0e-3,
hiclose=1.0e-6,
mxiter=250,
iter1=600,
iprsms=1,
nonlinmeth=0,
linmeth=1,
iacl=1,
norder=0,
level=7,
north=14,
iredsys=0,
rrctol=0.0,
idroptol=1,
epsrn=1.0e-3,
)
[9]:
dtype = np.dtype(
[
("node", int),
("shead", np.float32),
("ehead", np.float32),
("SHEADFACT", np.float32),
("EHEADFACT", np.float32),
("CELLGRP", int),
("c01", np.float32),
]
)
lrcsc = {
0: [[0, 11.0, 11.0, 1.0, 1.0, -1, 1.0], [60, 10.0, 10.0, 1.0, 1.0, -1, 0.0]],
1: [[0, 11.0, 11.0, 1.0, 1.0, -1, 0.0], [60, 10.0, 10.0, 1.0, 1.0, -1, 0.0]],
}
chd = ModflowChd(mf, ipakcb=ipakcb, options=[], dtype=dtype, stress_period_data=lrcsc)
[10]:
diffnc = 0.0504576
bct = MfUsgBct(
mf,
itrnsp=3,
ipakcb=27,
itvd=0,
cinact=-999.9,
nseqitr=1,
diffnc=diffnc,
prsity=0.1,
bulkd=157.0,
anglex=anglex,
dl=0.0,
dt=0.0,
ifod=1,
fodrw=0.0561,
)
[11]:
lrcsc = {0: [0, 1, 1.0], 1: [0, 1, 0.0]}
pcb = MfUsgPcb(mf, ipakcb=27, stress_period_data=lrcsc)
[12]:
mdt = MfUsgMdt(
mf,
mdflag=4,
volfracmd=1.0,
pormd=0.01,
rhobmd=2.0,
difflenmd=1.0e-10,
tortmd=0.1,
decaymd=0.0561,
diffmd=diffnc,
)
[13]:
lrcsc = {
(0, 0): [
"DELTAT 0.1",
"TMAXAT 0.2",
"TMINAT 0.01",
"TADJAT 1.0",
"TCUTAT 1.0",
"SAVE HEAD",
"SAVE BUDGET",
"SAVE CONC",
],
(1, 0): [
"DELTAT 0.1",
"TMAXAT 0.2",
"TMINAT 0.01",
"TADJAT 1.0",
"TCUTAT 1.0",
"SAVE HEAD",
"SAVE BUDGET",
"SAVE CONC",
],
}
oc = MfUsgOc(
mf, atsa=1, npsteps=1, unitnumber=[14, 30, 31, 0, 0, 132], stress_period_data=lrcsc
)
[14]:
mf.write_input()
success, buff = mf.run_model(silent=True)
True
[15]:
concobj = HeadUFile(f"{mf.model_ws}/{mf.name}.con", text="conc")
simconc = np.squeeze(concobj.get_data())
[16]:
nnodes = disu.nodes
gridx = g.get_cellxy(nnodes)[:, 0]
[17]:
fig = plt.figure(figsize=(8, 5), dpi=150)
ax = fig.add_subplot(111)
ax.plot(gridx, simconc)
ax.set_xlabel("Distance from Source (m)")
ax.set_ylabel("Normalized concentration")
ax.set_title("MODFLOW-USG Transport Simulation Results for Discrete Fracture")
[17]:
Text(0.5, 1.0, 'MODFLOW-USG Transport Simulation Results for Discrete Fracture')
