Imported from lzwei196/KISS-Knowledge-Infrastructure-for-Scientific-Simulation (models/MODFLOW6/SKILL.md). Install upstream with npx skills add lzwei196/KISS-Knowledge-Infrastructure-for-Scientific-Simulation --skill MODFLOW6. Copyright stays with the author.

MANDATORY EXECUTION POLICY — READ BEFORE PROCEEDING

You MUST run the actual model binary or package described in this document. If the model fails to import, compile, or execute, you MUST:

1. Check diagnostics/triplets.yaml for a matching error pattern
2. Try the fix suggested in the triplet's remedy section
3. If still failing, report the error to the user with full details

You MUST NOT substitute a simplified Python formula, regression equation, or hand-coded approximation in place of the real model. Doing so produces scientifically invalid results and defeats the purpose of the KI.

Before starting, run: python preflight_check.py (in this KI directory) to verify that the model binary/package and required data are available.

DEBUGGING PROTOCOL — When something goes wrong, follow this order:

1. Check triplets — diagnostics/triplets.yaml may already cover this error
2. Read official docs — The model's own documentation for expected formats/units
3. Find working examples — Check outputs/ or the model's shipped test data
4. Fix the tool — With knowledge of what "correct" looks like

Do NOT write custom debug scripts. The answers are in the docs and examples.

KI map — what to read, and when

Projected 2026-08-17 from the KI's actual contents — 10 components present. Refresh: python3 ki_tools_common/generate_skill_map.py --ki_dir <this KI>.

Executable tool index (projected — complete by construction)

Every public tool in this KI, by exact path. What each is FOR lives in the human-written Tool Inventory above; --help on any of these prints its arguments.

28 public tools; _-prefixed helpers and packaging files excluded.

Data Preparation

Input data

Input Source: MODFLOW 6 takes recharge and boundary conditions, not raw met forcing.

• tools/s4/build_rch_package.py — Builds recharge package (from VIC drainage or direct estimation)
• tools/s3/assign_k_from_glhymps.py — Assigns hydraulic conductivity from GLHYMPS 2.0 dataset
• tools/s2/build_layers_from_global.py — Builds layer geometry from DTB China / global datasets
• tools/s4/build_riv_from_cama.py — Builds river boundary from CaMa-Flood network

Data Validation Reference: See data_ki/GLHYMPS/SKILL.md for hydrogeology data documentation. See data_ki/FanWTD/SKILL.md for water table depth documentation. See data_ki/GRACE/SKILL.md for GRACE TWS validation data.

MODFLOW 6 Knowledge Infrastructure — Agent Entry Point

Package: modflow6-knowledge-infrastructure v1.2.0 Model: MODFLOW 6 (USGS Modular Hydrologic Model) Domain: Groundwater hydrology, subsurface flow, solute transport Role in HydroCraft: Fills the groundwater modeling gap — couples with VIC (recharge), CaMa-Flood (river stage), and Lohmann routing (baseflow)

What This Enables

An AI agent can autonomously:

1. Build a MODFLOW 6 groundwater flow model for any basin using FloPy
2. Map VIC deep percolation to MODFLOW recharge (RCH package)
3. Couple CaMa-Flood river stages to MODFLOW river boundaries (RIV package)
4. Feed MODFLOW drain discharge back to routing models as baseflow
5. Extract water table depths for VIC soil moisture feedback
6. Run transient simulations driven by HydroCraft forcing data
7. Add solute transport (GWT) to any existing flow model (nitrate, chloride, etc.)
8. Analyze contaminant plumes, breakthrough curves, and mass budgets
9. Add streamflow routing (SFR) for physics-based stream-aquifer interaction
10. Analyze gaining/losing reaches, stream stage, and downstream flow routing

Quick Reference

6. Output Description

Source of truth: dag.yaml. The dag is the model identity for observable outputs; if this section and dag.yaml ever differ, dag.yaml wins.

Headline output (the dag's validation_rank: 1 variable -- the one this model is judged by):

hydraulic_head -- Hydraulic head per active cell; water table is the head in the topmost active non-dry cell per column. (m)

Other dag outputs currently declared by this KI:

8. Unit Conversion Table

Source of truth: model package expectations, the KI's stage tools, and the existing coupling notes in this skill document. MODFLOW 6 does not convert units internally; every row below must be enforced before writing package files or post-processing outputs.

8c. Sign Conventions and Output Units

11. Validated Results -- Convention Bars

Source of truth: docs/validation_convention.yaml. A metric without its convention band is not a verdict. Null convention bands are written as no cited threshold.

The dag's rank-1 output is hydraulic_head. The currently extracted convention bars provided for this KI are keyed to water_table_elevation:

No cited threshold is stated here for hydraulic_head because the extracted convention facts supplied for this update do not include a hydraulic_head convention band.

Pipeline Overview (10 Stages)

S1 Installation ──> S2 Grid/DIS ──> S3 NPF/STO ──> S4 Boundary Conditions
                                         │                    │
                                         v                    v
                                    S5 TDIS/IC ──────> S6 Solver/IMS
                                                            │
                                                            v
                                                    S7 Execute mf6
                                                            │
                                                            v
                                                S8 Extract Heads/Budget
                                                            │
                                                            v
                                                S9 Postprocess & Visualize
                                                            │
                                                            v
                                          S10 Transport (GWT) [OPTIONAL]
                                          Add GWT model ──> Run ──> Parse
                                                            │
                                                            v
                                          S11 Streamflow Routing (SFR) [OPTIONAL]
                                          Add SFR package ──> Run ──> Parse

Skill Documents (Read for Each Stage)

Tools Reference

Global Data Sources

Automated Setup Pipeline

Three new tools automate the most error-prone setup steps using the global datasets above:

1. build_layers_from_global (S2) — Builds TOP/BOTM/IDOMAIN/STRT from DEM + depth-to-bedrock + water table depth rasters. Eliminates manual layer elevation guesswork.
2. assign_k_from_glhymps (S3) — Assigns K/Sy from GLHYMPS 2.0 for deep layers + HWSD pedotransfer for the shallow soil layer. Handles CRS reprojection and logK decoding automatically.
3. build_riv_from_cama (S4) — Builds RIV package from CaMa-Flood river network + discharge. Auto-clamps rbot to prevent fatal MODFLOW errors.

Critical Domain Knowledge (Non-Obvious Facts)

1. Units are the user's responsibility. MODFLOW 6 does not convert units. If LENGTH_UNITS is "meters" and TIME_UNITS is "days", then K must be m/day, recharge must be m/day, well rates must be m3/day. There is no error if units are wrong — results will simply be silently wrong.

2. Recharge is a rate, not a volume. The RCH package expects recharge in LENGTH/TIME (e.g., m/day), not volume (m3/day). A common error is passing VIC deep percolation as mm/day without converting to m/day (divide by 1000).

3. HDRY sentinel for dry cells. When a cell goes dry (head drops below cell bottom), MODFLOW 6 sets head to 1.0e30 (HDRY). Reading this as a real head value corrupts statistics. Always filter HDRY before computing means or plotting.

4. Layer numbering is 0-indexed in FloPy, 1-indexed in MODFLOW input files. FloPy layer 0 = MODFLOW layer 1. Cellid tuples in FloPy are (layer, row, col) with 0-based indices. In raw MODFLOW input files, they are 1-based.

5. ICELLTYPE matters for water table problems. If the top layer is confined (ICELLTYPE=0) but the water table is in that layer, MODFLOW will not correctly compute saturated thickness — heads can rise above cell top with no error. Set ICELLTYPE=1 for convertible (unconfined) layers.

6. Newton-Raphson formulation is essential for drying/rewetting. Standard formulation can oscillate when cells repeatedly wet and dry. Use NEWTON UNDER_RELAXATION in the GWF model options for unconfined problems.

7. Budget percent discrepancy > 1% indicates a problem. Typical well-posed models have < 0.01% discrepancy. Large discrepancy usually means convergence was not achieved or DVCLOSE/RCLOSE are too loose.

8. The listing file is the primary diagnostic. Budget summaries, convergence iterations, dry cell warnings, and error messages are all in <model_name>.lst. Always read it after a run.

9. FloPy arrays must be 3D (nlay x nrow x ncol). Common error: saving 2D arrays for strt or sy. FloPy will raise "Unable to set data layer 0. Data is not in a valid format" if you pass a 2D (nrow x ncol) array instead of 3D (nlay x nrow x ncol). Always use np.broadcast_to() or np.stack() to ensure 3D shape.

10. Layer 1 K: use VIC/HWSD soil Ksat, NOT GLHYMPS bedrock K. GLHYMPS represents deep bedrock permeability (log-scale, often 1e-15 to 1e-12 m²). For the shallow water table layer (0-10m), VIC soil Ksat from HWSD pedotransfer is more appropriate (typically 0.1-10 m/day for alluvial basins). Tested on Wangjiaba: VIC K=1.0 m/d + Sy=0.15 gave GRACE r=0.56; GLHYMPS K=0.002 m/d + Sy=0.01 gave r=0.41. Use GLHYMPS for Layer 2-3 (deeper bedrock) only. The assign_k_from_glhymps.py tool does this by default — Layer 1 from HWSD, deeper from GLHYMPS — but verify the L1 values are reasonable (>0.1 m/day for alluvial, >0.01 for clay plains).

11. GLHYMPS CRS is Cylindrical Equal Area (meters), not WGS84. Use geopandas bbox reprojection when doing spatial joins: reproject your WGS84 grid centroids to GLHYMPS CRS, or reproject GLHYMPS to WGS84. Direct coordinate comparison will fail silently.

12. GLHYMPS logK values are stored x100. The logK_Ferr_ field contains logK(m^2) multiplied by 100. For example, -1180 means logK = -11.80 m^2. Decode: K_m2 = 10**(logK_field / 100), then convert to m/day: K_m_day = K_m2 * rho * g / mu * 86400.

13. RIV rbot must be above cell BOTM. If river bottom elevation is below the MODFLOW cell bottom, MODFLOW aborts with "RIVER BOTTOM IS LESS THAN CELL BOTTOM". Always clamp: rbot = max(rbot, botm + 0.5).

Error Handling

When errors occur, look up symptoms in diagnostics/triplets.yaml (17 triplets across 7 failure domains). Key triplets:

• dt_mf6_001: Convergence failure — solver did not converge within iteration limit
• dt_mf6_002: Dry cells in unconfined simulation — cells going dry causing oscillation
• dt_mf6_003: Budget imbalance — percent discrepancy > 0.1%
• dt_mf6_004: Recharge unit error — mm/day vs m/day mismatch (silent error)
• dt_mf6_008: Silent wrong recharge — RCH applied to wrong layer or IRCH not set

HydroCraft Coupling

See docs/model_couplings.yaml for detailed coupling specifications:

• VIC deep percolation -> MODFLOW RCH: VIC baseflow_out (mm/day) / 1000 -> m/day recharge
• MODFLOW DRN flux -> Routing baseflow: Drain discharge (m3/day) -> routing input
• CaMa-Flood river stage -> MODFLOW RIV: CaMa sfcelv (m) -> RIV package stage
• MODFLOW water table -> VIC feedback: Water table depth modifies VIC soil moisture capacity

Minimal Working Example (FloPy)

import flopy

# 1. Create simulation
sim = flopy.mf6.MFSimulation(sim_name="test", sim_ws="./workspace")

# 2. Temporal discretization (1 steady-state period)
tdis = flopy.mf6.ModflowTdis(sim, nper=1, perioddata=[(1.0, 1, 1.0)])

# 3. Solver
ims = flopy.mf6.ModflowIms(sim, complexity="SIMPLE")

# 4. Groundwater flow model
gwf = flopy.mf6.ModflowGwf(sim, modelname="gwf", newtonoptions="NEWTON")

# 5. Discretization (10x10 grid, 1 layer, 100m cells)
dis = flopy.mf6.ModflowGwfdis(gwf, nlay=1, nrow=10, ncol=10,
                                delr=100.0, delc=100.0,
                                top=50.0, botm=[0.0])

# 6. Node properties (K=10 m/day, unconfined)
npf = flopy.mf6.ModflowGwfnpf(gwf, icelltype=1, k=10.0)

# 7. Initial conditions
ic = flopy.mf6.ModflowGwfic(gwf, strt=40.0)

# 8. Constant head on left boundary
chd_data = [[(0, i, 0), 45.0] for i in range(10)]
chd = flopy.mf6.ModflowGwfchd(gwf, stress_period_data=chd_data)

# 9. Recharge
rch = flopy.mf6.ModflowGwfrcha(gwf, recharge=0.001)  # 1 mm/day = 0.001 m/day

# 10. Output control
oc = flopy.mf6.ModflowGwfoc(gwf,
        head_filerecord="gwf.hds",
        budget_filerecord="gwf.cbc",
        saverecord=[("HEAD", "ALL"), ("BUDGET", "ALL")])

# 11. Write and run
sim.write_simulation()
sim.run_simulation()

# 12. Read results
import flopy.utils.binaryfile as bf
hds = bf.HeadFile("workspace/gwf.hds")
heads = hds.get_data()

Validated Results — Bengbu Basin (Re-validated 2026-03-22)

Validation Protocol: 3-step (VALIDATION_PROTOCOL.md) Status: production_validated FloPy Version: 3.10.0

Step 1: Binary Test

• MODFLOW 6.6.1 binary works (compiled Intel Fortran, double precision)
• FloPy 3.10.0: run_simulation(exe_name=...) removed, use subprocess.run() instead
• Binary output readable with precision='double'
• Newton + SIMPLE solver = ERROR (CG incompatible with asymmetric matrix) -> use MODERATE/COMPLEX

Step 2: Progressive Data Replacement

Step 3: Full HydroCraft Run Results

Key Lessons Learned

1. HWSD soil Ksat != aquifer K: HWSD gives soil K (~0.01 m/day), aquifer K is ~1-10 m/day. Scale by ~100x for Layer 1.
2. Newton requires MODERATE/COMPLEX solver: SIMPLE uses CG which fails on asymmetric Newton matrices.
3. DRN conductance must use realistic dimensions: C = K * drain_width * cell_width / bed_thickness, NOT K * delr * delc.
4. FloPy 3.10 API changes: run_simulation() no longer accepts exe_name, inner_rclose -> rcloserecord, STO period data format changed.
5. Deep layer dry cells are normal: For unconfined aquifer with WT in Layer 1, Layers 2-3 will be dry. Only check Layer 1 dry cells for validation.
6. Initial heads must be within layer bounds: Clamping WTD to 1-5m prevents starting heads below cell bottom.

Validation Files

• Script: outputs/bengbu_modflow6_revalidation/step2_3_full_bengbu.py
• Results JSON: outputs/bengbu_modflow6_revalidation/validation_results.json
• Plot: outputs/bengbu_modflow6_revalidation/modflow6_bengbu_validation.png
• MODFLOW workspace: outputs/bengbu_modflow6_revalidation/bengbu_workspace/

Stage S10: Groundwater Transport (GWT)

Added: 2026-04-03 KI Version: 1.1.0 Validated: Minimal 1D advection-dispersion test (coupled GWF-GWT)

What GWT Does

The GWT Model simulates three-dimensional transport of a single solute species in flowing groundwater. It solves the advection-dispersion equation using numerical methods and can represent:

1. Advective transport with flowing groundwater
2. Hydrodynamic mechanical dispersion and chemical diffusion
3. Sorption (linear, Freundlich, Langmuir isotherms)
4. First-order and zero-order solute decay/production
5. Mass transfer between mobile and immobile zones (dual-domain)
6. Mixing from GWF stress package sources/sinks

GWT Workflow

Existing GWF Model ──> S10a: configure_gwt.py ──> S10b: Run mf6 ──> S10c: parse_gwt_output.py
     (flow)              (add GWT packages)        (coupled run)       (concentrations, plumes)

Two coupling modes:

1. COUPLED mode (recommended): GWF and GWT in same simulation. The GWF-GWT Exchange automatically passes flows. Both models share the same TDIS.

2. STANDALONE mode: GWT reads previously saved GWF head (.hds) and budget (.cbc) files via the FMI (Flow Model Interface) package. Useful for testing multiple transport scenarios without re-running flow.

Required GWT Packages

Optional GWT Packages

Key Transport Parameters

Contaminant Presets

The configure_gwt.py tool includes documentation-derived parameter presets:

GWF-GWT Exchange (Coupled Mode)

In the simulation name file (mfsim.nam), add:

BEGIN models
  gwf6  gwf.nam  gwf
  gwt6  gwt.nam  gwt
END models

BEGIN exchanges
  GWF6-GWT6  gwf-gwt.gwfgwt  gwf  gwt
END exchanges

BEGIN solutiongroup 1
  ims6  gwf.ims  gwf
  ims6  gwt.ims  gwt
END solutiongroup 1

The exchange file (gwf-gwt.gwfgwt) is empty -- MODFLOW 6 handles the coupling automatically.

GWT Solver Requirements

• MUST use BICGSTAB (not CG) for linear acceleration -- the advection term makes the matrix asymmetric
• Typical convergence: OUTER_DVCLOSE=1e-6, INNER_DVCLOSE=1e-6
• RCLOSE with "strict" option recommended
• Separate IMS for GWT (do not share with GWF IMS)

SSM (Source/Sink Mixing) -- How Concentrations Enter

Three methods to assign concentrations to GWF stress packages:

1. Default: Sources have concentration 0, sinks withdraw at cell concentration
2. Auxiliary variable: Add CONCENTRATION auxiliary to GWF package (e.g., WEL with aux CONCENTRATION), reference in SSM SOURCES block
3. SPC6 file: Separate file with time-varying concentrations per stress period

Example SSM with WEL concentration:

# In GWF: well with concentration auxiliary
wel = flopy.mf6.ModflowGwfwel(gwf, auxiliary="CONCENTRATION",
    stress_period_data=[[(0, 5, 5), -100.0, 50.0]])  # Q=-100, C=50

# In GWT SSM:
ssm = flopy.mf6.ModflowGwtssm(gwt,
    sources=[("WEL-1", "AUX", "CONCENTRATION")])

Tools Reference

Example Usage

# 1. Coupled mode -- add GWT to existing GWF simulation
from tools.s10_transport.configure_gwt import build_coupled_gwt
result = build_coupled_gwt(
    gwf_sim_ws="/path/to/gwf/workspace",
    contaminant="nitrate",
    source_concentration=50.0,       # mg/L
    source_cells=[(0, 5, 5)],        # layer 0, row 5, col 5
    advection_scheme="UPSTREAM",
    run_model=True,
)

# 2. Parse results
from tools.s10_transport.parse_gwt_output import GWTOutputParser
gp = GWTOutputParser(
    concentration_file=result["concentration_file"],
    listing_file=result["listing_file"],
    delr=100.0, delc=100.0,
)
summary = gp.summary()
plume = gp.plume_extent(threshold=1.0, layer=0)
btc = gp.breakthrough_curve(layer=0, row=10, col=10)

CLI Usage

# List contaminant presets
python configure_gwt.py --list-presets

# Build coupled GWT from existing GWF
python configure_gwt.py --gwf-dir /path/to/gwf --contaminant nitrate \
    --source-conc 50.0 --source-cell 0,5,5

# Parse results
python parse_gwt_output.py gwt.ucn --summary --listing-file gwt.lst
python parse_gwt_output.py gwt.ucn --plume 1.0 --delr 100 --delc 100
python parse_gwt_output.py gwt.ucn --breakthrough 0,10,10
python parse_gwt_output.py gwt.ucn --plume-evolution 1.0 --json
python parse_gwt_output.py gwt.ucn --export-csv concentrations.csv

Critical Domain Knowledge for GWT

1. GWT "flows" are mass flows, not water flows. Budget terms in the GWT listing file are in mass/time units, not volume/time.

2. GWT concentration files (.ucn) use text="CONCENTRATION", not "HEAD". When reading with FloPy's HeadFile, pass text="CONCENTRATION".

3. The GWT solver MUST use BICGSTAB, not CG. The advection term produces an asymmetric matrix. CG will fail or produce wrong results.

4. Dispersion XT3D is expensive. Use xt3d_off=True for 1D problems or when flow aligns with the grid. Remove for 2D/3D problems where diagonal flow is important.

5. SSM is required if the GWF model has ANY stress packages, even if all concentrations are zero. Without SSM, MODFLOW 6 will error.

6. GWT has no steady-state option. To find steady-state concentrations, run a long transient simulation until concentrations stabilize.

7. Multiple species require multiple GWT models. Each GWT model simulates one species. Add multiple GWT entries in mfsim.nam, each with its own exchange.

8. Dispersivity is scale-dependent. Field-scale ALH is typically 0.1x to 1x the plume length scale (Gelhar et al., 1992). Lab-scale values are much smaller.

9. MODFLOW 6 GWT does NOT have particle tracking or MOC. Only finite-difference methods (upstream, central, TVD). For sharp fronts, use TVD scheme or refine the grid.

10. HYPE-MODFLOW coupling for nutrients: HYPE nutrient leachate (N, P) can be mapped to MODFLOW GWT via the SSM package. HYPE deep percolation concentration -> RCH auxiliary CONCENTRATION, or via SPC6 files for time-varying concentrations.

HydroCraft Coupling for Transport

Validation Results (Minimal Test)

Stage S11: Streamflow Routing (SFR)

Added: 2026-04-03 KI Version: 1.2.0 Validated: Self-contained test model (10-reach diagonal stream, 10x10 grid)

What SFR Does

The SFR Package simulates 1D open-channel flow in stream networks. Unlike the simpler RIV package (which requires user-specified stream stage), SFR computes stream stage from Manning's equation based on actual stream geometry and flow. This provides:

1. Physically realistic stream-aquifer interaction based on computed stream depth
2. Downstream flow routing through connected reaches
3. Diversions (irrigation canals, water management)
4. Tracking of gaining/losing reaches along the network
5. Proper handling of dry reaches (no-flow conditions)

When to Use SFR vs RIV

Rule of thumb: Use RIV when you have CaMa-Flood or observed stream stages. Use SFR when you need to compute stream stage from flow, route streamflow, or analyze stream-aquifer interaction in detail.

SFR Workflow

Existing GWF Model ──> S11a: configure_sfr.py ──> S11b: Run mf6 ──> S11c: parse_sfr_output.py
     (flow)              (add SFR package)         (simulation)       (stages, exchange, budget)

SFR Key Concepts

1. Reach: A segment of stream within a single model cell. Each reach has static properties (length, width, slope, Manning's n, streambed K, streambed thickness) defined in PACKAGEDATA.

2. Connection: Topology defining which reaches connect to which. Upstream connections are positive; downstream connections are negative in the CONNECTIONDATA block.

3. Period data: Time-varying settings per reach (inflow, rainfall, evaporation, runoff, diversion, status).

4. Manning's equation: Q = (1/n) * A * R^(2/3) * S^(1/2). SFR uses this to compute stream stage from flow. Requires unit conversion factors when using meters + days.

Required SFR Parameters

Stream Type Presets

The configure_sfr.py tool includes documentation-derived presets:

Tools Reference

Example Usage

# 1. Simple linear stream -- auto-builds reaches from cell list
from tools.s11_sfr.configure_sfr import build_simple_sfr
result = build_simple_sfr(
    gwf_sim_ws="/path/to/gwf/workspace",
    stream_cells=[(0,0,0), (0,1,1), (0,2,2), (0,3,3), (0,4,4)],
    stream_elevations=[100.0, 98.0, 96.0, 94.0, 92.0],
    stream_preset="lowland_river",
    inflow_rate=1000.0,     # m3/day
    run_model=True,
)

# 2. Full control -- custom reach data + connections
from tools.s11_sfr.configure_sfr import build_sfr_package
reach_data = [
    {"ifno": 1, "cellid": (0,0,5), "rlen": 500, "rtp": 50.0, "ncon": 1},
    {"ifno": 2, "cellid": (0,1,5), "rlen": 500, "rtp": 49.5, "ncon": 2},
    {"ifno": 3, "cellid": (0,2,5), "rlen": 500, "rtp": 49.0, "ncon": 1},
]
connection_data = [
    [1, -2],      # reach 1 flows downstream to reach 2
    [2, 1, -3],   # reach 2 connects to upstream 1 and downstream 3
    [3, 2],       # reach 3 connects to upstream 2
]
result = build_sfr_package(
    gwf_sim_ws="/path/to/gwf",
    reach_data=reach_data,
    connection_data=connection_data,
    stream_preset="generic",
    period_data={0: [[1, "INFLOW", 500.0]]},
    run_model=True,
)

# 3. Parse results
from tools.s11_sfr.parse_sfr_output import SFROutputParser
sp = SFROutputParser(
    stage_file=result["stage_file"],
    budget_file=result["budget_file"],
    listing_file=result["listing_file"],
)
summary = sp.summary()
exchange = sp.get_stream_aquifer_exchange()
ts = sp.stage_time_series(reach_number=1)

CLI Usage

# List stream presets
python configure_sfr.py --list-presets

# Build self-contained test model
python configure_sfr.py --test --workspace /tmp/sfr_test

# Add SFR to existing GWF model from JSON
python configure_sfr.py --gwf-dir /path/to/gwf --reach-json reaches.json \
    --preset lowland_river

# Parse output
python parse_sfr_output.py --stage gwf.sfr.stg --budget gwf.sfr.cbc \
    --listing gwf.lst --summary
python parse_sfr_output.py --stage gwf.sfr.stg --reach-ts 5
python parse_sfr_output.py --budget gwf.sfr.cbc --exchange
python parse_sfr_output.py --stage gwf.sfr.stg --export-csv stages.csv
python parse_sfr_output.py --stage gwf.sfr.stg --summary --json

Critical Domain Knowledge for SFR

1. SFR uses Manning's equation for stage calculation. Stream stage is NOT user-specified (unlike RIV). SFR computes depth from flow using Q = (1/n) * A * R^(2/3) * S^(1/2), where A = width * depth, R = A / (width + 2*depth).

2. Unit conversion is required for Manning's equation. For meters + days, use length_conversion=1.0 and time_conversion=86400.0. The unit_conversion keyword is deprecated since v6.4.2. Without proper conversion, computed stages will be wildly wrong.

3. Stream gradient (rgrd) must be positive. MODFLOW 6 will terminate if any reach has rgrd <= 0. When building from DEM, clamp minimum slope to ~1e-6.

4. Streambed top (rtp) minus streambed thickness (rbth) must be above cell bottom. Otherwise MODFLOW 6 will abort. Always validate: rtp - rbth > botm[lay, row, col].

5. Connection data uses signed reach numbers. Negative values = downstream, positive = upstream. The first entry in each connection list is the reach number itself.

6. SFR reaches must be connected sequentially. Unlike RIV (independent cells), SFR requires proper network topology. Incorrect connections cause mass balance errors or convergence failures.

7. Dry reaches get DNODATA sentinel (-1e30) in stage output. When a reach goes dry (no flow), MODFLOW sets stage to -1e30. Always filter these before statistics or plotting.

8. SFR is an "advanced" stress package. Unlike simple packages (RIV, CHD), SFR period data persists across stress periods unless explicitly changed. An empty PERIOD block does NOT remove SFR stresses.

9. SFR and RIV should NOT be used in the same cell. Having both packages in the same cell causes double-counting of stream-aquifer exchange. Use one or the other per cell.

10. CaMa-Flood integration: CaMa-Flood can provide river discharge as SFR upstream inflow (INFLOW keyword in period data), or stream widths for reach width. However, CaMa-Flood does not provide cross-sections or Manning's n -- those must come from literature or field data.

HydroCraft Coupling for SFR

Validation Results (Test Model)

This knowledge infrastructure was built using the Knowledge Dissection Toolkit v1.0 (Jianyun Zhang Research Group, Hohai University).