Skip to content

Commit

Permalink
fix(uzf): uzf not routing properly (#1082)
Browse files Browse the repository at this point in the history
* uzf did not route water to water table unless there was an underlying uzf cell
* Close #1075
* add test for uzf fix
* update release notes for next release
  • Loading branch information
langevin-usgs authored Nov 16, 2022
1 parent 18deae7 commit a38723b
Show file tree
Hide file tree
Showing 3 changed files with 312 additions and 10 deletions.
298 changes: 298 additions & 0 deletions autotest/test_gwf_uzf05.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,298 @@
"""
Test uzf for case where uzf is only in top cell. There was a
bug with this in the past in which case UZF would not send
water to water table unless there was a uzf in each cell.
"""

import os

import numpy as np
import pytest

try:
import pymake
except:
msg = "Error. Pymake package is not available.\n"
msg += "Try installing using the following command:\n"
msg += " pip install https://github.com/modflowpy/pymake/zipball/master"
raise Exception(msg)

try:
import flopy
except:
msg = "Error. FloPy package is not available.\n"
msg += "Try installing using the following command:\n"
msg += " pip install flopy"
raise Exception(msg)

from framework import testing_framework
from simulation import Simulation

ex = ["gwf_uzf05a"]
exdirs = []
for s in ex:
exdirs.append(os.path.join("temp", s))
ddir = "data"
nlay, nrow, ncol = 3, 1, 1

thts = 0.30 # saturated water content
thtr = 0.05 # residual water content
thti = 0.10 # initial water content
strt = 15.0


def build_model(idx, dir):

perlen = [1.0]
nper = len(perlen)
nstp = [1]
tsmult = nper * [1.0]
delr = 1.0
delc = 1.0
delv = 30.0
top = 90.0
botm = [top - (k + 1) * delv for k in range(nlay)]
laytyp = 1
ss = 1.0e-5
sy = 0.3

# unsat props
hk = 10.0
infiltration_rate = 10.
evapotranspiration_rate = 0.0
evt_extinction_depth = 2.0
brooks_corey_epsilon = 3.5 # brooks corey exponent

tdis_rc = []
for id in range(nper):
tdis_rc.append((perlen[id], nstp[id], tsmult[id]))

name = ex[idx]

# build MODFLOW 6 files
ws = dir
sim = flopy.mf6.MFSimulation(
sim_name=name, version="mf6", exe_name="mf6", sim_ws=ws
)

# create tdis package
tdis = flopy.mf6.ModflowTdis(
sim, time_units="DAYS", nper=nper, perioddata=tdis_rc
)

# create gwf model
gwfname = name
newtonoptions = "NEWTON UNDER_RELAXATION"
gwf = flopy.mf6.ModflowGwf(
sim,
modelname=gwfname,
newtonoptions=newtonoptions,
save_flows=True,
)

# create iterative model solution and register the gwf model with it
nouter, ninner = 100, 10
hclose, rclose, relax = 1.5e-6, 1e-6, 0.97
imsgwf = flopy.mf6.ModflowIms(
sim,
print_option="SUMMARY",
outer_dvclose=hclose,
outer_maximum=nouter,
under_relaxation="DBD",
under_relaxation_theta=0.7,
inner_maximum=ninner,
inner_dvclose=hclose,
rcloserecord=rclose,
linear_acceleration="BICGSTAB",
scaling_method="NONE",
reordering_method="NONE",
relaxation_factor=relax,
filename=f"{gwfname}.ims",
)
sim.register_ims_package(imsgwf, [gwf.name])

dis = flopy.mf6.ModflowGwfdis(
gwf,
nlay=nlay,
nrow=nrow,
ncol=ncol,
delr=delr,
delc=delc,
top=top,
botm=botm,
idomain=np.ones((nlay, nrow, ncol), dtype=int),
)

# initial conditions
ic = flopy.mf6.ModflowGwfic(gwf, strt=strt)

# node property flow
npf = flopy.mf6.ModflowGwfnpf(
gwf, save_flows=False, icelltype=laytyp, k=hk
)
# storage
sto = flopy.mf6.ModflowGwfsto(
gwf,
save_flows=False,
iconvert=laytyp,
ss=ss,
sy=sy,
steady_state={0: False},
transient={0: True},
)

# ghb
ghbspdict = {
0: [[(nlay - 1, 0, 0), 15.0, hk / (0.5 * delv)]],
}
ghb = flopy.mf6.ModflowGwfghb(
gwf,
print_input=True,
print_flows=True,
stress_period_data=ghbspdict,
save_flows=False,
)

# note: for specifying uzf number, use fortran indexing!
uzf_obs = {
name
+ ".uzf.obs.csv": [
(f"wc{k + 1}", "water-content", 1, depth)
for k, depth in enumerate(np.linspace(1, 20, 15))
]
}

surfdep = 1.0e-5
uzf_pkdat = [
[
0,
(0, 0, 0),
1,
-1,
surfdep,
hk,
thtr,
thts,
thti,
brooks_corey_epsilon,
"uzf01",
],
]
uzf_spd = {
0: [
[
0,
infiltration_rate,
evapotranspiration_rate,
evt_extinction_depth,
thtr,
0.0,
0.0,
0.0,
],
]
}
uzf = flopy.mf6.ModflowGwfuzf(
gwf,
print_input=True,
print_flows=True,
save_flows=True,
boundnames=True,
simulate_et=False,
unsat_etwc=True,
simulate_gwseep=True,
ntrailwaves=15,
nwavesets=40,
nuzfcells=len(uzf_pkdat),
packagedata=uzf_pkdat,
perioddata=uzf_spd,
budget_filerecord=f"{name}.uzf.bud",
wc_filerecord=f"{name}.uzf.bin",
observations=uzf_obs,
pname="uzf-1",
filename=f"{name}1.uzf",
)

# output control
oc = flopy.mf6.ModflowGwfoc(
gwf,
budget_filerecord=f"{gwfname}.bud",
head_filerecord=f"{gwfname}.hds",
headprintrecord=[("COLUMNS", 10, "WIDTH", 15, "DIGITS", 6, "GENERAL")],
saverecord=[("HEAD", "ALL"), ("BUDGET", "ALL")],
printrecord=[("HEAD", "LAST"), ("BUDGET", "ALL")],
)

obs_lst = []
obs_lst.append(["obs1", "head", (0, 0, 0)])
obs_dict = {f"{gwfname}.obs.csv": obs_lst}
obs = flopy.mf6.ModflowUtlobs(
gwf, pname="head_obs", digits=20, continuous=obs_dict
)

return sim, None


def eval_flow(sim):
print("evaluating flow...")

name = ex[sim.idxsim]
ws = exdirs[sim.idxsim]

fname = os.path.join(ws, f"{name}.uzf.bin")
wobj = flopy.utils.HeadFile(fname, text="WATER-CONTENT")
wc = wobj.get_alldata()

fname = os.path.join(ws, f"{name}.hds")
wobj = flopy.utils.HeadFile(fname)
head = wobj.get_alldata()

bpth = os.path.join(ws, name + ".uzf.bud")
bobj = flopy.utils.CellBudgetFile(bpth, precision="double")
flowtogwf = bobj.get_data(text="GWF")[0]

print("Checking to make sure UZF cell 1 sends water to GWF cell 1")
node, node2, q, flow_area = flowtogwf[0]
assert node == 1, "uzf node should be 1"
assert node2 == 1, "GWF node should be 1"
assert np.isclose(q, -4.), "Flow from UZF to node 1 should be -4."

return


# - No need to change any code below
@pytest.mark.parametrize(
"idx, dir",
list(enumerate(exdirs)),
)
def test_mf6model(idx, dir):
# initialize testing framework
test = testing_framework()

# build the model
test.build_mf6_models(build_model, idx, dir)

# run the test model
test.run_mf6(Simulation(dir, exfunc=eval_flow, idxsim=idx))


def main():
# initialize testing framework
test = testing_framework()

# run the test model
for idx, dir in enumerate(exdirs):
test.build_mf6_models(build_model, idx, dir)
sim = Simulation(dir, exfunc=eval_flow, idxsim=idx)
test.run_mf6(sim)

return


if __name__ == "__main__":
# print message
print(f"standalone run of {os.path.basename(__file__)}")

# run main routine
main()
20 changes: 12 additions & 8 deletions doc/ReleaseNotes/ReleaseNotes.tex
Original file line number Diff line number Diff line change
Expand Up @@ -173,7 +173,7 @@ \section{Release History}
6.2.1 & February 18, 2021 \\
6.2.2 & July 30, 2021 \\
6.3.0 & March 4, 2022 \\
6.x.x & Month xx, 202x (this release) \\
6.4.0 & November 30, 2022 (this release) \\
\hline
\label{tab:releases}
\end{tabular*}
Expand All @@ -188,32 +188,33 @@ \section{Changes Introduced in this Release}

\begin{itemize}

\item Version mf6.x.x--Month xx, 202x
%\item Version mf6.x.x--Month xx, 202x
\item Version mf6.4.0--November 30, 2022

\underline{NEW FUNCTIONALITY}
\begin{itemize}
\item A new Viscosity (VSC) package for the Groundwater Flow (GWF) Model is introduced in this release. The effects of viscosity are accounted for by updates to intercell conductance, as well as the conductance between the aquifer and head-dependent boundaries, based on simulated concentrations and\/or temperatures. The viscosity package is activated by specifying ``VSC6'' as the file type in a GWF name file. Changes to the code and input may be required in response to user needs and testing.
\item A new Viscosity (VSC) package for the Groundwater Flow (GWF) Model is introduced in this release. The effects of viscosity are accounted for by updates to intercell conductance, as well as the conductance between the aquifer and head-dependent boundaries, based on simulated concentrations and\/or temperatures. The VSC Package is activated by specifying ``VSC6'' as the file type in a GWF name file. Changes to the code and input may be required in the future in response to user needs and testing. Implementation details for the VSC Package are described in the Supplemental Technical Information guide, which is included with the MODFLOW 6 distribution.
% \item xxx
% \item xxx
\end{itemize}

\underline{EXAMPLES}
\begin{itemize}
\item A new example called ex-gwt-stallman was added. This new problem uses the GWT Model as a surrogate for simulating heat flow. The example represents heat conduction in subsurface with a periodic temperature boundary condition at the surface.
\item A new example called ex-gwt-stallman was added. This new problem uses the GWT Model as a surrogate for simulating heat flow. The example represents one-dimensional heat convection and conduction in the subsurface in response to a periodic temperature boundary condition imposed at land surface. Results from the MODFLOW 6 simulation are in good agreement with an analytical solution.
% \item xxx
% \item xxx
\end{itemize}

\textbf{\underline{BUG FIXES AND OTHER CHANGES TO EXISTING FUNCTIONALITY}} \\
\underline{BASIC FUNCTIONALITY}
\begin{itemize}
\item Corrected programming error in XT3D functionality that could occur when running coupled flow models or transport models. The XT3D code would result in a memory access error when a child model with a much larger level of refinement was coupled to a coarser parent model. The XT3D code was generalized to handle this situation.
\item Corrected programming error in XT3D functionality that could affect coupled flow models or coupled transport models. The XT3D code would result in a memory access error when a child model with a much larger level of refinement was coupled to a coarser parent model. The XT3D code was generalized to handle this situation.
\item Corrected a programming error in which the final message would be written twice to the screen and twice to mfsim.lst when the simulation terminated prematurely.
\item Terminate with error if METHOD or METHODS not specified in time series input files. Prior to this change, the program would continue without an interpolated value for one or more time series records.
\item When a GWF Model and a corresponding GWT model are solved in the same simulation, the GWF Model must be solved before the corresponding GWT model. The GWF Model must also be solved by a different IMS than the GWT Model. There was not a check for this in previous versions and if these conditions were not met, the solution would often not converge or it would give erroneous results.
\item The DISV Package would not raise an error if a model cell was defined as a line. The program was modified to check for the case where the calculated cell area is equal to zero. If the calculated cell area is equal to zero, the program terminates with an error.
\item When searching for a required block in an input file, the program would not terminate with a sensible error message if the end of file was found instead of the required block. Program now indicates that the required block was not found.
\item This release contains a first step toward implementation of generic input routines to read input files. The new input routines were implemented for the DIS, DISV, and DISU Packages of the GWF and GWT Models, for the NPF Package of the GWF Model, and the DSP Package of the GWT Model. Output summaries written to the GWF and GWT Model listing files are different from summaries written using previous versions. For packages that use the new input data model, the IPRN capability of the READARRAY utility (described in mf6io.pdf) is no longer supported as a way to write input arrays to the model listing file.
\item This release contains a first step toward implementation of generic input routines to read input files. The new input routines were implemented for the DIS, DISV, and DISU Packages of the GWF and GWT Models, for the NPF Package of the GWF Model, and the DSP Package of the GWT Model. Output summaries written to the GWF and GWT Model listing files are different from summaries written using previous versions of MODFLOW 6. For packages that use the new input data model, the IPRN capability of the READARRAY utility (described in mf6io.pdf) is no longer supported as a way to write input arrays to the model listing file. The IPRN capability may not be supported in future versions as the new generic input routines are implemented for other packages.
\end{itemize}

\underline{INTERNAL FLOW PACKAGES}
Expand All @@ -225,7 +226,7 @@ \section{Changes Introduced in this Release}

\underline{STRESS PACKAGES}
\begin{itemize}
\item The Evapotranspiration (EVT) Package was modified to include a new error check if the segmented evapotranspiration capability is active. If the number of ET segments is greater than 1, then the user must specify values for PXDP (as well as PETM). For a cell, PXDP is a one-dimensional array of size NSEG - 1. Values in this array must be greater than zero and less than one. Furthermore, the values in PXDP must increase monotonically. The program now checks for these conditions and terminates with an error if these conditions are not met. The segmented ET capability can used be used for list-based EVT Package input. Provided that the PXDP conditions are met, this new error check should have no effect on simulation results.
\item The Evapotranspiration (EVT) Package was modified to include a new error check if the segmented evapotranspiration capability is active. If the number of ET segments is greater than 1, then the user must specify values for PXDP (as well as PETM). For a cell, PXDP is a one-dimensional array of size NSEG - 1. Values in this array must be greater than zero and less than one. Furthermore, the values in PXDP must increase monotonically. The program now checks for these conditions and terminates with an error if these conditions are not met. The segmented ET capability can be used for list-based EVT Package input. Provided that the PXDP conditions are met, this new error check should have no effect on simulation results.
\item The Evapotranspiration (EVT) Package would throw an index error when SURF\_RATE\_SPECIFIED was specified in the OPTIONS block and NSEG was set equal to 1. The code now supports this combination of input options.
% \item xxx
% \item xxx
Expand All @@ -251,7 +252,7 @@ \section{Changes Introduced in this Release}

\underline{EXCHANGES}
\begin{itemize}
\item The GWT-GWT Exchange did not work when the XT3D\_OFF option was specified. Fixed the program so that the XT3D dispersion terms can be shut off for a GWT-GWT Exchange. GWT-GWT Exchange keywords were renamed from ADVSCHEME, XT3D\_OFF, and XT3D\_RHS to ADV\_SCHEME, DSP\_XT3D\_OFF, and DSP\_XT3D\_OFF, respectively, to more clearly indicate how the keywords relate to the underlying processes.
\item The GWT-GWT Exchange did not work when the XT3D\_OFF option was specified. The program was fixed so that the XT3D dispersion terms can be shut off for a GWT-GWT Exchange. GWT-GWT Exchange keywords were renamed from ADVSCHEME, XT3D\_OFF, and XT3D\_RHS to ADV\_SCHEME, DSP\_XT3D\_OFF, and DSP\_XT3D\_RHS, respectively, to more clearly indicate how the keywords relate to the underlying processes.
% \item xxx
% \item xxx
\end{itemize}
Expand Down Expand Up @@ -286,6 +287,9 @@ \section{Known Issues and Incompatibilities}
\item
If a GWF-GWF Exchange is activated with the XT3D option, then the two connected GWF Models cannot have BUY Packages active.

\item
The Time-Variable Hydraulic Conductivity (TVK) Package is incompatible with the Horizontal Flow Barrier (HFB) Package if the TVK Package is used to change hydraulic properties of cells near horizontal flow barriers.

\end{enumerate}

In addition to the issues shown here, a comprehensive and up-to-date list is available under the issues tab at \url{https://github.com/MODFLOW-USGS/modflow6}.
Expand Down
4 changes: 2 additions & 2 deletions src/Model/ModelUtilities/UzfCellGroup.f90
Original file line number Diff line number Diff line change
Expand Up @@ -708,9 +708,9 @@ subroutine solve(this, thiswork, jbelow, icell, totfluxtot, ietflag, &
this%surfluxbelow(icell) = DZERO
if (this%ivertcon(icell) > 0) then
this%finf(jbelow) = DZERO
if (this%watab(icell) < this%celbot(icell)) &
this%watab(icell) = this%celbot(icell)
end if
if (this%watab(icell) < this%celbot(icell)) &
this%watab(icell) = this%celbot(icell)
!
! -- initialize derivative variables
deriv1 = DZERO
Expand Down

0 comments on commit a38723b

Please sign in to comment.