fix(sfr): fix groundwater discharge to reach issue (#1717)

* add test that evaluates groundwater discharge with reach inflow equal to zero

.vscode/build_vscode.py

@@ -1,6 +1,7 @@
 import subprocess
 import os
 import argparse
+import platform
 import shutil
 import shlex
 
@@ -11,7 +12,7 @@
 args = parser.parse_args()
 
 os.environ["FC"] = args.compiler
-builddir = f"builddir_{args.compiler}_{args.buildtype}"
+builddir = f"builddir_{platform.system()}_{args.compiler}_{args.buildtype}"
 
 arg_parallel = "-Dparallel=false"
 if os.getenv("BUILD_PARALLEL_MF6") is not None:

autotest/test_gwe_sfe_strmbedcond.py

@@ -1,7 +1,6 @@
 # Test conduction between an advanced package feature, in this case stream
 # reaches with varying channel geometries and the host GWE gw cells.
 # This test should include:
-#  - no gw-sw interaction
 #  - with gw-sw interaction
 #  - hot gw cell warming an upstream reach
 #  - thermally hot stream water warming host gw cells
@@ -623,15 +622,9 @@ def check_output(idx, test):
             + str(j)
             + "in stress period 1 does not match explicitly-calculated answer"
         )
-        assert np.isclose(wa, shared_area[0, j], atol=1e-4), msg
-
-    msg = (
-        "Wetted streambed area of all reaches should be zero in stess "
-        "period 2"
-    )
-    for val in list(sfrstg[1])[1:]:
-        assert val == 0.0, msg
 
+        assert np.isclose(wa, shared_area[0, j], atol=1e-4), msg
+
     # Sub-scenario checks
     # initialize search term
     srchStr = "SFE-1 BUDGET FOR ENTIRE MODEL AT END OF TIME STEP    1, STRESS PERIOD   1"

autotest/test_gwf_sfr_gwdischarge.py

@@ -0,0 +1,120 @@
+# Test evap in SFR reaches (no interaction with gwf)
+
+import math
+import pathlib as pl
+
+import flopy
+import numpy as np
+import pytest
+
+from framework import TestFramework
+
+cases = ["sfr-gwfout"]
+
+
+def build_models(idx, test):
+    # Base simulation and model name and workspace
+    ws = test.workspace
+    name = cases[idx]
+
+    length_units = "m"
+    time_units = "sec"
+
+    nrow = 1
+    ncol = 1
+    nlay = 1
+    delr = delc = 1.0
+
+    sim = flopy.mf6.MFSimulation(
+        sim_name=name, sim_ws=ws, exe_name="mf6", version="mf6"
+    )
+    flopy.mf6.ModflowTdis(sim, time_units=time_units)
+    flopy.mf6.ModflowIms(
+        sim,
+        inner_dvclose=1e-5,
+        inner_hclose=1e-6,
+    )
+
+    gwf = flopy.mf6.ModflowGwf(
+        sim,
+        modelname=name,
+    )
+    flopy.mf6.ModflowGwfdis(
+        gwf,
+        length_units=length_units,
+        nlay=nlay,
+        nrow=nrow,
+        ncol=ncol,
+        delr=delr,
+        delc=delc,
+        top=0.0,
+        botm=-100.0,
+    )
+    flopy.mf6.ModflowGwfnpf(
+        gwf,
+        icelltype=1,  # >0 means saturated thickness varies with computed head
+    )
+    flopy.mf6.ModflowGwfic(gwf, strt=1.0)
+    flopy.mf6.ModflowGwfghb(gwf, stress_period_data=[((0, 0, 0), 1.0, 1e6)])
+
+    # sfr data
+    # <ifno> <cellid(ncelldim)> <rlen> <rwid> <rgrd> <rtp> <rbth> <rhk> <man> <ncon> <ustrf> <ndv>
+    package_data = [
+        (0, (0, 0, 0), delr, 1.0, 1e-3, 0.0, 1.0, 1.0, 0.001, 0, 0.0, 0)
+    ]
+    connection_data = [(0)]
+
+    sfr_obs = {
+        f"{name}.sfr.obs.csv": [
+            ("gwf", "sfr", (0,)),
+            ("outflow", "ext-outflow", (0,)),
+            ("depth", "depth", (0,)),
+        ],
+        "filename": name + ".sfr.obs",
+    }
+
+    flopy.mf6.ModflowGwfsfr(
+        gwf,
+        save_flows=True,
+        print_stage=True,
+        print_flows=True,
+        print_input=True,
+        length_conversion=1.0,
+        time_conversion=1.0,
+        nreaches=1,
+        packagedata=package_data,
+        connectiondata=connection_data,
+        observations=sfr_obs,
+    )
+
+    return sim, None
+
+
+def check_output(idx, test):
+    answer = np.array(
+        [
+            1.0,
+            -0.92094535738673577,
+            -0.92094535738673577,
+            0.79053721667952215e-1,
+        ]
+    )
+    obs_pth = pl.Path(f"{test.workspace}/{cases[idx]}.sfr.obs.csv")
+    sim_data = flopy.utils.Mf6Obs(obs_pth).get_data()
+    data_names = sim_data.dtype.names
+    for idx, name in enumerate(data_names):
+        assert np.allclose(
+            sim_data[name][0], answer[idx]
+        ), f"simulated sfr {name} results do not match 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_output(idx, t),
+    )
+    test.run()

doc/ReleaseNotes/develop.tex

@@ -43,6 +43,7 @@
 	\begin{itemize}
 		\item A divide by zero error would occur in the Streamflow Routing package when reaches were deactivated during a simulation. This bug was fixed by checking if the downstream reach is inactive before calculating the flow to the downstream reach.
 		\item When using the mover transport (MVT) package with UZF and UZT, rejected infiltration was paired with the calculated concentration of the UZF object rather than the user-specified concentration assigned to the infiltration.  This bug was fixed by instead pairing the rejected infiltration transferred by the MVR and MVT packages with the user-specified concentration assigned in UZTSETTING with the keyword ``INFILTRATION.''  With this change, MODFLOW 6 simulations that use UZF with the ``SIMULATE\_GWSEEP'' option will not transfer this particular source of water with the correct concentration.  Instead, the DRN package should be used to simulate groundwater discharge to land surface.  By simulating groundwater discharge to land surface with the DRN package and rejected infiltration with the UZF package, the correct concentrations will be assigned to the water transferred by the MVR package.
+		\item The Streamflow Routing package would not calculate groundwater discharge to a reach in cases where the groundwater head is above the top of the reach and the inflow to the reach from upstream reaches, specified inflows, rainfall, and runoff is zero. This bug has been fixed by eliminating the requirement that the conductance calculated based on the initial calculated stage is greater than zero in order to solve for groundwater . Instead, .As a result, differences in groundwater and surface-water exchange and groundwater heads in existing models may occur.  
 	%	\item xxx
 	\end{itemize}
 

src/Model/GroundWaterFlow/gwf-sfr.f90

@@ -184,6 +184,7 @@ module SfrModule
     procedure, private :: sfr_solve
     procedure, private :: sfr_update_flows
     procedure, private :: sfr_calc_qgwf
+    procedure, private :: sfr_gwf_conn
     procedure, private :: sfr_calc_cond
     procedure, private :: sfr_calc_qman
     procedure, private :: sfr_calc_qd
@@ -3372,7 +3373,6 @@ subroutine sfr_solve(this, n, h, hcof, rhs, update)
     real(DP) :: tp
     real(DP) :: bt
     real(DP) :: hsfr
-    real(DP) :: cstr
     real(DP) :: qd
     real(DP) :: en1
     real(DP) :: en2
@@ -3499,15 +3499,11 @@ subroutine sfr_solve(this, n, h, hcof, rhs, update)
     qd = MAX(qsrc, DZERO)
     qgwf = DZERO
     !
-    ! -- calculate reach conductance for a unit depth of water
-    !    if equal to zero will skip iterations
-    call this%sfr_calc_cond(n, d1, cstr, hsfr, hgwf)
-    !
     ! -- set flag to skip iterations
     isolve = 1
     if (hsfr <= tp .and. hgwf <= tp) isolve = 0
     if (hgwf <= tp .and. qc < DEM30) isolve = 0
-    if (cstr < DEM30) isolve = 0
+    if (this%sfr_gwf_conn(n) == 0) isolve = 0
     if (this%iboundpak(n) < 0) isolve = 0
     !
     ! -- iterate to achieve solution
@@ -4100,6 +4096,28 @@ subroutine sfr_calc_qgwf(this, n, depth, hgwf, qgwf, gwfhcof, gwfrhs)
     return
   end subroutine sfr_calc_qgwf
 
+  !> @brief Determine if a reach is connected to a gwf cell
+    !!
+    !! Function to determine if a reach is connected to a gwf cell. If connected,
+    !! the return value is 1. Otherwise, the return value is 0.
+    !!
+  !<
+  function sfr_gwf_conn(this, n)
+    ! -- return variable
+    integer(I4B) :: sfr_gwf_conn !< flag indicating if reach is connected to a gwf cell
+    ! -- dummy variables
+    class(SfrType) :: this !< SfrType object
+    integer(I4B), intent(in) :: n !< reach number
+    ! -- local variables
+    integer(I4B) :: node
+
+    sfr_gwf_conn = 0
+    node = this%igwfnode(n)
+    if (node > 0 .and. this%hk(n) > DZERO) then
+      sfr_gwf_conn = 1
+    end if
+  end function sfr_gwf_conn
+
   !> @brief Calculate reach-aquifer conductance
     !!
     !! Method to calculate the reach-aquifer conductance for a SFR package reach.
@@ -4125,6 +4143,7 @@ subroutine sfr_calc_cond(this, n, depth, cond, hsfr, htmp)
     vscratio = DONE
     !
     ! -- calculate conductance if GWF cell is active
+    !    rch-gwf flow will not occur if reach connected to an constant head cell
     node = this%igwfnode(n)
     if (node > 0) then
       if (this%ibound(node) > 0) then