refactor(LAK): specify no lake bed using DNODATA (#1321)

- [x] deprecate use of NONE keyword to define lake-gwf conductance only a function of aquifer properties - [x] update LAKE tests that used NONE keyword to use DNODATA - [x] add new LAK test to evaluate use of DNODATA, NONE, and negative number (failure) for bedleak - [x] update mf6io - [x] update release notes - [x] update mwt and ssm tex files - a result of running `python mf6ivar.py`
MODFLOW-USGS · Aug 3, 2023 · b2ee6cf · b2ee6cf
1 parent 5ad8543
commit b2ee6cf
Show file tree

Hide file tree

Showing 14 changed files with 273 additions and 57 deletions.
diff --git a/autotest/simulation.py b/autotest/simulation.py
@@ -17,6 +17,7 @@
 from flopy.utils.compare import compare_heads
 from modflow_devtools.misc import is_in_ci
 
+DNODATA = 3.e+30
 sfmt = "{:25s} - {}"
 extdict = {
     "hds": "head",

diff --git a/autotest/test_gwf_auxvars.py b/autotest/test_gwf_auxvars.py
@@ -5,7 +5,7 @@
 import numpy as np
 import pytest
 from framework import TestFramework
-from simulation import TestSimulation
+from simulation import TestSimulation, DNODATA
 
 ex = ["aux01"]
 auxvar1 = 101.0
@@ -183,8 +183,8 @@ def build_model(idx, dir):
     ]
     # <lakeno> <iconn> <cellid(ncelldim)> <claktype> <bedleak> <belev> <telev> <connlen> <connwidth>
     connectiondata = [
-        [0, 0, (0, 1, 1), "vertical", "none", 0.0, 0.0, 0.0, 0.0],
-        [1, 0, (0, 2, 2), "vertical", "none", 0.0, 0.0, 0.0, 0.0],
+        [0, 0, (0, 1, 1), "vertical", DNODATA, 0.0, 0.0, 0.0, 0.0],
+        [1, 0, (0, 2, 2), "vertical", DNODATA, 0.0, 0.0, 0.0, 0.0],
     ]
     lak = flopy.mf6.ModflowGwflak(
         gwf,

diff --git a/autotest/test_gwf_lak_bedleak.py b/autotest/test_gwf_lak_bedleak.py
@@ -0,0 +1,193 @@
+import os
+import sys
+
+import flopy
+import numpy as np
+import pytest
+from framework import TestFramework
+from simulation import TestSimulation, DNODATA
+
+ex = ["bedleak", "bedleak_fail", "bedleak_none"]
+
+
+def build_model(idx, dir):
+    nlay, nrow, ncol = 1, 10, 10
+    nper = 1
+    perlen = [
+        1.0,
+    ]
+    nstp = [
+        1,
+    ]
+    tsmult = [
+        1.0,
+    ]
+
+    lenx = 300.0
+    delr = delc = lenx / float(nrow)
+    strt = 100.0
+
+    nouter, ninner = 100, 300
+    hclose, rclose, relax = 1e-9, 1e-3, 0.97
+
+    tdis_rc = []
+    for i in range(nper):
+        tdis_rc.append((perlen[i], nstp[i], tsmult[i]))
+
+    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 iterative model solution and register the gwf model with it
+    ims = flopy.mf6.ModflowIms(
+        sim,
+        print_option="SUMMARY",
+        outer_dvclose=hclose,
+        outer_maximum=nouter,
+        under_relaxation="DBD",
+        inner_maximum=ninner,
+        inner_dvclose=hclose,
+        rcloserecord=rclose,
+        linear_acceleration="BICGSTAB",
+        scaling_method="NONE",
+        reordering_method="NONE",
+        relaxation_factor=relax,
+    )
+
+    # create gwf model
+    gwf = flopy.mf6.ModflowGwf(sim, modelname=name)
+
+    dis = flopy.mf6.ModflowGwfdis(
+        gwf,
+        nlay=nlay,
+        nrow=nrow,
+        ncol=ncol,
+        delr=delr,
+        delc=delc,
+        top=90.0,
+        botm=0.0,
+    )
+
+    # initial conditions
+    ic = flopy.mf6.ModflowGwfic(gwf, strt=strt)
+
+    # node property flow
+    npf = flopy.mf6.ModflowGwfnpf(
+        gwf, save_flows=True, icelltype=1, k=1.0, k33=0.01
+    )
+    # storage
+    sto = flopy.mf6.ModflowGwfsto(
+        gwf,
+        save_flows=True,
+        iconvert=1,
+        ss=0.0,
+        sy=0.1,
+        steady_state={0: True},
+    )
+
+    # chd files
+    chdlist0 = []
+    chdlist0.append([(0, 0, 0), 100.0])
+    chdlist0.append([(0, nrow - 1, ncol - 1), 95.0])
+
+    chdspdict = {0: chdlist0}
+    chd = flopy.mf6.ModflowGwfchd(
+        gwf,
+        stress_period_data=chdspdict,
+        save_flows=False,
+    )
+
+    # lak package
+    if "fail" in name:
+        bedleak = -100.0
+    elif "none" in name:
+        bedleak = "none"
+    else:
+        bedleak = DNODATA
+
+    # <lakeno> <strt> <nlakeconn> [<aux(naux)>] [<boundname>]
+    packagedata = [
+        [0, 100.0, 1, "lake1"],
+        [1, 100.0, 1, "lake2"],
+    ]
+    # <lakeno> <iconn> <cellid(ncelldim)> <claktype> <bedleak> <belev> <telev> <connlen> <connwidth>
+    connectiondata = [
+        [0, 0, (0, 1, 1), "vertical", bedleak, 0.0, 0.0, 0.0, 0.0],
+        [1, 0, (0, 2, 2), "vertical", bedleak, 0.0, 0.0, 0.0, 0.0],
+    ]
+    lak = flopy.mf6.ModflowGwflak(
+        gwf,
+        boundnames=True,
+        surfdep=1.0,
+        print_input=True,
+        print_stage=True,
+        print_flows=True,
+        save_flows=True,
+        budget_filerecord=f"{name}.lak.bud",
+        nlakes=len(packagedata),
+        packagedata=packagedata,
+        connectiondata=connectiondata,
+    )
+    # lak.remove()
+
+    # output control
+    oc = flopy.mf6.ModflowGwfoc(
+        gwf,
+        budget_filerecord=f"{name}.cbc",
+        headprintrecord=[("COLUMNS", 10, "WIDTH", 15, "DIGITS", 6, "GENERAL")],
+        saverecord=[("BUDGET", "ALL")],
+        printrecord=[("HEAD", "ALL"), ("BUDGET", "ALL")],
+    )
+
+    return sim, None
+
+
+def eval_model(sim):
+    print("evaluating model...")
+
+    name = ex[sim.idxsim]
+
+    # lak budget
+    if "fail" not in name:
+        fpth = os.path.join(sim.simpath, f"{name}.lak.bud")
+        bobj = flopy.utils.CellBudgetFile(fpth, precision="double")
+        bobj.list_unique_records()
+        records = bobj.get_data(text="GWF")
+        for r in records:
+            assert np.allclose(r["q"][0], -4.79616347e-12)
+            assert np.allclose(r["q"][1], -6.19237994e-12)
+
+
+@pytest.mark.parametrize(
+    "idx, name",
+    list(enumerate(ex)),
+)
+def test_mf6model(idx, name, function_tmpdir, targets):
+    ws = str(function_tmpdir)
+    if "fail" in ws:
+        require_fail = True
+    else:
+        require_fail = False
+    test = TestFramework()
+    test.build(build_model, idx, ws)
+    test.run(
+        TestSimulation(
+            name=name,
+            exe_dict=targets,
+            exfunc=eval_model,
+            idxsim=idx,
+            require_failure=require_fail,
+        ),
+        str(function_tmpdir),
+    )
diff --git a/autotest/test_gwt_lkt01.py b/autotest/test_gwt_lkt01.py
@@ -9,7 +9,7 @@
 import numpy as np
 import pytest
 from framework import TestFramework
-from simulation import TestSimulation
+from simulation import TestSimulation, DNODATA
 
 ex = ["lkt_01"]
 
@@ -134,13 +134,13 @@ def build_model(idx, dir):
     con_data = []
     # con_data=(lakeno,iconn,(cellid),claktype,bedleak,belev,telev,connlen,connwidth )
     con_data.append(
-        (0, 0, (0, 0, 1), "HORIZONTAL", "None", 10, 10, connlen, connwidth)
+        (0, 0, (0, 0, 1), "HORIZONTAL", DNODATA, 10, 10, connlen, connwidth)
     )
     con_data.append(
-        (0, 1, (0, 0, 3), "HORIZONTAL", "None", 10, 10, connlen, connwidth)
+        (0, 1, (0, 0, 3), "HORIZONTAL", DNODATA, 10, 10, connlen, connwidth)
     )
     con_data.append(
-        (0, 2, (0, 0, 2), "VERTICAL", "None", 10, 10, connlen, connwidth)
+        (0, 2, (0, 0, 2), "VERTICAL", DNODATA, 10, 10, connlen, connwidth)
     )
     p_data = [
         (0, "STATUS", "CONSTANT"),

diff --git a/autotest/test_gwt_lkt02.py b/autotest/test_gwt_lkt02.py
@@ -7,7 +7,7 @@
 import numpy as np
 import pytest
 from framework import TestFramework
-from simulation import TestSimulation
+from simulation import TestSimulation, DNODATA
 
 ex = ["lkt_02"]
 
@@ -136,21 +136,21 @@ def build_model(idx, dir):
     # con_data=(lakeno,iconn,(cellid),claktype,bedleak,belev,telev,connlen,connwidth )
     # lake 1
     con_data.append(
-        (0, 0, (0, 0, 1), "HORIZONTAL", "None", 10, 10, connlen, connwidth)
+        (0, 0, (0, 0, 1), "HORIZONTAL", DNODATA, 10, 10, connlen, connwidth)
     )
     con_data.append(
-        (0, 1, (0, 0, 2), "VERTICAL", "None", 10, 10, connlen, connwidth)
+        (0, 1, (0, 0, 2), "VERTICAL", DNODATA, 10, 10, connlen, connwidth)
     )
     # lake 2
     con_data.append(
-        (1, 0, (0, 0, 3), "VERTICAL", "None", 10, 10, connlen, connwidth)
+        (1, 0, (0, 0, 3), "VERTICAL", DNODATA, 10, 10, connlen, connwidth)
     )
     # lake 3
     con_data.append(
-        (2, 0, (0, 0, 4), "VERTICAL", "None", 10, 10, connlen, connwidth)
+        (2, 0, (0, 0, 4), "VERTICAL", DNODATA, 10, 10, connlen, connwidth)
     )
     con_data.append(
-        (2, 1, (0, 0, 5), "HORIZONTAL", "None", 10, 10, connlen, connwidth)
+        (2, 1, (0, 0, 5), "HORIZONTAL", DNODATA, 10, 10, connlen, connwidth)
     )
 
     p_data = [

diff --git a/autotest/test_gwt_lkt03.py b/autotest/test_gwt_lkt03.py
@@ -7,7 +7,7 @@
 import numpy as np
 import pytest
 from framework import TestFramework
-from simulation import TestSimulation
+from simulation import TestSimulation, DNODATA
 
 ex = ["lkt_03"]
 
@@ -134,15 +134,15 @@ def build_model(idx, dir):
     # con_data=(lakeno,iconn,(cellid),claktype,bedleak,belev,telev,connlen,connwidth )
     # lake 1
     con_data.append(
-        (0, 0, (0, 0, 2), "VERTICAL", "None", 10, 10, connlen, connwidth)
+        (0, 0, (0, 0, 2), "VERTICAL", DNODATA, 10, 10, connlen, connwidth)
     )
     # lake 2
     con_data.append(
-        (1, 0, (0, 0, 3), "VERTICAL", "None", 10, 10, connlen, connwidth)
+        (1, 0, (0, 0, 3), "VERTICAL", DNODATA, 10, 10, connlen, connwidth)
     )
     # lake 3
     con_data.append(
-        (2, 0, (0, 0, 4), "VERTICAL", "None", 10, 10, connlen, connwidth)
+        (2, 0, (0, 0, 4), "VERTICAL", DNODATA, 10, 10, connlen, connwidth)
     )
 
     p_data = [

diff --git a/autotest/test_gwt_lkt04.py b/autotest/test_gwt_lkt04.py
@@ -10,7 +10,7 @@
 import numpy as np
 import pytest
 from framework import TestFramework
-from simulation import TestSimulation
+from simulation import TestSimulation, DNODATA
 
 ex = ["lkt_04"]
 
@@ -134,13 +134,13 @@ def build_model(idx, dir, exe):
     con_data = []
     # con_data=(lakeno,iconn,(cellid),claktype,bedleak,belev,telev,connlen,connwidth )
     con_data.append(
-        (0, 0, (0, 0, 1), "HORIZONTAL", "None", 10, 10, connlen, connwidth)
+        (0, 0, (0, 0, 1), "HORIZONTAL", DNODATA, 10, 10, connlen, connwidth)
     )
     con_data.append(
-        (0, 1, (0, 0, 3), "HORIZONTAL", "None", 10, 10, connlen, connwidth)
+        (0, 1, (0, 0, 3), "HORIZONTAL", DNODATA, 10, 10, connlen, connwidth)
     )
     con_data.append(
-        (0, 2, (0, 0, 2), "VERTICAL", "None", 10, 10, connlen, connwidth)
+        (0, 2, (0, 0, 2), "VERTICAL", DNODATA, 10, 10, connlen, connwidth)
     )
     p_data = [
         (0, "STATUS", "ACTIVE"),

diff --git a/doc/ReleaseNotes/develop.tex b/doc/ReleaseNotes/develop.tex
@@ -41,8 +41,8 @@
 
 	\underline{ADVANCED STRESS PACKAGES}
 	\begin{itemize}
-		\item Added functionality to support zero values for each grid dimension when specifying the CELLID for SFR reaches that are not connected to an underlying groundwater grid cell. For example, for a DIS grid a CELLID of 0 0 0 should be specified for unconnected reaches. Warning messages will be issued if NONE is specified for unconnected reaches. Specifying a CELLID of NONE will eventually be deprecated and  will cause MODFLOW 6 to terminate with an error.
-	%	\item xxx
+		\item Added functionality to support zero values for each grid dimension when specifying the CELLID for SFR reaches that are not connected to an underlying groundwater grid cell. For example, for a DIS grid a CELLID of 0 0 0 should be specified for unconnected reaches. Warning messages will be issued if NONE is specified for unconnected reaches. Specifying a CELLID of NONE will eventually be deprecated and will cause MODFLOW 6 to terminate with an error.
+	\item Added functionality to support specification of a DNODATA (3.0E+30) BEDLEAK value for LAK package connections to identify lake-GWF connections where conductance is solely a function of aquifer properties in the connected GWF cell and lakebed sediments are assumed to be absent. Warning messages will be issued if NONE is specified for LAK package connections. Specifying a BEDLEAK value equal to NONE will eventually be deprecated and will cause MODFLOW 6 to terminate with an error.
 	%	\item xxx
 	\end{itemize}
 

diff --git a/doc/mf6io/mf6ivar/dfn/gwf-lak.dfn b/doc/mf6io/mf6ivar/dfn/gwf-lak.dfn
@@ -470,7 +470,7 @@ tagged false
 in_record true
 reader urword
 longname bed leakance
-description character string or real value that defines the bed leakance for the lake-GWF connection. BEDLEAK must be greater than or equal to zero or specified to be NONE. If BEDLEAK is specified to be NONE, the lake-GWF connection conductance is solely a function of aquifer properties in the connected GWF cell and lakebed sediments are assumed to be absent.
+description real value or character string that defines the bed leakance for the lake-GWF connection. BEDLEAK must be greater than or equal to zero, equal to the DNODATA value (3.0E+30), or specified to be NONE. If DNODATA or NONE is specified for BEDLEAK, the lake-GWF connection conductance is solely a function of aquifer properties in the connected GWF cell and lakebed sediments are assumed to be absent. Warning messages will be issued if NONE is specified. Eventually the ability to specify NONE will be deprecated and cause MODFLOW 6 to terminate with an error.
 
 block connectiondata
 name belev

diff --git a/doc/mf6io/mf6ivar/md/mf6ivar.md b/doc/mf6io/mf6ivar/md/mf6ivar.md
@@ -658,7 +658,7 @@
 | GWF | LAK | CONNECTIONDATA | ICONN | INTEGER | integer value that defines the GWF connection number for this lake connection entry. ICONN must be greater than zero and less than or equal to NLAKECONN for lake LAKENO. |
 | GWF | LAK | CONNECTIONDATA | CELLID | INTEGER (NCELLDIM) | is the cell identifier, and depends on the type of grid that is used for the simulation.  For a structured grid that uses the DIS input file, CELLID is the layer, row, and column.   For a grid that uses the DISV input file, CELLID is the layer and CELL2D number.  If the model uses the unstructured discretization (DISU) input file, CELLID is the node number for the cell. |
 | GWF | LAK | CONNECTIONDATA | CLAKTYPE | STRING | character string that defines the lake-GWF connection type for the lake connection. Possible lake-GWF connection type strings include:  VERTICAL--character keyword to indicate the lake-GWF connection is vertical  and connection conductance calculations use the hydraulic conductivity corresponding to the $K_{33}$ tensor component defined for CELLID in the NPF package. HORIZONTAL--character keyword to indicate the lake-GWF connection is horizontal and connection conductance calculations use the hydraulic conductivity corresponding to the $K_{11}$ tensor component defined for CELLID in the NPF package. EMBEDDEDH--character keyword to indicate the lake-GWF connection is embedded in a single cell and connection conductance calculations use the hydraulic conductivity corresponding to the $K_{11}$ tensor component defined for CELLID in the NPF package. EMBEDDEDV--character keyword to indicate the lake-GWF connection is embedded in a single cell and connection conductance calculations use the hydraulic conductivity corresponding to the $K_{33}$ tensor component defined for CELLID in the NPF package. Embedded lakes can only be connected to a single cell (NLAKECONN = 1) and there must be a lake table associated with each embedded lake. |
-| GWF | LAK | CONNECTIONDATA | BEDLEAK | STRING | character string or real value that defines the bed leakance for the lake-GWF connection. BEDLEAK must be greater than or equal to zero or specified to be NONE. If BEDLEAK is specified to be NONE, the lake-GWF connection conductance is solely a function of aquifer properties in the connected GWF cell and lakebed sediments are assumed to be absent. |
+| GWF | LAK | CONNECTIONDATA | BEDLEAK | STRING | real value or character string that defines the bed leakance for the lake-GWF connection. BEDLEAK must be greater than or equal to zero, equal to the DNODATA value (3.0E+30), or specified to be NONE. If DNODATA or NONE is specified for BEDLEAK, the lake-GWF connection conductance is solely a function of aquifer properties in the connected GWF cell and lakebed sediments are assumed to be absent. Warning messages will be issued if NONE is specified. Eventually the ability to specify NONE will be deprecated and cause MODFLOW 6 to terminate with an error. |
 | GWF | LAK | CONNECTIONDATA | BELEV | DOUBLE PRECISION | real value that defines the bottom elevation for a HORIZONTAL lake-GWF connection. Any value can be specified if CLAKTYPE is VERTICAL, EMBEDDEDH, or EMBEDDEDV. If CLAKTYPE is HORIZONTAL and BELEV is not equal to TELEV, BELEV must be greater than or equal to the bottom of the GWF cell CELLID. If BELEV is equal to TELEV, BELEV is reset to the bottom of the GWF cell CELLID. |
 | GWF | LAK | CONNECTIONDATA | TELEV | DOUBLE PRECISION | real value that defines the top elevation for a HORIZONTAL lake-GWF connection. Any value can be specified if CLAKTYPE is VERTICAL, EMBEDDEDH, or EMBEDDEDV. If CLAKTYPE is HORIZONTAL and TELEV is not equal to BELEV, TELEV must be less than or equal to the top of the GWF cell CELLID. If TELEV is equal to BELEV, TELEV is reset to the top of the GWF cell CELLID. |
 | GWF | LAK | CONNECTIONDATA | CONNLEN | DOUBLE PRECISION | real value that defines the distance between the connected GWF CELLID node and the lake for a HORIZONTAL, EMBEDDEDH, or EMBEDDEDV lake-GWF connection. CONLENN must be greater than zero for a HORIZONTAL, EMBEDDEDH, or EMBEDDEDV lake-GWF connection. Any value can be specified if CLAKTYPE is VERTICAL. |

diff --git a/doc/mf6io/mf6ivar/tex/gwf-lak-desc.tex b/doc/mf6io/mf6ivar/tex/gwf-lak-desc.tex
@@ -91,7 +91,7 @@
 
 \item \texttt{claktype}---character string that defines the lake-GWF connection type for the lake connection. Possible lake-GWF connection type strings include:  VERTICAL--character keyword to indicate the lake-GWF connection is vertical  and connection conductance calculations use the hydraulic conductivity corresponding to the $K_{33}$ tensor component defined for CELLID in the NPF package. HORIZONTAL--character keyword to indicate the lake-GWF connection is horizontal and connection conductance calculations use the hydraulic conductivity corresponding to the $K_{11}$ tensor component defined for CELLID in the NPF package. EMBEDDEDH--character keyword to indicate the lake-GWF connection is embedded in a single cell and connection conductance calculations use the hydraulic conductivity corresponding to the $K_{11}$ tensor component defined for CELLID in the NPF package. EMBEDDEDV--character keyword to indicate the lake-GWF connection is embedded in a single cell and connection conductance calculations use the hydraulic conductivity corresponding to the $K_{33}$ tensor component defined for CELLID in the NPF package. Embedded lakes can only be connected to a single cell (NLAKECONN = 1) and there must be a lake table associated with each embedded lake.
 
-\item \texttt{bedleak}---character string or real value that defines the bed leakance for the lake-GWF connection. BEDLEAK must be greater than or equal to zero or specified to be NONE. If BEDLEAK is specified to be NONE, the lake-GWF connection conductance is solely a function of aquifer properties in the connected GWF cell and lakebed sediments are assumed to be absent.
+\item \texttt{bedleak}---real value or character string that defines the bed leakance for the lake-GWF connection. BEDLEAK must be greater than or equal to zero, equal to the DNODATA value (3.0E+30), or specified to be NONE. If DNODATA or NONE is specified for BEDLEAK, the lake-GWF connection conductance is solely a function of aquifer properties in the connected GWF cell and lakebed sediments are assumed to be absent. Warning messages will be issued if NONE is specified. Eventually the ability to specify NONE will be deprecated and cause MODFLOW 6 to terminate with an error.
 
 \item \texttt{belev}---real value that defines the bottom elevation for a HORIZONTAL lake-GWF connection. Any value can be specified if CLAKTYPE is VERTICAL, EMBEDDEDH, or EMBEDDEDV. If CLAKTYPE is HORIZONTAL and BELEV is not equal to TELEV, BELEV must be greater than or equal to the bottom of the GWF cell CELLID. If BELEV is equal to TELEV, BELEV is reset to the bottom of the GWF cell CELLID.