MODFLOW-USGS
Trivial model does not converge or crashes (floating overflow) with buoyancy
Describe the bug A colleague came across a problem where MODFLOW 6 did not converge, whereas our Deltares SEAWAT version iMOD-WQ did (and very quickly). Our suspicion is that the problem lies with river boundaries where the stage equals the bottom, in combination with the buoyancy correction (since a similar issue required changes in iMOD-WQ). In this case, the River basically acts as the Drainage package. Using a drainage package makes the convergence issues disappear. However, it seems like a river should be allowed to fall dry during a simulation; the input validation doesn't disallow it either and the physical meaning of stage equals bottom is clear.
I've been trying to make a minimal working example:
- 3 rows
- 3 columns
- 1 river boundary in the center
- constant recharge which is forced to the river boundary
- my experiment consists of replacing the river boundary by a functionally equivalent drainage boundary
- model defines no storage, so it's steady-state, single stress period
And I seem to run into a number of connected issues, which I think are primarily due to buoyancy since running everything with 0 concentrations works. I also get float overflow crashes, and more so with BICGSTAB instead of CG. (In all cases I use the hhformulation_rhs=True.)
There are some constants in the script for easy variation.
- drainage, no transport no buoyancy, CG: converges.
- river, no transport no buoyancy, CG: converges
- river, no transport no buoyancy, BICGSTAB: converges
- drainage, buoyancy, CG: converges
- river, buoyancy all fresh, CG: converges
- river, buoyancy all fresh, BICGSTAB: converges
- river, buoyancy all salt (initial concentration=16.0), BICGSTAB: floating overflow crash
- river, buoyancy all salt (16.0), CG: convergence failure
- drainage, buoyancy all salt (16.0), CG: converges
- drainage, buoyancy all salt (16.0), BICGSTAB: floating overflow crash
- river, buoyancy all salt (16.0), stage 0.5 m above bottom, CG: floating overflow crash
- river, buoyancy all salt (16.0), stage 0.5 m above bottom, BICGSTAB: floating overflow crash
- river, buoyancy all fresh (0.0), stage 0.5 m above bottom, BICGSTAB: floating overflow crash
- river, buoyancy all fresh (0.0, stage 0.5 m above bottom, CG: converges
To summarize: it's the combination of the river with nonzero concentration buoyancy that appears to be problematic. I reckon there's a bug in the river buoyancy formulation. I also find it suspicous that CG converges in some cases whereas BICGSTAB crashes, which points to a bug in the BICSTAB implementation as well.
To reproduce
Run the following script. I have defined some upper cased constants to allow easy variation for the scenarios above.
# %%
import flopy
# %%
USE_RIVER = False
STAGE_ABOVE_BOTTOM = 0.0
USE_BUOYANCY = False
FLOW_LINEAR_ACCELERATION="CG"
INITIAL_HEAD = 0.0
INITIAL_CONCENTRATION = 0.0 # Initial concentration (unitless)
# %%
# Model units
length_units = "centimeters"
time_units = "seconds"
# Model parameters
nper = 1 # Number of periods
nstp = 1 # Number of time steps
perlen = 1.0 # Simulation time length ($d$)
nlay = 1 # Number of layers
nrow = 3 # Number of rows
ncol = 3 # Number of columns
delr = 10.0 # Column width ($m$)
delc = 10.0 # Row width ($m$)
delv = 1.0 # Layer thickness
top = 1.0 # Top of the model ($m$)
hydraulic_conductivity = 1.0 # Hydraulic conductivity ($m d^{-1}$)
porosity = 0.35 # porosity (unitless)
diffusion_coefficient = 0.57024 # diffusion coefficient ($m^2/d$)
botm = [top - k * delv for k in range(1, nlay + 1)]
nouter, ninner = 100, 300
hclose, rclose, relax = 1e-10, 1e-6, 0.97
name = "flow"
sim = flopy.mf6.MFSimulation(sim_name="test", sim_ws="test", exe_name="mf6")
tdis_ds = ((perlen, nstp, 1.0),)
flopy.mf6.ModflowTdis(sim, nper=nper, perioddata=tdis_ds, time_units=time_units)
gwf = flopy.mf6.ModflowGwf(sim, modelname=name, save_flows=True)
ims = flopy.mf6.ModflowIms(
sim,
print_option="ALL",
outer_dvclose=hclose,
outer_maximum=nouter,
under_relaxation="NONE",
inner_maximum=ninner,
inner_dvclose=hclose,
rcloserecord=rclose,
linear_acceleration=FLOW_LINEAR_ACCELERATION,
scaling_method="NONE",
reordering_method="NONE",
relaxation_factor=relax,
filename=f"{gwf.name}.ims",
)
sim.register_ims_package(ims, [gwf.name])
flopy.mf6.ModflowGwfdis(
gwf,
length_units=length_units,
nlay=nlay,
nrow=nrow,
ncol=ncol,
delr=delr,
delc=delc,
top=top,
botm=botm,
)
flopy.mf6.ModflowGwfnpf(
gwf,
save_specific_discharge=True,
icelltype=0,
k=hydraulic_conductivity,
)
flopy.mf6.ModflowGwfic(gwf, strt=INITIAL_HEAD)
stage = top + STAGE_ABOVE_BOTTOM
bot = top
cond = 10.0
rivconc = 0.0
if USE_RIVER:
flopy.mf6.ModflowGwfriv(
gwf,
stress_period_data=[[(0, 1, 1), stage, cond, bot, rivconc]],
pname="RIV-1",
auxiliary="CONCENTRATION",
)
else:
flopy.mf6.ModflowGwfdrn(
gwf,
stress_period_data=[[(0, 1, 1), stage, cond, rivconc]],
pname="RIV-1",
auxiliary="CONCENTRATION",
)
head_filerecord = f"{name}.hds"
budget_filerecord = f"{name}.bud"
flopy.mf6.ModflowGwfoc(
gwf,
head_filerecord=head_filerecord,
budget_filerecord=budget_filerecord,
saverecord=[("HEAD", "ALL"), ("BUDGET", "ALL")],
)
flopy.mf6.ModflowGwfrcha(
gwf,
recharge=0.001,
)
if USE_BUOYANCY:
gwt = flopy.mf6.ModflowGwt(sim, modelname="trans")
imsgwt = flopy.mf6.ModflowIms(
sim,
print_option="ALL",
outer_dvclose=hclose,
outer_maximum=nouter,
under_relaxation="NONE",
inner_maximum=ninner,
inner_dvclose=hclose,
rcloserecord=rclose,
linear_acceleration="BICGSTAB",
scaling_method="NONE",
reordering_method="NONE",
relaxation_factor=relax,
filename=f"{gwt.name}.ims",
)
sim.register_ims_package(imsgwt, [gwt.name])
pd = [(0, 1.3, 0.0, "trans", "concentration")]
flopy.mf6.ModflowGwfbuy(gwf, packagedata=pd, hhformulation_rhs=True)
flopy.mf6.ModflowGwtmst(gwt, porosity=porosity)
flopy.mf6.ModflowGwtic(gwt, strt=INITIAL_CONCENTRATION)
flopy.mf6.ModflowGwtadv(gwt, scheme="UPSTREAM")
flopy.mf6.ModflowGwtdsp(gwt, xt3d_off=True, diffc=diffusion_coefficient)
sourcerecarray = [
("RIV-1", "AUX", "CONCENTRATION"),
]
flopy.mf6.ModflowGwtssm(gwt, sources=sourcerecarray)
flopy.mf6.ModflowGwtoc(
gwt,
budget_filerecord=f"{gwt.name}.cbc",
concentration_filerecord=f"{gwt.name}.ucn",
concentrationprintrecord=[("COLUMNS", 10, "WIDTH", 15, "DIGITS", 6, "GENERAL")],
saverecord=[("CONCENTRATION", "ALL")],
printrecord=[("CONCENTRATION", "LAST"), ("BUDGET", "LAST")],
)
flopy.mf6.ModflowGwfgwt(
sim, exgtype="GWF6-GWT6", exgmnamea=gwf.name, exgmnameb=gwt.name
)
flopy.mf6.ModflowGwtdis(
gwt,
length_units=length_units,
nlay=nlay,
nrow=nrow,
ncol=ncol,
delr=delr,
delc=delc,
top=top,
botm=botm,
)
sim.write_simulation()
# %%
sim.run_simulation()
# %%
Expected behavior
- Both CG and BICGSTAB should be able to solve this (the linear problem is simple and small!)
- A River should be capable of acting identical to a Drainage boundary when stage equals bottom; it should not give convergence issues.
Screenshots This is the typical head result we'd expect:
Environment mf6.6.0_win64
Just had a thought: obviously if the initial stage is higher than the initial head estimate (for a steady-state computation), a linearized formulation is ill-posed since there's only fixed fluxes (recharge and river infiltration over stage - bottom) coming in and nothing coming out.
I should add a constant head cell or better a very low conductance GHB boundary, and run the scenarios again.
@Huite We you able to try adding a head-based boundary condition as you had proposed, and if so, did it resolve the problem?
