fix(npf): variablecv dewatered fix (#2500)

* update test_gwf_ats02.py
* exclude test028_sfr_rewet until next release
* update release notes

Author: jdhughes-dev

autotest/test_external_models.py

@@ -19,6 +19,8 @@
 EXCLUDE = [
     "alt_model",
     "test205_gwtbuy-henrytidal",
+    # todo reinstate after 6.7.0 release
+    "test028_sfr_rewet",
     # todo reinstate after 6.5.0 release
     "test001d_Tnewton",
     # remove tests with nwt usg conductance weighting

autotest/test_gwf_ats02.py

@@ -25,7 +25,7 @@
 def build_models(idx, test):
     perlen = [10, 10]
     nper = len(perlen)
-    nstp = [5, 5]
+    nstp = [1, 1]
     tsmult = nper * [1.0]
     delr = 1.0
     delc = 1.0
@@ -75,7 +75,7 @@ def build_models(idx, test):
     )
 
     # create iterative model solution and register the gwf model with it
-    nouter, ninner = 20, 5
+    nouter, ninner = 50, 5
     hclose, rclose, relax = 1e-6, 1e-6, 0.97
     imsgwf = flopy.mf6.ModflowIms(
         sim,
@@ -216,7 +216,7 @@ def check_output(idx, test):
     except:
         assert False, f'could not load data from "{fpth}"'
     # ensure layer 1 is dry with the DRY value
-    assert np.max(tc["OBS1"][:201]) == -1.0e30, "layer 1 should be dry for this period"
+    assert np.max(tc["OBS1"][:9]) == -1.0e30, "layer 1 should be dry for this period"
 
 
 @pytest.mark.parametrize("idx, name", enumerate(cases))

autotest/test_gwf_npf_vcv.py

@@ -0,0 +1,194 @@
+"""
+Test the dependent_variable_scale option for gwt for a one-dimensional model grid
+of square cells.
+"""
+
+import flopy
+import numpy as np
+import pytest
+from framework import TestFramework
+
+dewatered = [False, True, False, True]
+perched = [False, False, True, True]
+chd1_lst = (-0.5, 0.0)
+cases = [f"vcv{n:02d}_02" for n in range(len(dewatered))]
+
+
+def build_model(idx, ws, vhfb=False):
+    nlay, nrow, ncol = 2, 1, 2
+    nper = 2
+    perlen = [1.0] * nper
+    nstp = [1] * nper
+    tsmult = [1.0] * nper
+    delr = 1.0
+    delc = 1.0
+    top = 1.0
+    botm = [0.0, -1.0]
+
+    chd0 = 0.5
+
+    strt = chd0
+
+    icelltype = 1
+    hk = 100.0  # 50000.0 # 100.0
+    area = delc * delr
+    vk = 0.001
+    if vhfb:
+        if dewatered[idx]:
+            v0 = 1.0 / ((area * vk) / (0.5 * (chd0 - botm[0])))
+            # v1 should be 0.0
+            v1 = 0.0
+        else:
+            v0 = 1.0 / ((area * vk) / (0.5 * (chd0 - botm[0])))
+            v1 = 1.0 / ((area * vk) / (0.5 * (botm[0] - botm[1])))
+        vcont = 1.0 / (v0 + v1)
+    else:
+        pass
+
+    tdis_rc = []
+    for i in range(nper):
+        tdis_rc.append((perlen[i], nstp[i], tsmult[i]))
+
+    name = cases[idx]
+
+    sim = flopy.mf6.MFSimulation(
+        sim_name=name, version="mf6", exe_name="mf6", sim_ws=ws
+    )
+
+    tdis = flopy.mf6.ModflowTdis(sim, time_units="DAYS", nper=nper, perioddata=tdis_rc)
+
+    ims = flopy.mf6.ModflowIms(
+        sim,
+        print_option="SUMMARY",
+        complexity="simple",
+        outer_dvclose=1e-8,
+        outer_maximum=200,
+        inner_dvclose=1e-9,
+        inner_maximum=100,
+    )
+
+    gwf = flopy.mf6.ModflowGwf(
+        sim,
+        modelname=name,
+        print_input=True,
+    )
+
+    dis = flopy.mf6.ModflowGwfdis(
+        gwf,
+        length_units="feet",
+        nlay=nlay,
+        nrow=nrow,
+        ncol=ncol,
+        delr=delr,
+        delc=delc,
+        top=top,
+        botm=botm,
+    )
+
+    ic = flopy.mf6.ModflowGwfic(gwf, strt=strt)
+
+    if dewatered[idx]:
+        cvoptions = "variablecv dewatered"
+    else:
+        cvoptions = "variablecv"
+
+    npf = flopy.mf6.ModflowGwfnpf(
+        gwf,
+        print_flows=True,
+        save_flows=True,
+        cvoptions=cvoptions,
+        perched=perched[idx],
+        icelltype=icelltype,
+        k=hk,
+        k33=vk,
+    )
+
+    chd_spd = {}
+    for n in range(nper):
+        spd = [((0, 0, 0), chd0)]
+        spd += [((1, 0, 0), chd1_lst[n])]
+        chd_spd[n] = spd
+    chd = flopy.mf6.ModflowGwfchd(
+        gwf,
+        maxbound=len(spd),
+        stress_period_data=chd_spd,
+    )
+
+    obs_data = {
+        "flowja.obs.csv": [
+            ("h0", "head", (0, 0, 1)),
+            ("h1", "head", (1, 0, 1)),
+            ("vflow", "FLOW-JA-FACE", (0, 0, 1), (1, 0, 1)),
+        ]
+    }
+    obs = flopy.mf6.ModflowUtlobs(gwf, continuous=obs_data)
+
+    return sim
+
+
+def build_models(idx, test):
+    ws = test.workspace
+    sim = build_model(idx, ws, vhfb=True)
+
+    return sim, None
+
+
+def check_results(idx, test):
+    ws = test.workspace
+    sim = flopy.mf6.MFSimulation.load(sim_ws=ws, verbosity_level=0)
+    gwf = sim.get_model()
+
+    nper = sim.tdis.nper.array
+    delr = gwf.dis.delr.array[1]
+    delc = gwf.dis.delc.array[0]
+    area = delr * delc
+
+    cellids = [(0, 0, 1), (1, 0, 1)]
+    botm0 = gwf.dis.botm.array[cellids[0]]
+    botm1 = gwf.dis.botm.array[cellids[1]]
+    vk = gwf.npf.k33.array[cellids[0]]
+
+    obs = gwf.obs.output.obs().get_data()
+    h0 = obs["H0"]
+    h1 = obs["H1"]
+
+    answer = obs["VFLOW"]
+
+    is_dewatered = dewatered[idx]
+    is_perched = perched[idx]
+
+    vflow = []
+    for n in range(nper):
+        if is_dewatered and h1[n] < botm0:
+            v0 = 1.0 / ((area * vk) / (0.5 * (h0[n] - botm0)))
+            v1 = 0.0
+        else:
+            v0 = 1.0 / ((area * vk) / (0.5 * (h0[n] - botm0)))
+            v1 = 1.0 / ((area * vk) / (0.5 * (botm0 - botm1)))
+        vcont = 1.0 / (v0 + v1)
+        vflow.append(vcont * (h1[n] - h0[n]))
+
+        if is_perched and h1[n] < botm0:
+            qcorr = vcont * (h1[n] - botm0)
+            vflow[n] -= qcorr
+
+    vflow = np.array(vflow)
+
+    print(f"H0: {h0}")
+    print(f"H1: {h1}")
+    print(f"Vertical flow: {answer}")
+    print(f"Answer:        {vflow}")
+
+    assert np.allclose(vflow, answer), "simulated results not equal to the answer"
+
+
+@pytest.mark.parametrize("idx, name", enumerate(cases))
+def test_mf6model(idx, name, function_tmpdir, targets):
+    test = TestFramework(
+        name=name,
+        workspace=function_tmpdir,
+        targets=targets,
+        build=lambda t: build_models(idx, t),
+        check=lambda t: check_results(idx, t),
+    )
+    test.run()

doc/ReleaseNotes/develop.toml

@@ -36,4 +36,9 @@ description = "Added HIGHEST\\_SATURATED option to the GWT SRC package. When act
 [[items]]
 section = "features"
 subsection = "basic"
-description = "Added a coordinate reference system option CRS to discretization packages. For example, an EPSG integer code, authority string, or Open Geospatial Consortium Well-Known Text (WKT) specification. This option does not affect simulation results, but is written to the binary grid file so that postprocessors can locate the grid in space. A new binary grid file version number (2) is introduced to indicate that CRS information may be present."
+description = "Added a coordinate reference system option CRS to discretization packages. For example, an EPSG integer code, authority string, or Open Geospatial Consortium Well-Known Text (WKT) specification. This option does not affect simulation results, but is written to the binary grid file so that postprocessors can locate the grid in space. A new binary grid file version number (2) is introduced to indicate that CRS information may be present."
+
+[[items]]
+section = "fixes"
+subsection = "internal"
+description = "Fixed the DEWATERED VARIABLECV option. Previously, the saturation and vertical hydraulic conductivity of the underlying cell were included in the vertical conductance calculation when the water level in the underlying cell was below the top of the cell. This is not what is described in \\citet[][equation 4-44]{modflow6gwf}. The code has been modified to only include the saturation and vertical hydraulic conductivity of the overlying cell when calculating the vertical conductance when the water-level in the underlying cell is below the top of the cell."

src/Model/GroundWaterFlow/gwf-npf.f90

@@ -1249,8 +1249,10 @@ subroutine log_options(this, found)
     if (found%ivarcv) &
       write (this%iout, '(4x,a)') 'Vertical conductance varies with water table.'
     if (found%idewatcv) &
-      write (this%iout, '(4x,a)') 'Vertical conductance accounts for dewatered &
-                                  &portion of an underlying cell.'
+      write (this%iout, '(4x,a)') 'Vertical conductance is calculated using &
+                                  &only the saturated thickness and properties &
+                                  &of the overlying cell if the head in the &
+                                  &underlying cell is below its top.'
     if (found%ixt3d) write (this%iout, '(4x,a)') 'XT3D formulation is selected.'
     if (found%ixt3drhs) &
       write (this%iout, '(4x,a)') 'XT3D RHS formulation is selected.'

src/Model/ModelUtilities/GwfConductanceUtils.f90

@@ -204,11 +204,21 @@ function vcond(ibdn, ibdm, ictn, ictm, inewton, ivarcv, idewatcv, &
       satntmp = satn
       satmtmp = satm
       if (idewatcv == 0) then
-      if (botn > botm) then
-        satmtmp = DONE
+        if (botn > botm) then
+          satmtmp = DONE
+        else
+          satntmp = DONE
+        end if
       else
-        satntmp = DONE
-      end if
+        if (botn > botm) then
+          if (satmtmp < DONE) then
+            satmtmp = DZERO
+          end if
+        else
+          if (satntmp < DONE) then
+            satntmp = DZERO
+          end if
+        end if
       end if
       bovk1 = satntmp * (topn - botn) * DHALF / vkn
       bovk2 = satmtmp * (topm - botm) * DHALF / vkm