feat(sfr): kinematic wave approximation option

autotest/test_gwf_sfr_kinematic01.py

@@ -0,0 +1,185 @@
+"""
+This test sfr kinematic wave approximation using example
+in Ponce (1989) Example 9-3 part 1. In Ponce (1989) the
+Courant number was specified, was invariant, and was used
+directly in the kinematic wave equation. As a result, the
+of the sfr result is approximate since the sfr kinematic
+wave approximation is a numerical method where the Courant
+number is a function of the solution.
+
+The test has 1 sfr reach that is not connected to the
+groundwater model.
+
+The test evaluates EXT-OUTFLOW for reaches 1 against
+the known numerical solution for storage weight values
+set to 0.5 and 1.0.
+
+Ponce, V. M. (1989). Engineering Hydrology, Principles and Practices.
+"""
+
+import flopy
+import numpy as np
+import pytest
+from framework import TestFramework
+
+paktest = "sfr"
+cases = ("sfr-kwe01", "sfr-kwe02")
+storage_weights = (0.5, 1.0)
+
+
+def build_models(idx, test):
+    name = cases[idx]
+
+    dt = 60.0 * 60.0
+    flows = np.array(
+        [
+            0.0,
+            30.0,
+            60.0,
+            90.0,
+            120.0,
+            150.0,
+            120.0,
+            90.0,
+            60.0,
+            30.0,
+            0.0,
+            0.0,
+            0.0,
+            0.0,
+        ]
+    )
+    times = np.arange(0.0, (flows.shape[0]) * dt, dt)
+
+    # static model data
+    # temporal discretization
+    nper = times.shape[0] - 1
+    nstp = 1
+    perioddata = [(dt, nstp, 1.0) for idx in range(nper)]
+
+    # spatial discretization data
+    nlay, nrow, ncol = 1, 1, 6
+    delr, delc = 100.0, 100.0
+    dx = 7200.0
+    top = 100.0
+    botm = 0.0
+
+    sim = flopy.mf6.MFSimulation(
+        sim_name=name,
+        sim_ws=test.workspace,
+        version="mf6",
+        exe_name="mf6",
+    )
+    tdis = flopy.mf6.ModflowTdis(
+        sim, time_units="seconds", nper=nper, perioddata=perioddata
+    )
+    ims = flopy.mf6.ModflowIms(sim)
+
+    gwf = flopy.mf6.ModflowGwf(sim, modelname=name)
+    dis = flopy.mf6.ModflowGwfdis(
+        gwf,
+        length_units="meters",
+        nlay=nlay,
+        nrow=nrow,
+        ncol=ncol,
+        delr=dx,
+        delc=dx,
+        top=top,
+        botm=botm,
+    )
+    npf = flopy.mf6.ModflowGwfnpf(gwf, icelltype=1)
+    ic = flopy.mf6.ModflowGwfic(gwf, strt=top)
+    sto = flopy.mf6.ModflowGwfsto(
+        gwf, iconvert=1, ss=1e-6, sy=0.2, transient={0: True}
+    )
+    ghb = flopy.mf6.ModflowGwfghb(
+        gwf, stress_period_data=[(0, 0, 0, top, 5.0)]
+    )
+    oc = flopy.mf6.ModflowGwfoc(gwf, printrecord=[("budget", "all")])
+
+    # sfr file
+    nreaches = 1
+    slope = 1.0 / dx
+    roughness = 0.03574737676661647
+
+    # <ifno> <cellid(ncelldim)> <rlen> <rwid> <rgrd> <rtp> <rbth> <rhk> <man> <ncon> <ustrf> <ndv> [<aux(naux)>] [<boundname>]
+    pak_data = [
+        (
+            ifno,
+            -1,
+            -1,
+            -1,
+            dx,
+            10.0,
+            slope,
+            top,
+            1.0,
+            0.0,
+            roughness,
+            0,
+            0.0,
+            0,
+        )
+        for ifno in range(nreaches)
+    ]
+    sfr_spd = {idx: [(0, "inflow", q)] for idx, q in enumerate(flows[1:])}
+    sfr = flopy.mf6.ModflowGwfsfr(
+        gwf,
+        print_flows=True,
+        storage=True,
+        storage_weight=storage_weights[idx],
+        nreaches=nreaches,
+        packagedata=pak_data,
+        connectiondata=[(0,)],
+        perioddata=sfr_spd,
+        pname="sfr-1",
+    )
+    fname = f"{name}.sfr.obs"
+    sfr_obs = {
+        f"{fname}.csv": [
+            ("inflow", "ext-inflow", (0,)),
+            ("outflow", "ext-outflow", (0,)),
+        ]
+    }
+    sfr.obs.initialize(filename=fname, print_input=True, continuous=sfr_obs)
+
+    return sim
+
+
+def check_output(idx, test):
+    print("Checking sfr external outflow")
+    name = cases[idx]
+    obs_values = flopy.utils.Mf6Obs(
+        test.workspace / f"{name}.sfr.obs.csv"
+    ).get_data()
+    # fmt: off
+    test_values = {storage_weights[0]: np.array(
+        [ 
+         0.        ,  -10.68920803,  -58.58988171,  -92.59647491, -125.23616408, 
+         -144.76383592, -117.40352509,  -89.57181956, -64.81693279,  -53.46056541,  
+         -11.22910318,   -4.95517783, -2.72890377
+         ]
+        ),
+        storage_weights[1]: np.array(
+        [ 
+         0.        ,  -15.3918031 ,  -58.58988171,  -92.59647491, -125.23616408, 
+         -144.76383592, -117.40352509,  -88.84252125, -62.90365264,  -47.01932528,  
+         -19.87588927,   -9.96585666, -5.63155787
+         ]
+        ),
+    }
+    assert np.allclose(
+            obs_values["OUTFLOW"], test_values[storage_weights[idx]]
+        ), f"failed comparison for '{name}' observation"
+
+
+@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

@@ -1,13 +1,14 @@
-% Use this template for starting initializing the release notes
-% after a release has just been made.
+% Use this template for starting initializing the release notes after a release
+% has just been made.
 
 	\item \currentmodflowversion
 
 	\underline{NEW FUNCTIONALITY}
 	\begin{itemize}
 		\item A new adaptive time stepping (ATS) capability was added to the Advection (ADV) Package of the Groundwater Transport (GWT) Model.  A new input option, called ATS\_PERCEL, specifies the fractional cell distance that a particle of water can travel within one time step.  When ATS\_PERCEL is specified by the user, and the ATS utility is activated in the TDIS Package, the ADV Package will calculate the largest time step that will meet this fractional cell distance constraint, and will submit this time step to the ATS utility.  This option may improve time stepping for solute transport models and for variable-density flow and transport models by allowing step lengths to be calculated as a function of the flow system rather than being specified as input by the user.
 		\item Added the capability to write sorbate concentrations to a binary output file.  A new SORBATE option is now available in the Mass Storage and Transfer (MST) Package of the GWT Model to provide the name of the binary output file for sorbate concentrations. Sorbate concentrations will be written to the binary output file whenever concentrations for the GWT model are saved, as determined by the GWT Output Control option.
-	%	\item xxx
+		\item Add kinematic-wave routing option for the Streamflow Routing (SFR) Package. Prior to this change, the SFR Package could only simulate unidirectional, steady, uniform flow. With kinematic-wave routing, unidirectional waves can now propagate through the SFR network by explicitly including a storage term in the reach continuity equation. The kinematic-wave routing option is based on the ``TRANSROUTE'' option available in the SFR Package in MODFLOW-NWT. The kinematic-wave routing option is enabled by specifying ``STORAGE'' in the SFR Package OPTIONS block.
+		%	\item xxx
 	\end{itemize}
 
 	\underline{EXAMPLES}

doc/mf6io/gwf/sfr.tex

@@ -28,6 +28,10 @@ \subsubsection{Structure of Blocks}
 \noindent \textit{IF ndv IS GREATER THAN ZERO FOR ANY REACH}
 \lstinputlisting[style=blockdefinition]{./mf6ivar/tex/gwf-sfr-diversions.dat}
 
+\vspace{5mm}
+\noindent \textit{INITIALSTAGES BLOCK IS OPTIONAL}
+\lstinputlisting[style=blockdefinition]{./mf6ivar/tex/gwf-sfr-initialstages.dat}
+
 \vspace{5mm}
 \noindent \textit{FOR ANY STRESS PERIOD}
 \lstinputlisting[style=blockdefinition]{./mf6ivar/tex/gwf-sfr-period.dat}

doc/mf6io/ims.tex

@@ -12,7 +12,7 @@ \subsection{Explanation of Variables}
 \input{./mf6ivar/tex/sln-ims-desc.tex}
 \end{description}
 
-\noindent \textbf{IMS Default and Specified Complexity Values}
+\subsection{IMS Default and Specified Complexity Values}
 
 The values that are assigned to the \texttt{nonlinear} and \texttt{linear} variables for the the \texttt{simple}, \texttt{moderate}, and \texttt{complex} complexity options are summarized in table~\ref{table:imsopt}. The values defined for the \texttt{simple} complexity option are assigned if the \texttt{COMPLEXITY} keyword is not specified in the \texttt{OPTIONS} block.
 

doc/mf6io/mf6ivar/dfn/gwf-sfr.dfn

@@ -2,6 +2,22 @@
 # flopy multi-package
 # package-type advanced-stress-package
 
+block options
+name storage
+type keyword
+reader urword
+optional true
+longname activate reach storage
+description keyword that activates storage contributions to the stream-flow routing package continuity equation.
+
+block options
+name storage_weight
+type double precision
+reader urword
+optional true
+longname reach storage time weighting
+description real number value that defines the time weighting factor used to calculate the change in channel storage. STORAGE\_WEIGHT must have a value between 0.5 and 1. Default STORAGE\_WEIGHT value is 0.5.
+
 block options
 name auxiliary
 type string
@@ -641,6 +657,41 @@ reader urword
 longname iprior code
 description character string value that defines the the prioritization system for the diversion, such as when insufficient water is available to meet all diversion stipulations, and is used in conjunction with the value of FLOW value specified in the STRESS\_PERIOD\_DATA section. Available diversion options include:  (1) CPRIOR = `FRACTION', then the amount of the diversion is computed as a fraction of the streamflow leaving reach IFNO ($Q_{DS}$); in this case, 0.0 $\le$ DIVFLOW $\le$ 1.0.  (2) CPRIOR = `EXCESS', a diversion is made only if $Q_{DS}$ for reach IFNO exceeds the value of DIVFLOW. If this occurs, then the quantity of water diverted is the excess flow ($Q_{DS} -$ DIVFLOW) and $Q_{DS}$ from reach IFNO is set equal to DIVFLOW. This represents a flood-control type of diversion, as described by Danskin and Hanson (2002). (3) CPRIOR = `THRESHOLD', then if $Q_{DS}$ in reach IFNO is less than the specified diversion flow DIVFLOW, no water is diverted from reach IFNO. If $Q_{DS}$ in reach IFNO is greater than or equal to DIVFLOW, DIVFLOW is diverted and $Q_{DS}$ is set to the remainder ($Q_{DS} -$ DIVFLOW)). This approach assumes that once flow in the stream is sufficiently low, diversions from the stream cease, and is the `priority' algorithm that originally was programmed into the STR1 Package (Prudic, 1989).  (4) CPRIOR = `UPTO' -- if $Q_{DS}$ in reach IFNO is greater than or equal to the specified diversion flow DIVFLOW, $Q_{DS}$ is reduced by DIVFLOW. If $Q_{DS}$ in reach IFNO is less than DIVFLOW, DIVFLOW is set to $Q_{DS}$ and there will be no flow available for reaches connected to downstream end of reach IFNO.
 
+# --------------------- gwf initial stages ---------------------
+
+block initialstages
+name initialstages
+type recarray ifno initialstage
+shape (maxbound)
+valid
+optional false
+reader urword
+longname
+description
+
+block initialstages
+name ifno
+type integer
+shape
+tagged false
+in_record true
+optional false
+reader urword
+longname reach number for this entry
+description integer value that defines the feature (reach) number associated with the specified initial stage. Initial stage data must be specified for every reach or the program will terminate with an error. The program will also terminate with a error if IFNO is less than one or greater than NREACHES. 
+numeric_index true
+
+block initialstages
+name initialstage
+type double precision
+shape
+tagged false
+in_record true
+optional false
+reader urword
+longname initial reach stage
+description real value that defines the initial stage for the reach. The program will terminate with an error if INITIALSTAGE is less than the RTP value for reach IFNO defined in the PACKAGEDATA block.  INITIALSTAGE data are used only if STORAGE is specified in the Options block and the first stress period is transient or for reaches defined to use the SIMPLE STATUS in the Period block.
+
 
 # --------------------- gwf sfr period ---------------------
 

make/makefile

@@ -505,6 +505,8 @@ $(OBJDIR)/sparsekit.o \
 $(OBJDIR)/rcm.o \
 $(OBJDIR)/blas1_d.o \
 $(OBJDIR)/Iunit.o \
+$(OBJDIR)/GwfSfrCommon.o \
+$(OBJDIR)/gwf-sfr-transient.o \
 $(OBJDIR)/gwf-sfr-steady.o \
 $(OBJDIR)/gwf-sfr-constant.o \
 $(OBJDIR)/RectangularGeometry.o \

msvs/mf6core.vfproj

@@ -261,7 +261,8 @@
 		<Filter Name="GroundWaterFlow">
 		<Filter Name="submodules">
 		<File RelativePath="..\src\Model\GroundWaterFlow\submodules\gwf-sfr-constant.f90"/>
-		<File RelativePath="..\src\Model\GroundWaterFlow\submodules\gwf-sfr-steady.f90"/></Filter>
+		<File RelativePath="..\src\Model\GroundWaterFlow\submodules\gwf-sfr-steady.f90"/>
+	    <File RelativePath="..\src\Model\GroundWaterFlow\submodules\gwf-sfr-transient.f90"/></Filter>
 		<File RelativePath="..\src\Model\GroundWaterFlow\gwf-api.f90"/>
 		<File RelativePath="..\src\Model\GroundWaterFlow\gwf-buy.f90"/>
 		<File RelativePath="..\src\Model\GroundWaterFlow\gwf-chd.f90"/>
@@ -313,6 +314,7 @@
 		<File RelativePath="..\src\Model\ModelUtilities\GwfConductanceUtils.f90"/>
 		<File RelativePath="..\src\Model\ModelUtilities\GwfMvrPeriodData.f90"/>
 		<File RelativePath="..\src\Model\ModelUtilities\GwfNpfOptions.f90"/>
+		<File RelativePath="..\src\Model\ModelUtilities\GwfSfrCommon.f90"/>
 		<File RelativePath="..\src\Model\ModelUtilities\GwfStorageUtils.f90"/>
 		<File RelativePath="..\src\Model\ModelUtilities\GwfVscInputData.f90"/>
 		<File RelativePath="..\src\Model\ModelUtilities\GwtDspOptions.f90"/>