Question regarding SFR package definition

[sebmeng]

Hello Community!
Mybe this is a basic question, but I want to implement SFR package in "real" model so I started from a toy model to test its capabilities. I defined a simple model with 9 cells, and a stream that goes from the upper left cell to the lower right cell:

import flopy
import numpy as np
import math

modelname = 'SFRtest'

# Create the model
sim = flopy.mf6.MFSimulation(sim_name=modelname, version='mf6', exe_name='mf6')

# Create the groundwater flow (gwf) model
gwf = flopy.mf6.ModflowGwf(sim, modelname=modelname)

# Define the grid system
nlay, nrow, ncol = 1, 3, 3
delr = delc = 100.0  # cell size (m)
top = 3  # top of the model
botm = 0  # bottom of the model

# Create the discretization object
dis = flopy.mf6.ModflowGwfdis(gwf, nlay=nlay, nrow=nrow, ncol=ncol, delr=delr, delc=delc, top=top, botm=botm)

# Define the stress periods
nper = 3
perlen = [1.0 for _ in range(nper)]  # length of each stress period
nstp = [1 for _ in range(nper)]  # number of time steps in each stress period
steady = [False for _ in range(nper)]  # transient stress period

# Define the tsmult for each stress period
tsmult = [1.0 for _ in range(nper)]  # length of each stress period

# Create the time discretization object
tdis = flopy.mf6.ModflowTdis(sim, time_units='days', nper=nper, perioddata=[(perlen[i], nstp[i], tsmult[i]) for i in range(nper)])

# Create the initial conditions package
ic = flopy.mf6.ModflowGwfic(gwf, strt=top)

# Define the hydraulic conductivity
hk = 10.0  # m/d

# Create the node property flow package
npf = flopy.mf6.ModflowGwfnpf(gwf, save_flows=True, icelltype=1, k=hk)

# Define the storage properties
ss = 1e-5
sy = 0.05

# Create the storage package
sto = flopy.mf6.ModflowGwfsto(gwf, save_flows=True, iconvert=1, ss=ss, sy=sy, steady_state=[False for _ in range(nper)], transient=[True for _ in range(nper)])

# Create the iterative model solver
ims = flopy.mf6.ModflowIms(sim)


# Calculate the reach length as the diagonal of the cell
rlen = math.sqrt(2) * 100

# Define the other reach parameters
rwid = 1.0  # reach width
rgrd = 0.01  # stream gradient (slope)
rbth = 0.5  # thickness of the streambed
rhk = 1.0  # hydraulic conductivity of the streambed
man = 0.03  # Manning's roughness coefficient
ncon = 1  # number of reaches connected to the reach
ustrf = 1.0  # fraction of upstream flow
ndv = 0  # number of downstream diversions

# Define the top elevation for each cell
top_elevation = [1.0, 1.0, 1.0]  # replace with your actual values

packagedata = [
    [0, (0, 0, 0), rlen, rwid, rgrd, top_elevation[0], rbth, rhk, man, ncon, ustrf, ndv],
    [1, (0, 1, 1), rlen, rwid, rgrd, top_elevation[1], rbth, rhk, man, ncon, ustrf, ndv],
    [2, (0, 2, 2), rlen, rwid, rgrd, top_elevation[2], rbth, rhk, man, ncon, ustrf, ndv]
]

connectiondata = [
    [0, -1],
    [1, 0, -2],
    [2, 1]
]
perioddata = [
    [0, "status","active"],
    [1,"status","active"],
    [2, "status","active"]
]
# Create the ModflowGwfsfr object
sfr = flopy.mf6.ModflowGwfsfr(gwf, pname='sfr', nreaches=3, packagedata=packagedata, connectiondata=connectiondata, perioddata=perioddata)


# Write the simulation to new files and run the simulation
sim.write_simulation()
sim.run_simulation()

The problem is when i run the model, I get this error:

"ERROR REPORT:

1. Reach 3 is connected to reach 2 but reach 2 is not connected
   to reach 3.
   "

According to the docs this is the way to write it, so I don't know where the error could be.

Also an additional question: Can a reach have more than one cell? By the structure of the package i'm guessing not, but maybe I'm wrong.

thanks in advance for the help!
Regards
Sebastian

[emorway-usgs]

Hello @sebmeng,

Regarding your question about a reach having more than one cell: the answer is no. A reach may only be assigned to a single cell. However, a cell may host more than one reach.

For your first question, I think you just need to adjust nconn in the second entry of your sfr packagedata block. Note that the middle reach has two connections:

packagedata = [
    [0, (0, 0, 0), rlen, rwid, rgrd, top_elevation[0], rbth, rhk, man, ncon, ustrf, ndv],
    [1, (0, 1, 1), rlen, rwid, rgrd, top_elevation[1], rbth, rhk, man, 2, ustrf, ndv],
    [2, (0, 2, 2), rlen, rwid, rgrd, top_elevation[2], rbth, rhk, man, ncon, ustrf, ndv]
]

For others coming across this post, I will highlight the mf6 input/output guide which contains this small demonstration (see below). The table on the right of the image provides a nice demonstration. Note that for each reach, upstream reaches still need to be specified even when the connection has already been listed for an upstream reach. For example., note the first two lines in the table of the included image: reach 1 is connected to reach 2 and is preceded with a negative sign since it is a downstream connection. Reach 2 lists reach 1 as a positive value since it is an upstream connection.

[sebmeng]

Thanks for the guidance @emorway-usgs! It worked now!

[emorway-usgs]

Thanks for the small reproducible example, @sebmeng.