Getting different results from MF2005 and MF6 in the same example

[liuly12]

Hi! I try to run the unconfined-unsteady example (tutorial 2) from the modflow documentation website (https://flopy.readthedocs.io/en/3.3.2/_notebooks/tutorial02_mf.html) using both MF2005 and MF6 with the same model settings, but I got completely different results. The mf6 results seem to converge faster to the steady state than the mf2005. Here is my code below. Just wonder what would be the cause? Thank you!

`
import numpy as np
import flopy
import matplotlib.pyplot as plt
import flopy.utils.binaryfile as bf
from pathlib import Path

MF2005

Lx = 1000.0
Ly = 1000.0
ztop = 10.0
zbot = -50.0
nlay = 1
nrow = 10
ncol = 10
delr = Lx / ncol
delc = Ly / nrow
delv = (ztop - zbot) / nlay
botm = np.linspace(ztop, zbot, nlay + 1)
hk = 1.0
vka = 1.0
sy = 0.1
ss = 1.0e-4
laytyp = 1

ibound = np.ones((nlay, nrow, ncol), dtype=np.int32)
strt = 10.0 * np.ones((nlay, nrow, ncol), dtype=np.float32)

nper = 3
perlen = [1, 100, 100]
nstp = [1, 100, 100]
steady = [True, False, False]

modelname = "tutorial2_mf"
mf = flopy.modflow.Modflow(modelname, exe_name="mf2005")
dis = flopy.modflow.ModflowDis(
mf, nlay, nrow, ncol,
delr=delr, delc=delc,
top=ztop, botm=botm[1:],
nper=nper, perlen=perlen,
nstp=nstp, steady=steady,
)
bas = flopy.modflow.ModflowBas(mf, ibound=ibound, strt=strt)
lpf = flopy.modflow.ModflowLpf(
mf, hk=hk, vka=vka, sy=sy, ss=ss, laytyp=laytyp, ipakcb=53
)
pcg = flopy.modflow.ModflowPcg(mf)

GHB for SP1

stageleft = 10.0
stageright = 10.0
bound_sp1 = []
for il in range(nlay):
condleft = hk * (stageleft - zbot) * delc
condright = hk * (stageright - zbot) * delc
for ir in range(nrow):
bound_sp1.append([il, ir, 0, stageleft, condleft])
bound_sp1.append([il, ir, ncol - 1, stageright, condright])

GHB for SP2 and SP3

stageleft = 10.0
stageright = 0.0
condleft = hk * (stageleft - zbot) * delc
condright = hk * (stageright - zbot) * delc
bound_sp2 = []
for il in range(nlay):
for ir in range(nrow):
bound_sp2.append([il, ir, 0, stageleft, condleft])
bound_sp2.append([il, ir, ncol - 1, stageright, condright])

stress_period_data = {0: bound_sp1, 1: bound_sp2}
ghb = flopy.modflow.ModflowGhb(mf, stress_period_data=stress_period_data)

WEL

pumping_rate = -500.0
wel_sp1 = [[0, nrow / 2 - 1, ncol / 2 - 1, 0.0]]
wel_sp2 = [[0, nrow / 2 - 1, ncol / 2 - 1, 0.0]]
wel_sp3 = [[0, nrow / 2 - 1, ncol / 2 - 1, pumping_rate]]
stress_period_data = {0: wel_sp1, 1: wel_sp2, 2: wel_sp3}
wel = flopy.modflow.ModflowWel(mf, stress_period_data=stress_period_data)

OC

stress_period_data = {}
for kper in range(nper):
for kstp in range(nstp[kper]):
stress_period_data[(kper, kstp)] = [
"save head",
"save drawdown",
"save budget",
"print head",
"print budget",
]
oc = flopy.modflow.ModflowOc(
mf, stress_period_data=stress_period_data, compact=True
)

mf.write_input()
success, mfoutput = mf.run_model(silent=True, pause=False)
if not success:
raise Exception("MODFLOW did not terminate normally.")

headobj_mf2005 = bf.HeadFile(modelname + ".hds")
times_mf2005 = headobj_mf2005.get_times()
head_101_mf2005 = headobj_mf2005.get_data(totim=101.0)
head_201_mf2005 = headobj_mf2005.get_data(totim=201.0)

print(f"MF2005 - Day 101 well head: {head_101_mf2005[0, nrow//2-1, ncol//2-1]:.4f} m")
print(f"MF2005 - Day 201 well head: {head_201_mf2005[0, nrow//2-1, ncol//2-1]:.4f} m")
print(f"MF2005 - Day 101 head range: {head_101_mf2005.min():.4f} - {head_101_mf2005.max():.4f}")
print(f"MF2005 - Day 201 head range: {head_201_mf2005.min():.4f} - {head_201_mf2005.max():.4f}")

print("\nMF2005 Day 101 - central head in every column:")
for col in range(ncol):
print(f" 列{col}: {head_101_mf2005[0, nrow//2-1, col]:.4f} m")

ws = Path("./mf6_ghb_comparison")
ws.mkdir(exist_ok=True)
sim_name = "tutorial2_mf6"

sim = flopy.mf6.MFSimulation(
sim_name=sim_name,
sim_ws=str(ws),
exe_name="mf6",
version="mf6",
)

tdis = flopy.mf6.ModflowTdis(
sim,
time_units="DAYS",
nper=3,
perioddata=[
(1.0, 1, 1.0),
(100.0, 100, 1.0),
(100.0, 100, 1.0),
],
)

ims = flopy.mf6.ModflowIms(
sim,
print_option="SUMMARY",
complexity="simple",
linear_acceleration="bicgstab",
)

gwf = flopy.mf6.ModflowGwf(
sim,
modelname="tutorial2_mf6",
save_flows=True,
)

DIS

dis = flopy.mf6.ModflowGwfdis(
gwf,
nlay=1, nrow=10, ncol=10,
delr=100.0, delc=100.0, # Lx/ncol = 1000/10 = 100
top=10.0,
botm=[-50.0],
)

npf = flopy.mf6.ModflowGwfnpf(
gwf,
icelltype=1,
k=1.0,
k33=1.0,
save_specific_discharge=True,
save_saturation=True,
)

STO

sto = flopy.mf6.ModflowGwfsto(
gwf,
sy=0.1,
ss=1.0e-4,
steady_state={0: True, 1: False, 2: False},
transient={0: False, 1: True, 2: True},
)

ic = flopy.mf6.ModflowGwfic(gwf, strt=10.0)

SP1: both sides 10m

ghb_sp1 = []
stageleft = 10.0
stageright = 10.0
condleft = 1.0 * (10.0 - (-50.0)) * 100.0 # hk * (stage-zbot) * delc = 160100 = 6000
condright = 1.0 * (10.0 - (-50.0)) * 100.0
for ir in range(10):
ghb_sp1.append([(0, ir, 0), stageleft, condleft])
ghb_sp1.append([(0, ir, 9), stageright, condright])

SP2-3: left 10m，right0m

ghb_sp2 = []
stageleft = 10.0
stageright = 0.0
condleft = 1.0 * (10.0 - (-50.0)) * 100.0 # 6000
condright = 1.0 * (0.0 - (-50.0)) * 100.0 # 150100 = 5000
for ir in range(10):
ghb_sp2.append([(0, ir, 0), stageleft, condleft])
ghb_sp2.append([(0, ir, 9), stageright, condright])

print(f"\nGHB cond:")
print(f" SP1: 左cond={condleft}, 右cond={condright}")
print(f" SP2-3: 左cond={condleft}, 右cond={condright}")

ghb = flopy.mf6.ModflowGwfghb(
gwf,
stress_period_data={
0: ghb_sp1,
1: ghb_sp2,
2: ghb_sp2
},
)

WEL

wel = flopy.mf6.ModflowGwfwel(
gwf,
stress_period_data={
0: [[(0, 4, 4), 0.0]],
1: [[(0, 4, 4), 0.0]],
2: [[(0, 4, 4), -500.0]],
}
)

OC

oc = flopy.mf6.ModflowGwfoc(
gwf,
budget_filerecord="tutorial2_mf6.cbc",
head_filerecord="tutorial2_mf6.hds",
saverecord=[("HEAD", "ALL"), ("BUDGET", "ALL")],
)

run MF6

sim.write_simulation()
success, buff = sim.run_simulation()
if not success:
raise Exception("MODFLOW 6 did not terminate normally.")

headobj_mf6 = bf.HeadFile(str(ws / "tutorial2_mf6.hds"))
times_mf6 = headobj_mf6.get_times()
head_101_mf6 = headobj_mf6.get_data(totim=101.0)
head_201_mf6 = headobj_mf6.get_data(totim=201.0)

print(f"\nMF6 - Day 101 well head: {head_101_mf6[0, 4, 4]:.4f} m")
print(f"MF6 - Day 201 well head: {head_201_mf6[0, 4, 4]:.4f} m")
print(f"MF6 - Day 101 head range: {head_101_mf6.min():.4f} - {head_101_mf6.max():.4f}")
print(f"MF6 - Day 201 head range: {head_201_mf6.min():.4f} - {head_201_mf6.max():.4f}")

print("\nMF6 Day 101 - central head in every column:")
for col in range(10):
print(f" Column{col}: {head_101_mf6[0, 4, col]:.4f} m")

==================== Comparison ====================

print("\n" + "=" * 50)
print("Comparison")
print("=" * 50)
print(f"{'column':<5} {'MF2005':<12} {'MF6':<12} {'diff':<12}")
print("-" * 45)
for col in range(10):
diff = head_101_mf6[0, 4, col] - head_101_mf2005[0, 4, col]
print(f"{col:<5} {head_101_mf2005[0, 4, col]:<12.4f} {head_101_mf6[0, 4, col]:<12.4f} {diff:<12.4f}")

visual comparison

fig, axes = plt.subplots(2, 3, figsize=(15, 10))
levels = np.linspace(0, 10, 11)

MF2005 results

for i, (time, title) in enumerate([(1.0, "MF2005 Day 1"), (101.0, "MF2005 Day 101"), (201.0, "MF2005 Day 201")]):
ax = axes[0, i]
head = headobj_mf2005.get_data(totim=time)
cs = ax.contour(head[0], levels=levels, colors='blue')
ax.clabel(cs, inline=True, fontsize=8, fmt='%.1f')
ax.contourf(head[0], levels=np.linspace(0, 10, 21), cmap='Blues', alpha=0.3)
ax.plot(4, 4, 'ro', markersize=8)
ax.set_title(f"{title}\nWell: {head[0,4,4]:.2f}m")
ax.set_aspect('equal')
ax.grid(True, alpha=0.3)

MF results

for i, (time, title) in enumerate([(1.0, "MF6 Day 1"), (101.0, "MF6 Day 101"), (201.0, "MF6 Day 201")]):
ax = axes[1, i]
head = headobj_mf6.get_data(totim=time)
cs = ax.contour(head[0], levels=levels, colors='red')
ax.clabel(cs, inline=True, fontsize=8, fmt='%.1f')
ax.contourf(head[0], levels=np.linspace(0, 10, 21), cmap='Reds', alpha=0.3)
ax.plot(4, 4, 'ro', markersize=8)
ax.set_title(f"{title}\nWell: {head[0,4,4]:.2f}m")
ax.set_aspect('equal')
ax.grid(True, alpha=0.3)

plt.suptitle("MODFLOW-2005 vs MODFLOW 6 Comparison", fontsize=16)
plt.tight_layout()
plt.show()

timeseries compariosn

fig, ax = plt.subplots(figsize=(10, 6))
ts_mf2005 = headobj_mf2005.get_ts((0, 4, 4))
ts_mf6 = headobj_mf6.get_ts((0, 4, 4))
ax.plot(ts_mf2005[:, 0], ts_mf2005[:, 1], 'b-', linewidth=2, label='MF2005')
ax.plot(ts_mf6[:, 0], ts_mf6[:, 1], 'r--', linewidth=2, label='MF6')
ax.set_xlabel('Time (days)')
ax.set_ylabel('Head (m)')
ax.set_title('Head at Well Location (4,4)')
ax.legend()
ax.grid(True)
plt.show()
`

[liuly12]

These are the visual outputs from my code, showing the differences.

[christianlangevin]

I converted this to a discussion as it is unlikely an "issue" in either of the codes. There are a few minor differences between the two programs, particularly in terms of default behavior for the flow packages (lpf/npf/sto). Descriptions of the NPF keyword options in the input/output guide might be one place to start.

[liuly12]

Hi, Christian! Thank you for your reply. I re-scanned the files you mentioned, and found iconvert = 0 in the sto file, meaning that my mf6 setting is for confined simulation, which causes the differences. I just realised that the icelltype in npf and iconvert in sto are separately set for mf6 and should be set together at both sides to achieve the unconfined simulation. When I make the following revisions, the results are consistent from both sides. I post the results here.

# STO
sto = flopy.mf6.ModflowGwfsto(
    gwf,
    sy=0.1,
    ss=1.0e-4,
    steady_state={0: True, 1: False, 2: False},
    transient={0: False, 1: True, 2: True},
    iconvert=1,
)