fix: update SFR diversion flow to zero when its upstream downflow is …

…zero (#1467)

autotest/test_gwf_sfr_diversion.py

@@ -0,0 +1,167 @@
+import os
+import flopy
+import numpy as np
+import pytest
+from framework import TestFramework
+from simulation import TestSimulation
+
+ex = ["sfr_div"]
+inflows   = np.array([10,  0, 10,  0, 10])
+diversion = np.array([.5, .5, .5, .5, .0])
+
+def build_model(idx, ws, ):
+
+    # static model data
+    # temporal discretization
+    nper = len(inflows)
+    tdis_rc = [(1.0, 1, 1.0),] * nper
+
+    # spatial discretization data
+    nlay, nrow, ncol = 1, 1, 1
+    delr, delc = 100.0, 100.0
+    top = 0.0
+    botm = -10.0
+    strt = 0.0
+
+    # build MODFLOW 6 files
+    name = ex[idx]
+    sim = flopy.mf6.MFSimulation(
+        sim_name=name,
+        version="mf6",
+        exe_name="mf6",
+        sim_ws=ws,
+    )
+    sim.simulation_data.verify_data = False
+
+    # 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="ALL",
+    )
+
+    # 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=top,
+        botm=botm,
+    )
+
+    # initial conditions
+    ic = flopy.mf6.ModflowGwfic(gwf, strt=strt)
+
+    # node property flow
+    npf = flopy.mf6.ModflowGwfnpf(gwf, icelltype=0, )
+
+    # output control
+    oc = flopy.mf6.ModflowGwfoc(gwf, 
+                                budget_filerecord=name+'.cbb', 
+                                saverecord=[['BUDGET', 'ALL']])
+
+    # sfr file
+    cellid = (0,0,0)
+    nreaches = 3
+    rlen = 10.0
+    rwid = 10.0
+    roughness = 0.001
+    rbth = 1.0
+    rhk = 0.0
+    slope = 0.001
+
+    sfrrch_prop = [cellid, rlen, rwid, slope, top, rbth, rhk, roughness]
+    packagedata = [
+        [0,] + sfrrch_prop + [2, 1.0, 1],
+        [1,] + sfrrch_prop + [1, 0.0, 0],
+        [2,] + sfrrch_prop + [1, 1.0, 0],
+    ]
+    connectiondata = [
+        [0, -1, -2],
+        [1, 0],
+        [2, 0],
+    ]
+    diversiondata = [
+        [0, 0, 1, 'FRACTION'],
+    ]
+    perioddata = {
+        i: [[0, "inflow", qin], 
+            [0, "diversion", 0, qdiv]]
+        for i, (qin, qdiv) in enumerate(zip(inflows, diversion))}
+
+    sfr = flopy.mf6.ModflowGwfsfr(
+        gwf,
+        print_stage=True,
+        print_flows=True,
+        print_input=True,
+        budgetcsv_filerecord=name+'.sfr.csv',
+        budget_filerecord=name+'.sfr.cbb',
+        nreaches=nreaches,
+        packagedata=packagedata,
+        connectiondata=connectiondata,
+        diversions=diversiondata,
+        perioddata=perioddata,
+        pname="sfr-1",
+    )
+
+    return sim, None
+
+
+def eval_results(sim):
+    print("evaluating results...")
+
+    # check flow for indivdual reach
+    fname = os.path.join(sim.simpath, sim.name+'.sfr.cbb')
+    with flopy.utils.CellBudgetFile(fname) as cbb:
+        outflows = cbb.get_data(text='EXT-OUTFLOW')
+
+    # check outflow for reach 2 and 3
+    assert np.allclose([r.q[1] for r in outflows], -inflows*diversion), \
+           'Incorrect outflow for diversion reach'
+    assert np.allclose([r.q[2] for r in outflows], -inflows*(1-diversion)), \
+           'Incorrect outflow for outlet reach'
+
+    # load SFR budget CSV and check overall budget
+    with open(fname.replace('.cbb', '.csv')) as f:
+        header = f.readline().strip().split(',')
+        flux = np.loadtxt(f, delimiter=',')
+
+    assert np.allclose(flux[:, header.index('EXT-OUTFLOW_IN')], 0), \
+           'External flow IN larger than zero'
+    assert np.allclose(flux[:, header.index('PERCENT_DIFFERENCE')], 0), \
+           'Large mass balance error in SFR'
+
+
+@pytest.mark.parametrize(
+    "idx, name",
+    list(enumerate(ex)),
+)
+def test_mf6model(idx, name, function_tmpdir, targets):
+    ws = str(function_tmpdir)
+    test = TestFramework()
+    test.build(build_model, idx, ws)
+    test.run(
+        TestSimulation(
+            name=name, 
+            exe_dict=targets, 
+            exfunc=eval_results, 
+            idxsim=idx
+        ),
+        ws,
+    )
+

doc/ReleaseNotes/develop.tex

@@ -43,7 +43,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 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 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 SFR diversion would not be updated if the outflow of its upstream reach is zero. If diversion was not zero in the previous stress period, it would report mass balance error in the SFR budget. This bug was fixed by always updating the diversion.
 	%	\item xxx
 	\end{itemize}
 
@@ -59,4 +60,4 @@
 	%	\item xxx
 	%	\item xxx
 	%	\item xxx
-	%\end{itemize}
+	%\end{itemize}

src/Model/GroundWaterFlow/gwf3sfr8.f90

@@ -3836,6 +3836,10 @@ subroutine sfr_update_flows(this, n, qd, qgwf)
       do i = this%ia(n) + 1, this%ia(n + 1) - 1
         if (this%idir(i) > 0) cycle
         this%qconn(i) = DZERO
+        idiv = this%idiv(i)
+        if (idiv == 0) cycle
+        jpos = this%iadiv(n) + idiv - 1
+        this%divq(jpos) = DZERO
       end do
     end if
     !