From a38723ba4949103d5fbffd784e71eae36b589728 Mon Sep 17 00:00:00 2001 From: langevin-usgs Date: Wed, 16 Nov 2022 15:25:49 -0600 Subject: [PATCH] fix(uzf): uzf not routing properly (#1082) * 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 --- autotest/test_gwf_uzf05.py | 298 ++++++++++++++++++++++ doc/ReleaseNotes/ReleaseNotes.tex | 20 +- src/Model/ModelUtilities/UzfCellGroup.f90 | 4 +- 3 files changed, 312 insertions(+), 10 deletions(-) create mode 100644 autotest/test_gwf_uzf05.py diff --git a/autotest/test_gwf_uzf05.py b/autotest/test_gwf_uzf05.py new file mode 100644 index 00000000000..7258a5a03cf --- /dev/null +++ b/autotest/test_gwf_uzf05.py @@ -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() diff --git a/doc/ReleaseNotes/ReleaseNotes.tex b/doc/ReleaseNotes/ReleaseNotes.tex index 9d8aa6dafdb..65d5840ec22 100644 --- a/doc/ReleaseNotes/ReleaseNotes.tex +++ b/doc/ReleaseNotes/ReleaseNotes.tex @@ -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*} @@ -188,18 +188,19 @@ \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} @@ -207,13 +208,13 @@ \section{Changes Introduced in this Release} \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} @@ -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 @@ -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} @@ -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}. diff --git a/src/Model/ModelUtilities/UzfCellGroup.f90 b/src/Model/ModelUtilities/UzfCellGroup.f90 index 9ee7e236f2a..5e8cc879535 100644 --- a/src/Model/ModelUtilities/UzfCellGroup.f90 +++ b/src/Model/ModelUtilities/UzfCellGroup.f90 @@ -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