# nwsrfsr


R package interface to the NWSRFS (National Weather Service River
Forecast System) Fortran hydrologic models: SAC-SMA, SNOW-17, Unit
Hydrograph, Lag-K, Chanloss, and Consuse.

Provides both [low-level Fortran
wrappers](#running-individual-model-components) and a [high-level
orchestration layer](#quick-start) that reads NWRFC autocalibration
directories, auto-detects model components, and runs the full model
chain. Also includes an [AdjustQ module](#adjustq-preprocessing) for
preprocessing upstream flow.

[Full
documentation](https://noaa-nwrfc.github.io/nwsrfs-hydro-models/r/)

## Installation

Using [pixi](https://pixi.prefix.dev/) (recommended for development —
provides R, gfortran, and test dependencies):

``` bash
pixi run install-r
```

See the [top-level
README](https://github.com/NOAA-NWRFC/nwsrfs-hydro-models#development-environment-pixi)
for pixi setup details.

Or from R:

``` r
# install.packages("devtools")
devtools::install_github("NOAA-NWRFC/nwsrfs-hydro-models", subdir = "nwsrfs_r")
```

Or from a local clone:

``` bash
R CMD INSTALL nwsrfs_r
```

Requires `gfortran` for compiling the Fortran source.

## Quick Start

[`load_example()`](https://noaa-nwrfc.github.io/nwsrfs-hydro-models/r/reference/load_example.html)
runs the full model chain on bundled data and returns an `nwsrfs_run`
object:

``` r
library(nwsrfsr)

# NRKW1: Nooksack River — 2 SAC-SMA/SNOW17 zones + 3 Lag-K upstream tribs
run = load_example("NRKW1")
print(run)
#> NWSRFS Model Run
#>   Zones:     NRKW1-1, NRKW1-2
#>   Timestep:  6 hours
#>   Length:    62820 timesteps
#>   Components:
#>     SAC-SMA/SNOW17/UH: TRUE
#>     Lag-K:             TRUE
#>     Chanloss:          FALSE
#>     Consuse:           FALSE

# SFLN2: Salmon Falls Creek — 2 zones + CU zone, chanloss, consumptive use
run2 = load_example("SFLN2")
print(run2)
#> NWSRFS Model Run
#>   Zones:     SFLN2-1, SFLN2-2, SFLN2-CU
#>   Timestep:  6 hours
#>   Length:    62828 timesteps
#>   Components:
#>     SAC-SMA/SNOW17/UH: TRUE
#>     Lag-K:             FALSE
#>     Chanloss:          TRUE
#>     Consuse:           TRUE
```

The `$sim` element is the simulated streamflow vector (cfs, 6-hour
timestep):

``` r
summary(run$sim)
#>    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.
#>   561.6  1894.5  2888.4  3590.3  4302.7 52279.9

plot(run$sim, type = "l", ylab = "Flow (cfs)", main = "NRKW1 Simulated Flow")
```

You can also run from an NWRFC autocalibration directory on disk via
[`nwsrfs_run()`](https://noaa-nwrfc.github.io/nwsrfs-hydro-models/r/reference/nwsrfs_run.html):

``` r
# example, assuming you have downloaded the sameple data: https://doi.org/10.5281/zenodo.18829935
run = nwsrfs_run("nwsrfs-ac-sample-runs/1zone/FSSO3")
```

Sample autocalibration data for 5 basins is available from the
[autocalibration
repository](https://github.com/NOAA-NWRFC/nwsrfs-hydro-autocalibration)
and archived on [Zenodo (DOI:
10.5281/zenodo.18829935)](https://doi.org/10.5281/zenodo.18829935).

## Exploring Bundled Data

The package ships example data for two stations. The parameter tables
have columns `p_name` (unique ID), `name`, `type` (model component),
`zone`, and `value`:

``` r
head(nrkw1_pars)
#>            p_name      name type  zone       value
#> 1   init_co_MFNW1   init_co lagk MFNW1 247.4027145
#> 2   init_if_MFNW1   init_if lagk MFNW1 111.5698487
#> 3   init_of_MFNW1   init_of lagk MFNW1 391.9661841
#> 4 init_stor_MFNW1 init_stor lagk MFNW1 113.7500227
#> 5    ktbl_a_MFNW1    ktbl_a lagk MFNW1  -0.7160900
#> 6    ktbl_b_MFNW1    ktbl_b lagk MFNW1  -0.1514212
```

Forcing data is a named list of data.frames, one per zone, with 6-hour
timesteps:

``` r
names(nrkw1_forcing)
#> [1] "NRKW1-1" "NRKW1-2"

head(nrkw1_forcing[["NRKW1-1"]])
#>   year month day hour map_mm mat_degc ptps
#> 1 1979    10   1    0      0    9.278    0
#> 2 1979    10   1    6      0   12.759    0
#> 3 1979    10   1   12      0   17.602    0
#> 4 1979    10   1   18      0   10.296    0
#> 5 1979    10   2    0      0    7.491    0
#> 6 1979    10   2    6      0   10.907    0
```

Upstream flow data (NRKW1 only) is a named list of data.frames for each
Lag-K reach:

``` r
names(nrkw1_upflow)
#> [1] "MFNW1" "NFNW1" "NSSW1"

head(nrkw1_upflow[["MFNW1"]])
#>   year month day hour flow_cfs
#> 1 1979    10   1    0 209.3404
#> 2 1979    10   1    6 209.3404
#> 3 1979    10   1   12 202.4441
#> 4 1979    10   1   18 179.9215
#> 5 1979    10   2    0 181.8996
#> 6 1979    10   2    6 189.0886
```

## Running Individual Model Components

The model chain can be run step-by-step using the low-level wrappers.
This follows the same sequence the high-level
[`load_example()`](https://noaa-nwrfc.github.io/nwsrfs-hydro-models/r/reference/load_example.html)
uses internally: FA -\> SAC-SMA/SNOW17 -\> UH -\> Lag-K.

``` r
library(nwsrfsr)
dt_hours = 6

# 1. Apply forcing adjustments (monthly climatological corrections + PET)
forcing_adj = fa_nwrfc(dt_hours, nrkw1_forcing, nrkw1_pars)
head(forcing_adj[["NRKW1-1"]])
#>   year month day hour map_mm  mat_degc ptps    pet_mm    etd_mm
#> 1 1979    10   1    0      0  8.063764    0 0.5458438 0.5108825
#> 2 1979    10   1    6      0 11.544764    0 0.5458438 0.5108825
#> 3 1979    10   1   12      0 16.387764    0 0.5458438 0.5108825
#> 4 1979    10   1   18      0  9.081764    0 0.5458438 0.5108825
#> 5 1979    10   2    0      0  6.275488    0 0.6351838 0.5918346
#> 6 1979    10   2    6      0  9.691488    0 0.6351838 0.5918346

# 2. SAC-SMA / SNOW-17 — returns total channel inflow (mm) per zone
tci = sac_snow(dt_hours, forcing_adj, nrkw1_pars)
dim(tci)
#> [1] 62820     2
head(tci)
#>           [,1]      [,2]
#> [1,] 0.2462325 0.9752563
#> [2,] 0.2488699 0.9984781
#> [3,] 0.2476989 0.9733128
#> [4,] 0.2465621 0.9494683
#> [5,] 0.2454517 0.9268678
#> [6,] 0.2443759 0.9054463

# 3. Unit Hydrograph routing — converts TCI to streamflow (cfs)
local_flow = uh(dt_hours, tci, nrkw1_pars,
                sum_zones = TRUE, start_of_timestep = TRUE, backfill = TRUE)
head(local_flow)
#> [1] 485.3301 485.3301 629.9145 661.2546 661.6406 653.3936

# 4. Lag-K routing of upstream tributaries
lagk_flow = lagk(dt_hours, nrkw1_upflow, nrkw1_pars, sum_routes = TRUE)
head(lagk_flow)
#> [1] 887.7466 979.0988 672.8225 674.2631 635.2833 633.9531

# 5. Combine local + routed upstream flow
sim = local_flow + lagk_flow
head(sim)
#> [1] 1373.077 1464.429 1302.737 1335.518 1296.924 1287.347
```

See the API reference for each function:
[`fa_nwrfc()`](https://noaa-nwrfc.github.io/nwsrfs-hydro-models/r/reference/fa_nwrfc.html),
[`sac_snow()`](https://noaa-nwrfc.github.io/nwsrfs-hydro-models/r/reference/sac_snow.html),
[`uh()`](https://noaa-nwrfc.github.io/nwsrfs-hydro-models/r/reference/uh.html),
[`lagk()`](https://noaa-nwrfc.github.io/nwsrfs-hydro-models/r/reference/lagk.html),
[`chanloss()`](https://noaa-nwrfc.github.io/nwsrfs-hydro-models/r/reference/chanloss.html),
[`consuse()`](https://noaa-nwrfc.github.io/nwsrfs-hydro-models/r/reference/consuse.html)

<!-- TODO Add autocalb documentation
Updating Parameters
-------------------
&#10;[`update_pars()`](https://noaa-nwrfc.github.io/nwsrfs-hydro-models/r/reference/update_pars.html)
modifies parameters in an existing run and re-runs the model chain. Pass a data.frame
with `p_name` and `value` columns:
&#10;```r
run = load_example("NRKW1")
&#10;# See current uztwm (upper zone tension water max) values
nrkw1_pars[nrkw1_pars$name == "uztwm", ]
#>            p_name  name type    zone    value
#> 199 uztwm_NRKW1-1 uztwm  sac NRKW1-1 25.12025
#> 346 uztwm_NRKW1-2 uztwm  sac NRKW1-2 25.05820
&#10;# Update zone 1 uztwm and re-run
new_pars = data.frame(p_name = "uztwm_NRKW1-1", value = 120.0)
run2 = update_pars(run, new_pars)
&#10;mean(run$sim)   # 3590.3
mean(run2$sim)  # 3572.0
```
-->

## AdjustQ Preprocessing

[AdjustQ](https://noaa-nwrfc.github.io/nwsrfs-hydro-models/r/reference/adjustq.html)
replicates the CHPS FEWS AdjustQ tool, creating 6-hour upstream flow
time series from daily and instantaneous observations, optionally
blended with simulation output.

The quickest way to run it is with
[`adjustq_load_example()`](https://noaa-nwrfc.github.io/nwsrfs-hydro-models/r/reference/adjustq_load_example.html),
which uses bundled NRKW1 data:

``` r
# With simulation (runs load_example("NRKW1") internally)
aq = adjustq_load_example(sim = TRUE)
head(aq)
#>              datetime flow_cfs
#> 1 1979-10-01 00:00:00 1089.111
#> 2 1979-10-01 06:00:00 1161.571
#> 3 1979-10-01 12:00:00 1033.318
#> 4 1979-10-01 18:00:00 1059.320
#> 5 1979-10-02 00:00:00 1032.968
#> 6 1979-10-02 06:00:00 1029.565

# Without simulation (observation-only merge)
aq2 = adjustq_load_example(sim = FALSE)
```

To call
[`adjustq()`](https://noaa-nwrfc.github.io/nwsrfs-hydro-models/r/reference/adjustq.html)
manually with your own data:

``` r
run = load_example("NRKW1")

# Build a sim data.frame with year/month/day/hour/flow_cfs columns
sim_df = data.frame(
  run$forcings[[1]][, c("year", "month", "day", "hour")],
  flow_cfs = run$sim
)

result = adjustq(
  daily_flow = nrkw1_daily_flow,
  inst_flow = nrkw1_inst_flow,
  sim = sim_df
)
```

## Bundled Datasets

| Object | Description |
|----|----|
| [`nrkw1_pars`](https://noaa-nwrfc.github.io/nwsrfs-hydro-models/r/reference/nrkw1_pars.html) | NRKW1 optimal parameters (366 rows) |
| [`nrkw1_forcing`](https://noaa-nwrfc.github.io/nwsrfs-hydro-models/r/reference/nrkw1_forcing.html) | NRKW1 6-hour forcing, 2 zones |
| [`nrkw1_upflow`](https://noaa-nwrfc.github.io/nwsrfs-hydro-models/r/reference/nrkw1_upflow.html) | NRKW1 upstream flow, 3 Lag-K reaches (MFNW1, NFNW1, NSSW1) |
| [`nrkw1_daily_flow`](https://noaa-nwrfc.github.io/nwsrfs-hydro-models/r/reference/nrkw1_daily_flow.html) | NRKW1 observed daily streamflow |
| [`nrkw1_inst_flow`](https://noaa-nwrfc.github.io/nwsrfs-hydro-models/r/reference/nrkw1_inst_flow.html) | NRKW1 observed 6-hour instantaneous flow |
| [`sfln2_pars`](https://noaa-nwrfc.github.io/nwsrfs-hydro-models/r/reference/sfln2_pars.html) | SFLN2 optimal parameters (514 rows) |
| [`sfln2_forcing`](https://noaa-nwrfc.github.io/nwsrfs-hydro-models/r/reference/sfln2_forcing.html) | SFLN2 6-hour forcing, 3 zones (incl. CU) |
| [`sfln2_daily_flow`](https://noaa-nwrfc.github.io/nwsrfs-hydro-models/r/reference/sfln2_daily_flow.html) | SFLN2 observed daily streamflow |
| [`sfln2_inst_flow`](https://noaa-nwrfc.github.io/nwsrfs-hydro-models/r/reference/sfln2_inst_flow.html) | SFLN2 observed 6-hour instantaneous flow |

**NRKW1** — Nooksack River at North Cedarville, WA (USGS-12210700).
2-zone SAC-SMA/SNOW17/UH model with 3-reach Lag-K routing, no chanloss
or consuse.

**SFLN2** — Salmon Falls Creek NR San Jacinto, NV (USGS-13105000).
2-zone SAC-SMA/SNOW17/UH model with chanloss and consumptive use, no
upstream flow.

## Citation

``` r
citation("nwsrfsr")
```

This includes the software citation and the accepted 2026 JAWRA paper
reference. Repository-level citation metadata is in `/CITATION.cff`
(repo root).

## Testing & Building

``` bash
pixi run test-r          # run R tests (testthat)
pixi run check-r         # R CMD check (no manual)
pixi run build-r         # build R source tarball
pixi run check-r-cran    # CRAN submission check on built tarball
```

Or from R:

``` r
devtools::test("nwsrfs_r")
```

## Documentation

[**Live documentation
site**](https://noaa-nwrfc.github.io/nwsrfs-hydro-models/r/)

Build locally with pixi:

``` bash
pixi run build-docs-r          # build pkgdown site
pixi run build-docs-r-readme   # render README.qmd to HTML
pixi run build-docs            # build all docs (Python + R + R README)
```

Or from R:

``` r
pkgdown::build_site(pkg = "nwsrfs_r", install = FALSE, preview = FALSE)
```

## Legal

- Package license: `LICENSE`
- Package notice: `NOTICE`
