Data and code for: Endangered salmon juvenile synchrony signals risk of extirpation

Ryan, Rachael E.1 ; Carlson, Stephanie 1 ; Grantham, Theodore1; Ettlinger, Eric2; Reichmuth, Michael3; Ruhi, Albert 1

Published Sep 03, 2026 on Dryad. https://doi.org/10.5061/dryad.rfj6q57p6

Data files

Sep 03, 2026 version files 64.99 MB

lagunitas_data_2026-09-01.zip

64.92 MB
README.md

66.10 KB

Abstract

Understanding whether small, imperiled population complexes retain buffering capacity is central to advancing species recovery and ensuring persistence. Asynchrony among sub-populations can buffer population complexes against environmental disturbance, but it remains poorly understood whether asynchronous dynamics provide similar buffering within imperiled (depressed) population complexes. We examined this question using the case study of an endangered coho salmon population complex in the Lagunitas Basin in central California. We analyzed 18 years of fish and flow monitoring data across four streams in the watershed and assessed long-term fluctuations in juvenile densities, what these dynamics mean for watershed-scale extinction risks, and possible drivers. Multivariate autoregressive state-space (MARSS) models revealed that variation was 19 times higher on average in sub-populations in tributary streams compared to the mainstem; however, accounting for observation error revealed that temporal fluctuations in juvenile density were highly synchronous among the sub-populations. Simulations showed that as a result of these highly correlated dynamics across the small sub-populations, the watershed is at heightened risk of collapse. Thus, efforts to diversify dynamics, in addition to increasing abundance, would help reduce risk of simultaneous collapse. Finally, adult abundance and hydroclimate explained some variation in juvenile density, but most shared process variance was unexplained. Our findings underscore the critical need to understand the spatial patterns and drivers of population dynamics, even in small, endangered population complexes, to promote recovery and resilience. Small population complexes with a high degree of synchrony among sub-populations may be particularly susceptible to collapse due to simultaneous declines of sub-populations with already low abundances.

Related publication: Ryan et al. Endangered salmon juvenile synchrony signals risk of extirpation. Ecological Applications.

Dataset DOI: 10.5061/dryad.rfj6q57p6

Description of the data

This dataset supports time series analysis of juvenile coho salmon dynamics in the Lagunitas Creek watershed, Marin County, California, using long-term monitoring data from Marin Water and the National Park Service. The goal of the analysis is to evaluate subpopulation dynamics within this watershed to understand synchrony, investigate how synchrony contributes to quasi-extinction risk, and assess how spawner abundance and hydroclimatic variation contribute to juvenile density dynamics.

We used multivariate autoregressive state-space (MARSS) models, synchrony estimates, covariate analyses, and quasi-extinction analyses to ask three questions:

Do juvenile coho salmon subpopulations in the Lagunitas Creek watershed show evidence of population diversity, defined by differences in mean density, temporal variability, and synchrony among streams?
Given the estimated variance and covariance in juvenile density dynamics, what is the risk of simultaneous multi-stream declines over 10-, 15-, and 20-year horizons, and how much does each stream contribute to buffering watershed-scale decline risk?
How do spawner abundance and hydroclimatic variation influence juvenile density dynamics and inferred process synchrony among streams?

The response variable is late-summer juvenile coho salmon density, calculated as the number of juvenile fish per square metre of sampled pool habitat within consistently sampled index reaches from 2005 to 2022. The analysis also uses annual redd counts as a proxy for spawner input, daily streamflow and functional flow metrics to represent hydroclimatic conditions, and spatial watershed data to describe the focal streams.

Analytical approach

To answer these questions, we used long-term monitoring data from annual electrofishing surveys in the four streams that consistently support juvenile coho salmon rearing in the Lagunitas Creek watershed: Lagunitas Creek, the regulated mainstem below Peters Dam; and three tributaries, Devil's Gulch, San Geronimo Creek, and Olema Creek. Juvenile sampling was conducted by Marin Water in Lagunitas Creek, Devil's Gulch, and San Geronimo Creek, and by the National Park Service in Olema Creek. We focused on pool habitat sampled during August and September and used 100-300 m index reaches that were sampled consistently across years. Within each reach-year, juvenile coho counts from sampled pools were divided by the total sampled pool surface area to produce fish/m^2 of pool habitat.

We summarized annual coho "redd" (gravel nest) counts for each focal stream as a proxy for spawner abundance and egg production. Redd data were compiled from Marin Water and National Park Service spawning ground surveys and aligned to the juvenile time series by water year. We also compiled daily streamflow records and annual functional flow metrics to describe hydroclimatic variation. Annual functional flow metrics were reduced with principal components analysis, and the retained hydroclimatic covariate was used in the MARSS covariate models. Because Devil's Gulch does not have a long-term gauge, San Geronimo Creek flow was used as the hydroclimatic proxy for Devil's Gulch.

Briefly, our analysis did the following:

We used a multivariate autoregressive state-space (MARSS) model to estimate density dynamics of the four subpopulations, removing observation error to get process error variance and correlation. We tested 5 different structures of observation error, along with two different configurations for drift, and selected the best model.
We then estimated synchrony using Pearson correlations of raw data, modeled data, and correlations in process error.
Over 10,000 simulations, we forecasted the risk of simultaneous declines in density of 2+, 3+, and 4 streams across 10-, 15-, and 20-year time horizons using the base MARSS model process equation. The main decline threshold was 90%, with additional 80% and 95% thresholds used as sensitivity analyses.
Finally, we ran MARSS models with biological (redds) and environmental (hydroclimate) covariates to assess how these factors explain density dynamics (variance, synchrony).

Archive contents and reproducibility

The archive lagunitas_data_2026-09-01.zip contains five principal components and extracts directly into the destination folder selected by the user:

data/raw/: monitoring-agency extracts and source spatial or hydroclimatic data.
data/processed/: current derived datasets used by the manuscript analysis and figure code.
outputs/model/: current MARSS and PVA objects and tabular results.
outputs/tables/manuscript_2026-06-11/: CSV data underlying manuscript figures and tables.
scripts/: seven commented R Markdown scripts documenting processing, analysis, and manuscript-output generation.

Sensitive geographic data

Biological location data were generalized in accordance with data-sharing guidance for Near Threatened species. All published latitude and longitude values regarding fish locations have been rounded to 0.01 degrees. Public environmental base layers, including NHDPlus hydrography and catchments and the Marin County geology layer, are not species-occurrence records and remain at their source resolution.

Project organization

Data directories are numbered to match the script that first processes or uses them. data/raw/ contains monitoring-agency extracts and original spatial or hydroclimatic source files; data/processed/ contains cleaned or summarized analysis inputs. outputs/model/ contains fitted model objects and exported model or PVA results. outputs/tables/ contains tabular values underlying the manuscript's figures and tables. Published figure images are not deposited, but the code and data needed to regenerate them are retained.

(archive root)/
|- lagunitas-MARSS.Rproj
|- data/
|  |- raw/
|  |  |- 01_juveniles_data_raw/
|  |  |- 02_redds_data_raw/
|  |  |- 03_flow_data_raw/
|  |  |- 06_manuscript_figures_data_raw/
|  |  `- 07_stream_attributes_data_raw/
|  |     |- prism_ppt_annual/
|  |     `- prism_tmean_annual/
|  `- processed/
|     |- 01_juveniles_data_processed/
|     |- 02_redds_data_processed/
|     `- 03_flow_data_processed/
|        |- ffc_upload/
|        `- ffc_results/
|- outputs/
|  |- model/
|  |  `- pva/
|  `- tables/
|     `- manuscript_2026-06-11/
`- scripts/

Script order

scripts/01_juveniles_2026-05-05.Rmd: prepare juvenile counts, habitat areas, and reach-level densities.
scripts/02_redds_2026-05-05.Rmd: prepare annual streamwide redd counts.
scripts/03_flow_2026-05-05.Rmd: reconstruct missing San Geronimo discharge, process Functional Flows Calculator outputs, and calculate hydroclimatic principal components.
scripts/04_marss_model_2026-05-04.Rmd: fit and compare MARSS models and export model summaries.
scripts/05_pva_2026-05-05.Rmd: run joint-decline and buffering simulations.
scripts/06_manuscript_figures_2026-06-29.Rmd: generate manuscript figures and manuscript-ready CSV tables.
scripts/07_stream_attributes_2026-06-12.Rmd: calculate stream and watershed attributes.

Scripts 04–06 use analysis_date <- "2026-06-11" for the retained model, PVA, and manuscript-table outputs. The current processed CSVs are included as the exact inputs used for those analyses. In scripts 01–03, some write_csv() calls are commented to protect the deposited processed files from accidental overwriting; uncomment the applicable export calls to rebuild those files from the raw inputs.

Most inputs are included locally. The study-system map portion of script 06 additionally retrieves 2024 California county and water layers through the tigris R package, so that portion requires an internet connection and access to the U.S. Census Bureau web service.

Conventions used throughout the archive

Stream and organization abbreviations

DG: Devil's Gulch.
LG: Lagunitas Creek.
SG: San Geronimo Creek.
OL: Olema Creek.
RD: Redwood Creek; retained only in source NPS data and excluded from the focal MARSS analysis.
MW: Marin Water.
NPS: National Park Service.
juv: juvenile coho salmon.
redd: a salmonid gravel nest.
Efish orefish: electrofishing.
MARSS: multivariate autoregressive state-space model.
PVA: population viability analysis.
FFC: Functional Flows Calculator.
FFM or ffm: functional-flow metric.
PC1: first principal component.
Q: MARSS process-error covariance matrix.
R: MARSS observation-error covariance matrix.
U: MARSS drift vector.
C: MARSS covariate-effect matrix.
AICc: Akaike information criterion corrected for small sample size.
CI: confidence interval.
SE: standard error.
CV: coefficient of variation.
kNN: k-nearest-neighbor classification.
COMID or comid: unique National Hydrography Dataset Plus flowline identifier.
NHD+: National Hydrography Dataset Plus.
PRISM: Parameter-elevation Regressions on Independent Slopes Model climate dataset.
NLCD: National Land Cover Database.

Dates and time indices

FieldSeason is the calendar year of late-summer juvenile sampling.
WaterYear is the water year ending in the named calendar year; California water years run from October 1 through September 30.
SpawnerSeason is written as two calendar years, such as and maps to the water year identified by the second year.
In MARSS state files, t = 1 corresponds to 2005, t = 2 to 2006, and so forth through t = 18 for 2022.
year_ahead is the number of simulated years after the fixed PVA starting state.
Four-digit column headings in wide files are years or water years, as stated for each file group below. The cell is the value of the row variable for that year.

Measurement units and analytical scales

Juvenile length: millimetres (mm) where explicitly labeled Length coho (mm).
Fish weight: grams (g).
Habitat length and width: metres (m).
Habitat area: square metres (m²).
Juvenile linear density: fish per metre of sampled pool habitat (fish/m).
Juvenile areal density: fish per square metre of sampled pool habitat (fish/m²).
Daily discharge and FFC flow-magnitude metrics: cubic feet per second (cfs).
FFC timing metrics: water-year day, with October 1 equal to day 1.
FFC duration metrics: days.
FFC frequency metrics: occurrences per season.
FFC spring rate of change: percent change per day, reported by the FFC as percent.
Geographic Latitude and Longitude values in the NPS CSV files are decimal degrees rounded to 0.01 degrees. Inconsistent, unused CoordUnits and CoordSystem fields inherited from the NPS database export were omitted from the deposited CSVs to prevent the coordinates from being misinterpreted.
MARSS latent states are fitted as log(density + 1), where density is fish/m².
MARSS process and observation variances are on the squared log(density + 1) scale; covariances are on the corresponding cross-product scale. Correlations, probabilities, AICc weights, and coefficients of variation are unitless.
PVA risk fields in outputs/model/pva/ are proportions from 0 to 1. Manuscript-ready tables may format the same quantities as percentage strings.

Missing and special values

Depending on the source file, missing values may be blank, NA, lowercase na, or None. In the FFC annual-result matrices, None means the FFC could not calculate that metric for that water year. In processed redd covariates, a zero is a retained count of zero focal coho redds after the documented filtering rules, not a missing value.

Tabular data dictionary

Variables repeated across files are defined once below and apply wherever the same column name is used. Differences in capitalization or presentation labels are listed explicitly as aliases.

General identifiers, dates, and uncertainty fields

Variable	Description	Units or values
`Stream`, `stream`	Stream name or code associated with the observation or result.	Names or `DG`, `LG`, `SG`, `OL`; some raw NPS data also contain Redwood Creek.
`Reach`	Consistently identified juvenile index reach. Examples include `DG-1`, `LG-12`, `SG-2`, and `OL1`.	Identifier; no units.
`Site`	Row label for a stream-level or regional series in a wide summary.	Identifier; no units.
`FieldSeason`	Calendar year of juvenile field sampling.	Year.
`WaterYear`	Water year labeled by its ending calendar year.	Year; October 1–September 30.
`year`, `Year`	Calendar year or water year. In FFC annual-result matrices, the column named `Year` actually contains metric names; the four-digit columns contain water years.	Year, except for the FFC matrix row-label column.
`Date`, `date`	Observation date or daily-flow date. ISO `yyyy-mm-dd` is used in cleaned daily data; FFC upload files use `mm/dd/yyyy`.	Date.
`t`	Sequential MARSS time index.	Integer 1–18, corresponding to 2005–2022.
`year_ahead`	PVA simulation year after the starting state.	Years.
`horizon`	Maximum PVA simulation horizon.	10, 15, or 20 years.
`Model`, `model`	Model name or model-structure label.	Categorical text.
`estimate`	Estimated quantity. In state files it is latent `log(density + 1)`; in coefficient files it is the fitted MARSS parameter identified by `term` and `parameter`.	Context-dependent fitted scale.
`se`, `SE`, `std.error`	Standard error of the associated estimate.	Same scale as the estimate.
`conf.low`, `Lower 95% CI`	Lower bound of a 95% confidence interval.	Same scale as the estimate.
`conf.high`, `conf.up`, `Upper 95% CI`	Upper bound of a 95% confidence interval.	Same scale as the estimate.

Raw juvenile monitoring files

This dictionary covers:

data/raw/01_juveniles_data_raw/MW_juveniles_2005-2022.csv
data/raw/01_juveniles_data_raw/MW_juveniles_pool_2021_2022.csv
data/raw/01_juveniles_data_raw/mw_index_reach_habitat_modified.csv
data/raw/01_juveniles_data_raw/Olema_Redwood_2005_2022_juveniles.csv

Variable	Description	Units or values
`EventID`	NPS source-database identifier for a sampling event.	Integer identifier.
`LocationID`	NPS source-databaseidentifier for the specific sampling location within a stream.	Integer identifier.
`StreamID`	NPS source-database stream identifier. Scripts group IDs 5, 9, 16, and 17 as Olema Creek and IDs 34, 35, and 37 as Redwood Creek. Only Olema records are carried into the focal juvenile analysis.	Integer code.
`SectionNum`	NPS source-database monument-section number used to define repeatable stream sections; retained sections are converted to reach labels in script 01.	Integer identifier.
`Habitat`, `Habitat_Unit`, `Hab`	NPS source-database habitat description or standardized habitat grouping. `Hab` is the simplified value `pool` used in processing.	Categorical text.
`UnitType`	NPS source-database habitat-unit type code. Codes present are `BW` = backwater pool, `MC` = mid-channel/main-channel pool, `PL` = plunge pool, and `SC` = scour pool. The original NPS Access Database codebook was not available for deposit, so the `MC` wording records the unresolved naming variant; script 01 groups all four codes as pool habitat.	Categorical code.
`Length coho (mm)`	Fork or total length recorded for an individual coho salmon in the Marin Water source extract.	Millimetres.
`Weight Coho (g)`	Weight recorded for an individual coho salmon.	Grams.
`Length`	Length of the sampled NPS habitat unit.	Metres.
`AverageWidth`	Average width of the sampled NPS habitat unit.	Metres.
`Pool_length_m`	Length of an individual Marin Water pool habitat unit.	Metres.
`Pool_width_m`	Width of an individual Marin Water pool habitat unit.	Metres.
`total_length_m`	Total sampled pool length within a reach-year.	Metres.
`avg_width_m`	Average sampled pool width within a reach-year.	Metres.
`Pass`	Sequential electrofishing pass number.	Integer, normally 1–3.
`SpeciesCode`	Source fish-species code. The retained juvenile extract contains `CO` = coho salmon. Other NPS salmonid codes described in related redd data are `CH` = Chinook salmon and `SH` = steelhead trout.	Categorical code.
`LifeStage`	Life stage of the sampled fish. `YoY` means young-of-year.	Categorical text or code.
`Count`	Number of fish recorded for the NPS event, habitat unit, species, life stage, and electrofishing pass represented by the row.	Fish count.
`no_fish`	Total number of juvenile coho counted for the Marin Water reach-year row.	Fish count.

Processed juvenile files

This dictionary covers all CSV files retained in data/processed/01_juveniles_data_processed/.

Variable	Description	Units or values
`fish_m_pool_reach`	Observed or depletion-estimated juvenile coho abundance divided by total sampled pool length within the reach-year.	Fish/m.
`fish_m2_pool_reach`	Observed or depletion-estimated juvenile coho abundance divided by estimated sampled pool surface area within the reach-year. This is the reach-level response used to construct the MARSS input.	Fish/m².
`stream_fish`	Total juvenile coho counted or estimated across retained reaches for the stream-year.	Fish count.
`stream_area`	Total sampled pool surface area across retained reaches for the stream-year.	m².
`regional_fish`	Total juvenile coho counted or estimated across all four focal streams for the year.	Fish count.
`regional_area`	Total sampled pool surface area across all four focal streams for the year.	m².
`stream_density`	`stream_fish / stream_area`.	Fish/m².
`regional_density`	`regional_fish / regional_area`.	Fish/m².
`first_sample`	First juvenile sampling date in the stream and field season.	Date.
`num_yrs`	Number of field seasons in which the NPS section was sampled; used to identify repeat sections.	Years sampled, integer count.
Four-digit columns `2005`–`2022` in `juv_m2_wide_reduced.csv` and `juv_wide_m2.csv`	Reach-level `fish_m2_pool_reach` value for the named field season. `juv_m2_wide_reduced.csv` contains the ten reaches retained in the MARSS analysis.	Fish/m².

juvenile_data_all_tidy_Apr_2026.csv contains the count-based density index used for the MARSS response. Additionally, it contains depletion-estimated abundance densities used for manuscript sensitivity comparison. nps_fish_pool_abundance_Apr_2026.csv and nps_fish_pool_reach_final.csv are the NPS components of those summaries.

Raw redd files

This dictionary covers data/raw/02_redds_data_raw/Olema_Redwood_spawners.csv and the nonspatial Lagunitas-subbasin attribute table LG_redds_CLEAN_UPDATED_2005_2022.csv.

Variable	Description	Units or values
`EventID`	NPS source-database identifier for the spawning-ground survey event.	Integer identifier.
`tblEvents_StreamID`	NPS source-database stream identifier. The analysis explicitly uses ID 5 for Olema Creek and ID 34 for Redwood Creek when cleaning the NPS redd data; only Olema is retained in the focal covariate. Other source codes are not used by the current analysis.	Integer code.
`SpawnReddMeasurementID`	Unique redd record identifier within the NPS survey data fromthe NPS source-database.	Integer identifier.
`SpawnerSeason`	Spawning season written as starting and ending calendar years. The ending year defines `WaterYear`.	Text such as `2004-2005`.
`ReddConfidence`	NPS source-database surveyor confidence in the field species identification for the redd. Values are blank, `Definite`, or `Not Definite`.	Categorical text.
`ReddSpecies`	NPS source-database field-assigned species code for the redd: `CO` = coho salmon, `CH` = Chinook salmon, `SH` = steelhead trout, and `UNK` = unknown.	Categorical code.
`PredictedValue_kNN`	NPS source-database species predicted to have constructed the redd by the k-nearest-neighbor classification: `CO`, `CH`, or `SH`; `na` means no prediction.	Categorical code.
`LocationID`	NPS source-database identifier for the stream location containing the redd.	Integer identifier.
`SectionNum`	NPS source-database monument-section number.	Integer identifier.
`Latitude`, `Longitude`	Generalized NPS redd coordinates rounded to 0.01 degrees. These fields are not used by script 02.	Decimal degrees at 0.01-degree precision.
`Date`	Survey or redd-observation date in the Lagunitas-subbasin attribute table.	Date.
`Redd_ID`	Unique Lagunitas-subbasin redd identifier. Its prefix is used to assign `DG`, `SG`, or `LG`.	Text identifier.
`Species`	Species attribute from the Lagunitas-subbasin redd source data.	Source categorical text or code.

The NPS redd export contained an unused CoordSystem = UTM field that contradicted the longitude and latitude values. That field was omitted from the deposited CSV; the current redd-processing script does not use coordinates.

Precise point geometry from the Lagunitas-subbasin redd source was not deposited. The retained CSV contains all 2,773 attribute records used to calculate stream-year redd counts.

Processed redd files

This dictionary covers redds_streamwide_tidy_2005_2022.csv and updated_all_redds_wide_2005_2022.csv in data/processed/02_redds_data_processed/.

Variable	Description	Units or values
`n`	Number of focal coho redds retained for the stream-water-year after the filtering and duplicate-removal rules in script 02.	Redd count.
Four-digit columns `2005`–`2022`	Streamwide value of `n` for the named water year.	Redd count.

Daily discharge and FFC upload files

This dictionary covers the three daily discharge files retained in data/raw/03_flow_data_raw/ and the three files in data/processed/03_flow_data_processed/ffc_upload/.

Variable	Description	Units or values
`Mean_daily_discharge`	Mean discharge for the named calendar date.	cfs.
`flow`	Mean daily discharge formatted for upload to the Functional Flows Calculator.	cfs.

sg_ffc_upload_reconstructed_1981_2022.csvContains observed San Geronimo flow plus donor-model predictions for the missing period from September 14 through December 31, 2021, as documented in script 03. The FFC requires the exact upload headers date and flow, with flow in cfs (FFC User's Guide).

FFC annual-result matrices

This dictionary covers the three files in data/processed/03_flow_data_processed/ffc_results/.

The column labeled Year contains the metric name, not a year. Each four-digit column is a water year and contains that row's metric value. Metric definitions and units follow the California Environmental Flows Framework Functional Flows Calculator documentation.

Row value in `Year`	Full meaning	Units
`Std`	Standard deviation of daily discharge within the water year.	cfs.
`Avg`	Average daily discharge within the water year.	cfs.
`CV`	Coefficient of variation of daily discharge, calculated as standard deviation divided by mean.	Unitless.
`FA_Mag`	Fall-pulse magnitude: maximum daily discharge during the identified fall-pulse event.	cfs.
`FA_Tim`	Fall-pulse timing: water-year day of the fall-pulse peak.	Water-year day.
`FA_Dur`	Fall-pulse duration.	Days.
`Wet_BFL_Mag_10`	Wet-season baseflow magnitude: 10th percentile of daily flows during the wet season, including peak-flow events.	cfs.
`Wet_BFL_Mag_50`	Wet-season median flow: 50th percentile of daily flows during the wet season, including peak-flow events.	cfs.
`Wet_Tim_Water`	Start of the wet season.	Water-year day.
`Wet_BFL_Dur`	Wet-season baseflow duration from wet-season start to spring-recession start.	Days.
`Peak_2`, `Peak_5`, `Peak_10`	Magnitudes of the 2-, 5-, and 10-year recurrence-interval peak flows.	cfs.
`Peak_Dur_2`, `Peak_Dur_5`, `Peak_Dur_10`	Cumulative number of days in which the corresponding recurrence-interval peak-flow magnitude is exceeded.	Days.
`Peak_Fre_2`, `Peak_Fre_5`, `Peak_Fre_10`	Frequency of the corresponding recurrence-interval peak-flow event within a season.	Occurrences.
`SP_Mag`	Spring-recession magnitude: daily flow at the start of the spring-flow period, four days after the final wet-season peak.	cfs.
`SP_Tim`	Start of the spring recession.	Water-year day.
`SP_Dur`	Duration from spring-recession start to dry-season baseflow start.	Days.
`SP_ROC`	Spring-recession rate of change: median daily rate of change over decreasing-flow periods during the recession.	Percent.
`DS_Mag_50`	Dry-season median baseflow: 50th percentile of daily flow during the dry season.	cfs.
`DS_Mag_90`	Dry-season high baseflow: 90th percentile of daily flow during the dry season.	cfs.
`DS_Tim`	Start of dry-season baseflow.	Water-year day.
`DS_Dur_WS`	Dry-season duration from dry-season start to the following wet-season start.	Days.

Long-term functional-flow reference data

This dictionary covers data/raw/03_flow_data_raw/ffm_longterm_all_streams_observed.csv. This source dataset supplies the long-term functional-flow values reported in manuscript Table 1. Those values were compiled manuall,y and the file is not read programmatically by the current scripts. It is retained to document and permit verification of the reported stream attributes. This file is distinct from the annual FFC result files used for the hydroclimatic PCA in script 03. Field meanings follow the accompanying California Natural Flows Database metadata.

Variable	Description	Units or values
`comid`	Unique NHD+ flowline identifier for the stream reach.	Integer identifier.
`ffm`	Functional-flow metric code. Codes use lowercase versions of the FFC metrics defined above, such as `ds_dur_ws`, `ds_mag_50`, `fa_mag`, `peak_fre_2`, `sp_roc`, and `wet_tim`.	Categorical code.
`wyt`	Water-year type used to summarize the metric.	`all`, `wet`, `moderate`, or `dry`.
`p10`, `p25`, `p50`, `p75`, `p90`	10th, 25th, 50th, 75th, and 90th percentiles of the metric across the contributing water years. These describe interannual variability rather than uncertainty in an individual prediction.	Units given in `unit`.
`unit`	Unit of the functional-flow metric.	`cfs`, `days`, `water year day`, `occurrences`, or `percent`.
`source`	Source or method for the metric values.	`observed`, `model`, or `inferred`.
`gage_id`	Stream-gage identifier where applicable.	Text identifier.
`observed_years`	Number of observed water years contributing to the summary, where applicable.	Year count.
`observed_year_start`, `observed_year_end`	First and last observed water years contributing to the summary.	Water year.
`alteration`	Source classification of likely flow alteration.	Blank, `indeterminate`, `likely_altered_high`, `likely_altered_low`, or `likely_unaltered`.

Hydroclimatic PCA file

This dictionary covers data/processed/03_flow_data_processed/flow_pc1_wide.csv.

Variable	Description	Units or values
`metric`	Name of the PCA score. The retained value i`flow_PC1`, the first principal component of eight annual FFC metrics.	Categorical text.
Four-digit columns `2005`–`2022`	Stream-specific PC1 score for the named water year. The contributing metrics were standardized within stream before the pooled PCA.	Unitless principal-component score.

Higher PC1 scores represent more front-loaded hydrographs, with earlier seasonal transition, faster recession, and longer dry-season duration. San Geronimo scores are used as the Devil's Gulch hydroclimatic proxy because Devil's Gulch lacks a long-term gauge.

Mapping and spatial attribute tables

`efish_locations_nps_2023.csv`

Variable	Description	Units or values
`EventID`, `LocationID`, `SectionNum`	NPS sampling-event, location, and monument-section identifiers, as defined above.	Integer identifiers.
`Latitude`, `Longitude`	Generalized site coordinates interpreted by script 06 as WGS 84 geographic coordinates (`EPSG:4326`) and rounded to 0.01 degrees.	Decimal degrees at 0.01-degree precision.

The NPS source export also contained unused CoordUnits = m and CoordSystem = UTM fields that contradicted the decimal-degree values. Those fields were omitted from the deposited CSV; no retained script uses them.

`juvsites_generalized.dbf`

This is the attribute table for the generalized Marin Water juvenile-site shapefile used in the study-system map. It contains only the eight sites mapped by script 06.

Variable	Description	Units or values
`SITE`	Juvenile sampling reach/site code.	Text identifier.

The generalized shapefile uses WGS 84 (EPSG:4326). Its point coordinates were rounded to 0.01 degrees before deposit; precise source geometry and unused source attributes were not retained.

`Geology.dbf`

The geology layer schema matches Marin County's published Geology GIS service and is derived from regional geologic mapping cited by the source. Script 07 transforms the geometry to California Albers (EPSG:3310) before calculating square-metre intersection areas; it uses AgeOrder as the geology class in the percentage summaries.

Variable	Description	Units or values
`OBJECTID`	Source GIS feature identifier.	Integer identifier.
`PType`	Geologic map-unit symbol or polygon-type code used in the source map legend, such as `Tpr`.	Categorical code.
`MapUnitDes`	Full source description of the geologic map unit, including lithology and commonly its geologic age.	Text.
`AgeOrder`	Source integer used to order or group geologic units by age. It is an ordinal/source classification, not an age measured in years.	Integer code.
`RGB_Render`	Red-green-blue color string used to render the source map unit. This is the shapefile-truncated form of the source field. `RGB_Rendering`.	Text color code.
`SHAPESTAre`	Source polygon area calculated in the layer's native State Plane coordinate system.	Square US survey feet.
`SHAPESTLen`	Source polygon perimeter calculated in the layer's native State Plane coordinate system.	US survey feet.

Authoritative sources are the Marin County Geology GIS layer and the USGS regional map database, Blake et al. (2000), doi:10.3133/mf2337.

MARSS model-selection and coefficient files

This dictionary covers:

final_candidate_models_AICc_2026-06-11.csv
table_S1_base_model_comparisons_2026-06-11.csv
table_S2_unequal_drift_estimates_2026-06-11.csv
model_coefficients_2026-06-11.csv
their manuscript-ready counterparts in outputs/tables/manuscript_2026-06-11/

Variable	Description	Units or values
`AICc`	Akaike information criterion corrected for small sample size. Lower values indicate greater relative support within the stated candidate set.	Unitless.
`delta_AICc`, `Delta AICc`	Difference between a model's AICc and the minimum AICc in the candidate set.	Unitless.
`rel_lik`	Relative likelihood, calculated as `exp(-0.5 × delta_AICc)`.	Unitless.
`weight`, `AICc weight`	AICc weight normalized to sum to 1 within the candidate set.	Proportion, 0–1.
`Covariates`	Biological or hydroclimatic covariates included in the model.	Categorical text.
`C structure`	Constraint structure imposed on the MARSS covariate-effect matrix.	Categorical text.
`term`	Internal MARSS parameter name. Prefixes indicate `C` covariate effects, `Q` process variance/covariance, `R` observation variance, `U` drift, or `x0` initial state.	Categorical text.
`parameter`	Human-readable interpretation`term`, including the affected stream, reach, or parameter class.	Categorical text.
`coefficient_class`	Broad MARSS parameter class: `C`, `Q_var`, `Q_cov`, `R`, `U`, or `x0`.	Categorical code.
`Drift estimate`	Estimated annual change in the latent `log(density + 1)` state for the stream under an unequal-drift model.	Log-scale change per year.

Both redd counts and hydroclimatic PC1 scores are z-scored before fitting the covariate models. Their C coefficient estimates therefore represent the change in latent log(density + 1) per one-standard-deviation increase in the covariate.

Latent-state and stream-summary files

This dictionary covers states_log_density_2026-06-11.csv, states_raw_density_2026-06-11.csv, table_2_base_state_summary_2026-06-11.csv, and Table_2_stream_parameters_2026-06-11.csv.

Variable	Description	Units or values
`raw`	Back-transformed latent-state estimate, calculated as `exp(estimate) - 1`.	Fish/m².
`raw_low`, `raw_high`	Back-transformed lower and upper 95% confidence bounds.	Fish/m².
`mean_density`, `Modeled mean density (fish/m2)`	Mean back-transformed latent density across 2005–2022 for the retained base model.	Fish/m².
`max_density`, `Modeled max density (fish/m2)`	Maximum back-transformed latent density across 2005–2022 for the retained base model.	Fish/m².
`sd_density`	Standard deviation of back-transformed latent density across years.	Fish/m².
`cv_percent`	`100 × sd_density / mean_density`.	Perce	Presented
	ented as a ratio in the manuscript table.	Unitless ratio.
`Process variance, Q (95% CI)`	Formatted process-error variance estimate and 95% confidence interval for the stream.	Squared `log(density + 1)` scale.

Process-correlation and Q-decomposition files

This dictionary covers process_correlations_2026-06-11.csv, all retained q_*.csv model outputs, table_S4_q_variance_correlation_change_2026-06-11.csv, the Figure 2 correlation CSVs, Figure_S2_process_synchrony_across_R_structures_2026-06-11.csv, and Table_S3_basin_level_Q_correlation_change_2026-06-11.csv.

Variable	Description	Units or values
`ID_1`, `ID_2`	Numeric indices of the paired latent stream states. Indices 1–4 correspond to `DG`, `LG`, `SG`, and `OL`.	Integer index.
`Pairwise_Value`	Pairwise correlation derived from the MARSS process covariance matrix.	Correlation, -1 to 1.
`Var1`, `Var2`	First and second stream in a pairwise correlation table.	Stream code.
`Freq`	Correlation value produced when a correlation matrix is converted to long format.	Correlation, -1 to 1.
`Pair`, `pair`	Human-readable paired-stream label.	Text identifier.
`correlation_source`	Type of correlation: observed density, latent modeled density, or MARSS process error.	Categorical text.
`pair_group`	Group assigned to a stream pair: within the Lagunitas subbasin or an Olema cross-basin pair.	Categorical text.
`n_pairs`	Number of stream pairs summarized in the pair group.	Integer count.
`DG_var`, `LG_var`, `SG_var`, `OL_var`	Process-error variance for the named stream under the specified model.	Squared `log(density + 1)` scale.
`DG_LG_corr`, `DG_SG_corr`, `LG_SG_corr`, `DG_OL_corr`, `LG_OL_corr`, `SG_OL_corr`	Process-error correlation for the named stream pair.	Correlation, -1 to 1.
`mean_abs_corr`	Mean absolute value of the six off-diagonal process correlations.	Unitless correlation magnitude.
`DG_var_change`, `LG_var_change`, `SG_var_change`, `OL_var_change`	Named stream's process variance minus its value in the base model.	Squared log-scale difference.
`mean_abs_corr_change`	Model mean absolute process correlation minus the base-model value.	Correlation difference.
`metric_type`	Indicates whether a row summarizes stream variance, an overall correlation measure, a pairwise correlation, or a pair-group mean.	Categorical text.
`metric`	Name of the summarized measure, such as `mean_abs_corr`.	Categorical text.
`base_var`, `redds_var`, `pc1_var`	Process variance from the base, redd-only, and redd-plus-PC1 models.	Squared `log(density + 1)` scale.
`base_value`, `redds_value`, `pc1_value`	Overall correlation metric from the base, redd-only, and redd-plus-PC1 models.	Unitless.
`base_corr`, `redds_corr`, `pc1_corr`	Pairwise process correlation from the base, redd-only, and redd-plus-PC1 models.	Correlation, -1 to 1.
`mean_base_corr`, `mean_redds_corr`, `mean_pc1_corr`	Mean pairwise process correlation within `pair_group` for the three model stages.	Correlation, -1 to 1.
`redds_abs_change`	Redd-model value minus base-model value.	Same scale as the underlying variance or correlation.
`redds_pct_explained`	`100 × (base - redd model) / base`. Positive values indicate reduction relative to the base model.	Percent.
`climate_abs_change`	Redd-plus-PC1 value minus redd-model value.	Same scale as the underlying variance or correlation.
`climate_pct_explained_from_base`	`100 × (redd model - redd-plus-PC1 model) / base`.	Percent.
`climate_pct_explained_from_redds`	`100 × (redd model - redd-plus-PC1 model) / redd model`.	Percent.
`total_pct_explained`	`100 × (base - redd-plus-PC1 model) / base`.	Percent.
`mean_redds_abs_change`, `mean_climate_abs_change`	Pair-group mean of the corresponding absolute-change field.	Correlation difference.
`mean_redds_pct_explained`, `mean_climate_pct_explained_from_base`, `mean_climate_pct_explained_from_redds`, `mean_total_pct_explained`	Pair-group means of the corresponding percentage-change fields.	Percent.

PVA threshold, risk, and buffering files

This dictionary covers all retained CSVs in outputs/model/pva/ and their manuscript-ready derivatives.

Variable	Description	Units or values
`q_scenario`	Process-covariance simulation scenario. `Full Q` retains estimated among-stream covariance; `Diagonal Q` retains the same stream variances but sets covariance to zero.	Categorical text.
`decline_level`	Proportional decline from the long-term mean raw density. Raw output values are 0.80, 0.90, or 0.95; some manuscript tables format them as `80%`, `90%`, or `95%`.	Proportion or formatted percent.
`threshold_label`, `Decline threshold`	Display label for `decline_level`.	Percentage text.
`start_log`	Stream's fixed PVA starting state, equal to its mean modeled latent state over 2005–2022.	`log(density + 1)`.
`start_raw`	Back-transformed starting density, `exp(start_log) - 1`.	Fish/m².
`threshold_raw`	Raw-density threshold, `start_raw × (1 - decline_level)`.	Fish/m².
`threshold_log`	`log(threshold_raw + 1)`, used for comparison with simulated or modeled latent states.	`log(density + 1)`.
`event`, `Event`	Joint event being evaluated: at least one, at least two, at least three, or all four streams below their thresholds in the same year.	`1+ streams`, `2+ streams`, `3+ streams`, or `4 streams`.
`annual_risk`, `Mean same-year risk`	Fraction of 10,000 simulations in which the event occurs in the specified future year. The manuscript-ready field is formatted as a percentage.	Proportion or percent.
`cumulative_risk`, `10-year cumulative risk`	Fraction of simulations in which the event has occurred at least once by the specified year. The manuscript-ready field is formatted as a percentage.	Proportion or percent.
`risk_without_stream`	Cumulative probability that all three nonfocal streams cross their thresholds.	Proportion, 0–1.
`risk_with_stream`	Cumulative probability that all four streams, including the focal stream, cross their thresholds.	Proportion, 0–1.
`buffering`	`risk_without_stream - risk_with_stream`; the reduction in joint-decline risk attributable to requiring the focal stream also to decline.	Proportion-point difference.
`crossed_threshold`	Whether the modeled latent state is at or below `threshold_log` for the stream-year.	Boolean `TRUE` or `FALSE`.
`n_years_crossed`	Number of observed study years in which the stream's modeled latent state crossed the threshold.	Year count.
`proportion_years_crossed`	`n_years_crossed` divided by the 18-year study length.	Proportion, 0–1.
`lowest_log_state`	Minimum modeled latent state observed for the stream during 2005–2022.	`log(density + 1)`.
`n_years`, `Years observed`	Number of study years in which the specified joint event occurred.	Year count.
`proportion_years`, `Proportion of time series`	Fraction of study years in which the specified joint event occurred; the manuscript-ready field is formatted as a percentage.	Proportion or percent.
`ever_observed`, `Observed at least once`	Whether the specified joint event occurred at least once during 2005–2022.	Boolean.

Remaining manuscript-ready table variables

Variable	Description	Units or values
`n` in `Figure_S1_raw_count_abundance_summary_2026-06-11.csv`	Number of paired reach-year observations in the count-density versus depletion-estimated-density comparison.	Observation count.
`r`	Pearson correlation between count density and depletion-estimated density.	Correlation, -1 to 1.
`r_squared`	Coefficient of determination from the linear regression of depletion-estimated density on count density.	Proportion, 0–1.
`slope`	Regression slope for depletion-estimated density as a function of count density.	Unitless because both variables are fish/m².
`intercept`	Regression intercept.	Fish/m².

Non-tabular data

Shapefiles consist of .shp, .shx, .dbf, .prj, and associated sidecars; all components should remain together.
prism_tmean_annual/ and prism_ppt_annual/ contain the annual 30-year normal mean-temperature and precipitation .bil rasters and their supporting sidecar files. The original PRISM basenames were shortened solely for Windows path compatibility. Script 07 reads the renamed rasters directly and identifies temperature as degrees Celsius and precipitation as millimetres.
nlcd_raster.tif, NLCDclasses.pdf, and NLCD_Colour_Classification_Update.jpg provide the National Land Cover Database (NLCD) land-cover raster and local class-reference files. See the USGS NLCD documentation and land-cover class descriptions for product and class details.

Spatial RDS files

An .rds file stores one R object. The spatial files below contain sf or sfc objects; therefore, software may report each attribute column as a separate component. The files with more than 100 components are spatial NHDPlus data frames, not model or PVA outputs. Their attributes follow the USGS NHDPlus documentation.

File	Main contents and variables	Units or expected values	Use
`data/raw/06_manuscript_figures_data_raw/flowline_watershed.rds`	`sf` data frame with 129 NHDPlus flowline features and 139 columns. Main fields include `comid`, `gnis_name`, `lengthkm`, stream order and network identifiers, drainage-area fields, smoothed elevation, slope, annual/monthly flow and velocity attributes, and `geometry`.	Length in km; drainage area in km²; NHDPlus elevation values in cm; slope unitless; flow in cfs; velocity in feet/s; NAD83 line geometry.	Watershed flowlines in the study-system map.
`data/raw/06_manuscript_figures_data_raw/flowline_watershed_subset_updated.rds` and `data/raw/07_stream_attributes_data_raw/flowline_watershed_subset_updated.rds`	Identical copies of an `sf` data frame with 19 selected flowlines and 140 columns. It contains the same NHDPlus fields`flowline_watershed.rds`, plus `Stream`.	Same NHDPlus units as above; `Stream` is a categorical stream name; NAD83 line geometry.	Focal stream map and stream-attribute calculations.
`data/raw/06_manuscript_figures_data_raw/total_catchment_watershed.rds`	One `sfc_MULTIPOLYGON` object representing the combined watershed/catchment boundary; it has no separate attribute columns.	NAD83 polygon geometry.	Study-system map boundary.
`data/raw/06_manuscript_figures_data_raw/waterbody_watershed.rds`	`sf`DataFrame with six waterbody features and 23 columns. Main fields include `comid`, `gnis_name`, `areasqkm`, feature codes, depth and volume attributes, lake area, and `geometry`.	Area in km² or m² as named by the source field; depth in m; volume in m³; NAD83 polygon geometry.	Waterbodies in the study-system map.
`data/raw/07_stream_attributes_data_raw/simplified_catchment.rds`	`sf` Data frame with 111 catchment polygons and 147 columns. Main fields include `comid`, `gnis_name`, catchment area, joined NHDPlus flowline/network attributes, an `area` field, and `geometry`.	Source drainage area in km²; `area` in m²; other NHDPlus units as above; NAD83 polygon geometry.	Catchment summaries in script 07.

MARSS model RDS files

The model files contain fitted marssMLE objects created by script 04. Standard components within each fit include model and marss model specifications; par and coef parameter estimates; states and states.se smoothed latent states and their standard errors; ytT and ytT.se expected reach observations and their standard errors; logLik, AIC, and AICc fit statistics; and estimation metadata such as control, method, numIter, and convergence. These components follow the MARSS package documentation.

File	Top-level objects and interpretation	Units or expected values	Use
`outputs/model/candidate_models_pc1_2026-06-11.rds`	Named list of nine fitted models: `base_region`, `base_trib`, `base_agency`, `base_stream`, `base_reach`, `redds`, `pc1_stream`, `pc1_region`, and `pc1_trib`.	`states` is a 4-stream × 18-year matrix on the `log(density + 1)` scale; `ytT` is a 10-reach × 18-year matrix on that scale; coefficients use the scales of their named parameters; likelihood and information criteria are unitless.	Model comparison, PVA input, and manuscript figures/tables.
`outputs/model/base_model_drift_comparison_2026-06-11.rds`	Named list of ten fitted models: five `unequal_u_` models and five `zero_u_` models under regional, tributary/mainstem, agency, stream, and reach observation-error structures.	Same fitted-object scales as above; `U` drift estimates are annual changes in latent `log(density + 1)`.	Drift and observation-error model comparisons and supporting results.

Source-data references

U.S. Geological Survey. National Hydrography Dataset Plus: NHDPlus High Resolution documentation and NHDPlus High Resolution user guide and data dictionary.
U.S. Geological Survey. National Land Cover Database (NLCD): NLCD documentation and product overview and land-cover class descriptions.
National Park Service, San Francisco Bay Area Network coho salmon monitoring: NPS monitoring documents.
California Environmental Flows Framework Technical Team. Functional Flows Calculator: tools and documentation and FFC User's Guide.
The Nature Conservancy. California Natural Flows Database: data and API.
PRISM Climate Group: PRISM climate data.
Marin County: Geology GIS layer.
Blake, M.C., Graymer, R.W., Jones, D.L., and Soule, A. 2000. Geologic map and map database of parts of Marin, San Francisco, Alameda, Contra Costa, and Sonoma counties, California. USGS Miscellaneous Field Studies Map MF-2337. doi:10.3133/mf2337.

Code and software

Analyses were conducted in R using R Markdown. The deposited session information reports:

R version 4.5.2 (2025-10-31 ucrt) -- "[Not] Part in a Rumble"
Platform: x86_64-w64-mingw32/x64

Principal R packages used across the scripts include tidyverse, here, lubridate, MARSS, MASS, sf, terra, raster, exactextractr, nhdplusTools, tmap, ggplot2, and related plotting packages loaded in the individual scripts.