diff --git a/README.md b/README.md index 9c9a414..85ee578 100644 --- a/README.md +++ b/README.md @@ -1,19 +1,34 @@ -# Single-Column ParFlow-CLM ET Simulation +# Single-Column ParFlow-CLM Examples -Run PF-CLM single-column simulations at any Ameriflux site using improved -CLM physics (Medlyn stomata, CLM5 tanh interception, canopy clumping). -Compare modeled latent heat against tower observations. +Two parallel worked examples of single-column PF-CLM simulations with improved +CLM physics, organized by topic: + +- **`et/`** — ET example (Medlyn stomata, CLM5 tanh interception, canopy clumping) + compared against Ameriflux flux-tower observations. +- **`snow/`** — Snow example (fractional snow cover, albedo decay, shrubland land + cover) compared against SNOTEL SWE observations. + +Each topic is a three-notebook workflow that shares site-configuration and +forcing infrastructure through the top-level `helpers.py`. ## Prerequisites -### ParFlow (v3.14.1+) +### ParFlow (master) + +Both examples run on current `parflow/master`. The relevant CLM-physics PRs +are merged: -Requires ParFlow with CLM ET improvements from the `feature/clm_et` branch -(Medlyn stomata, CLM5 tanh interception, canopy clumping support). +| PR | Merged | What it adds | +|----|--------|--------------| +| #695 | 2026-01-28 | CLM snow parameterization options | +| #698 | 2026-02-02 | Extended CLM snow parameterization | +| #701 | 2026-02-11 | Additional CLM snow formulations | +| #709 | 2026-03-11 | SZA-modulated fractional snow cover | +| #712 | 2026-03-19 | CLM ET improvements (Medlyn, CLM5 interception, PFT photosynthesis) | Build from source: https://github.com/parflow/parflow -Set `PARFLOW_DIR` to your installation: +Set `PARFLOW_DIR` to your install: ```bash export PARFLOW_DIR=/path/to/parflow/install ``` @@ -39,133 +54,217 @@ import hf_hydrodata as hf hf.register_api_pin("your_email", "your_pin") ``` -## Quick Start +## Repository Layout -1. **`locate_plot_station_get_forcing.ipynb`** -- Pick a site, query CONUS2 subsurface (soil + geology), set water table depth, download CW3E forcing, prepare CLM input files. Writes `site_config.txt` for Notebook 2. +``` +single_column_PFCLM/ +├── README.md # this file +├── helpers.py # shared site-config, CW3E, IGBP/soil tables, +│ # SNOTEL catalog, vegm builder +├── clm_inputs/ +│ ├── drv_clmin_template.dat # shared CLM driver template (both topics) +│ ├── et/ +│ │ └── drv_vegp.dat # clm_et_improved defaults +│ └── snow/ +│ └── drv_vegp.dat # clm_snow_improved defaults +├── et/ # ET example (3 notebooks) +│ ├── locate_plot_station_get_forcing.ipynb +│ ├── Single_Column_PFCLM_netcdf.ipynb +│ └── pfclm_ameriflux_compare.ipynb +├── snow/ # Snow example (3 notebooks) +│ ├── locate_plot_station_get_forcing.ipynb +│ ├── Single_Column_PFCLM_netcdf.ipynb +│ └── pfclm_snotel_compare.ipynb +└── runs/ # run output, subfoldered by topic + ├── et/ + └── snow/ +``` + +## Quick Start -2. **`Single_Column_PFCLM_netcdf.ipynb`** -- Load site config, configure two-layer subsurface (soil over geology), run ParFlow-CLM with improved physics, verify output. +Pick a topic, then run the three notebooks in order from that topic's folder. -3. **`pfclm_ameriflux_compare.ipynb`** -- Load output, fetch Ameriflux observations, compute metrics (KGE, NSE, RMSE, bias), plot time series, scatter, ET components, soil moisture, and 5-variable forcing quality comparison. +### ET example (`et/`) -## Subsurface Structure +1. **`locate_plot_station_get_forcing.ipynb`** — Pick an Ameriflux site, query + CONUS2 subsurface + water-table depth, download CW3E forcing, write CLM + inputs and `site_config.txt`. +2. **`Single_Column_PFCLM_netcdf.ipynb`** — Configure and run PF-CLM with the + `clm_et_improved` physics defaults. +3. **`pfclm_ameriflux_compare.ipynb`** — Compare modeled latent heat against + Ameriflux observations (KGE, NSE, RMSE, bias, time series, ET components). -The single-column domain is 8m deep with 10 computational layers and two subsurface zones: +### Snow example (`snow/`) -| Zone | Layers | Physical depth | Source | -|------|--------|---------------|--------| -| Soil | 6-9 (top) | 2.0m (1.0 + 0.6 + 0.3 + 0.1m) | CONUS2.1 top-layer soil type | -| Geology | 0-5 (bottom) | 6.0m (6 x 1.0m) | CONUS2.1 layer-4 geology type | +1. **`locate_plot_station_get_forcing.ipynb`** — Pick a SNOTEL site, download + CW3E forcing, preview observed SWE, write CLM inputs and `site_config.txt`. +2. **`Single_Column_PFCLM_netcdf.ipynb`** — Configure and run PF-CLM with the + `clm_snow_improved` physics defaults (shrubland + slow albedo decay + high + fresh albedo). An SZA-scheme variant is commented in the same cell. +3. **`pfclm_snotel_compare.ipynb`** — Compare modeled SWE against SNOTEL + observations (NSE, RMSE, peak |bias|%, melt-out day Δ, time series, scatter, + CW3E-vs-SNOTEL precipitation check). -Soil hydraulic properties (Ksat, porosity, VG alpha/n) come from the CONUS2.1 soil -texture table (13 types). Geology uses per-type Ksat and porosity with domain-default -Van Genuchten parameters (alpha=0.5, n=2.5). +## Subsurface Structure (both topics) -Water table depth is queried from CONUS2 baseline and Ma et al. (2025) 30m products, -set by the user, and applied as an equilibrium bottom boundary condition. +Both examples use the same 8m, 10-layer single-column domain with two zones: -## Tested Sites +| Zone | Layers | Physical depth | Source | +|---------|---------|--------------------------------|---------------------------| +| Soil | 6-9 (top) | 2.0 m (1.0 + 0.6 + 0.3 + 0.1 m) | CONUS2.1 top-layer soil type | +| Geology | 0-5 (bot) | 6.0 m (6 × 1.0 m) | CONUS2.1 layer-4 geology type | -| Site ID | Name | IGBP | State | WY | KGE | -|---------|------|------|-------|----|-----| -| US-Slt | Silas Little | DBF | NJ | 2012 | 0.48 | -| US-xUK | NEON Univ. of Kansas | DBF | KS | 2024 | 0.45 | -| US-xBL | NEON Blandy Farm | DBF | VA | 2024 | 0.42 | -| US-GLE | GLEES | ENF | WY | 2020 | 0.15 | -| US-Ro5 | Rosemount I18 South | CRO | MN | 2024 | 0.68 | -| US-Ro1 | Rosemount G21 | CRO | MN | 2016 | 0.55 | -| US-Ro2 | Rosemount C7 | CRO | MN | 2016 | 0.61 | -| US-Fwf | Flagstaff Wildfire | GRA | AZ | 2008 | 0.47 | -| US-Mpj | Mountainair Pinyon-Juniper | WSA | NM | 2024 | 0.35 | -| US-Ton | Tonzi Ranch | WSA | CA | 2024 | 0.30 | -| US-xSR | NEON Santa Rita | OSH | AZ | 2024 | 0.40 | -| US-xJR | NEON Jornada | OSH | NM | 2024 | 0.38 | -| US-xDL | NEON Dead Lake | MF | AL | 2024 | 0.52 | +Soil hydraulic properties (Ksat, porosity, VG α/n) come from the CONUS2.1 soil +texture table (13 types). Geology uses per-type Ksat + porosity with domain- +default Van Genuchten parameters (α=0.5, n=2.5). -KGE values are from the `phase8_clump` configuration (best overall). +Water table depth is queried from CONUS2 baseline or Ma et al. (2025) 30m +products in the ET notebook; the snow notebook defaults to a deep (-10 m) BC +since mountain SNOTEL sites are snowmelt-driven, not groundwater-fed. ## Physics Configuration -The shipped CLM input files implement the `phase8_clump` configuration: +### `clm_et_improved` (shipped in `clm_inputs/et/`) + +- **Medlyn stomatal conductance** (`StomataScheme="Medlyn"`): VPD-based stomata + replacing Ball-Berry. PFT-dependent g1 parameters in `drv_vegp.dat`. +- **CLM5 tanh canopy interception** (`InterceptionScheme="CLM5Tanh"`): + `fpi = tanh(LAI+SAI)`, replacing the CLM3 `0.25·(1-exp(-0.5·LAI))` formula. +- **Canopy clumping** (He et al. 2012): PFT-dependent clumping indices in + `drv_vegp.dat`, reducing effective LAI for radiation transfer in forest + canopies. +- **CLM4.5 optical/structural parameters**: Updated leaf/stem reflectance and + transmittance. +- **PFT-dependent Vcmax**: Custom photosynthesis with PFT-specific `vcmx25`. + +### `clm_snow_improved` (shipped in `clm_inputs/snow/`) + +- **Fractional snow cover**: `FracSnoScheme="CLM"`, `FracSnoRoughness=1e-8`. + This was the foundation fix in the sensitivity study — reduced peak + |bias| from 24% to 13%. +- **Albedo decay (slow)**: `AlbedoDecayVis=0.3`, `AlbedoDecayNir=0.1` — slows + snow aging so fresh-snow albedo persists longer into spring. +- **Fresh-snow albedo (tuned high)**: `AlbedoVisNew=0.97`, `AlbedoNirNew=0.67`. +- **Shrubland land cover** (IGBP 7): forces a snow-friendly canopy radiative + transfer treatment — best-performing IGBP class in the low-bias study. +- **Optional SZA-modulated fSCA** (commented in Notebook 2): the interpolating + FracSnoScheme="SZA" formulation from PR #709, with `FracSnoRoughnessMin/Max`, + `FracSnoGammaSZA`, and `FracSnoAvgWindow`. + +## Tested Sites — ET (Ameriflux) + +These are the Ameriflux sites from the ET sensitivity study. KGE values are +for the `clm_et_improved` (phase8_clump) configuration. + +| Site ID | Name | IGBP | State | WY | KGE | +|---------|-----------------------------|------|-------|------|------| +| US-Slt | Silas Little | DBF | NJ | 2012 | 0.48 | +| US-xUK | NEON Univ. of Kansas | DBF | KS | 2024 | 0.45 | +| US-xBL | NEON Blandy Farm | DBF | VA | 2024 | 0.42 | +| US-GLE | GLEES | ENF | WY | 2020 | 0.15 | +| US-Ro5 | Rosemount I18 South | CRO | MN | 2024 | 0.68 | +| US-Ro1 | Rosemount G21 | CRO | MN | 2016 | 0.55 | +| US-Ro2 | Rosemount C7 | CRO | MN | 2016 | 0.61 | +| US-Fwf | Flagstaff Wildfire | GRA | AZ | 2008 | 0.47 | +| US-Mpj | Mountainair Pinyon-Juniper | WSA | NM | 2024 | 0.35 | +| US-Ton | Tonzi Ranch | WSA | CA | 2024 | 0.30 | +| US-xSR | NEON Santa Rita | OSH | AZ | 2024 | 0.40 | +| US-xJR | NEON Jornada | OSH | NM | 2024 | 0.38 | +| US-xDL | NEON Dead Lake | MF | AL | 2024 | 0.52 | + +## Showcase Sites — Snow (SNOTEL) + +The snow notebook ships a 6-site showcase picked for geographic diversity +and recent CW3E-covered water years: + +| Site | Region | State | Elev (m) | WY | Triplet | +|--------------------|-------------------|-------|----------|------|--------------| +| Paradise | Cascades | WA | 1564 | 2022 | 679:WA:SNTL | +| CSS_Lab | Sierra Nevada | CA | 2103 | 2022 | 428:CA:SNTL | +| Togwotee_Pass | Wyoming | WY | 2936 | 2020 | 822:WY:SNTL | +| Berthoud_Summit | Colorado Rockies | CO | 3536 | 2024 | 335:CO:SNTL | +| Brighton | Wasatch | UT | 2667 | 2024 | 366:UT:SNTL | +| Quemazon | Southern Rockies | NM | 2926 | 2023 | 708:NM:SNTL | + +### Extending to more sites + +The `snow_model_sensitivity` study identified 20 **low-bias site-years** +(|T bias| < 1°C and |P bias| < 10% at the SNOTEL gauge) that are reliable +for PF-CLM evaluation: + +| Region | Sites | Years | +|---------------------|------------------------------------------------|----------------------| +| Cascades | Stevens_Pass, Paradise | 2020, 2022, 2023 | +| Sierra Nevada | Leavitt_Lake, CSS_Lab, Donner_Summit | 2015, 2021, 2022, 2024 | +| Wyoming | Togwotee_Pass, Canyon | 1995, 2000, 2020, 2024 | +| Colorado Rockies | Berthoud_Summit, Schofield_Pass | 2024 | +| Northern Rockies | Northeast_Entrance | 1995 | +| Wasatch | Brighton | 2024 | +| Southern Rockies | Quemazon | 2020, 2023 | + +To use another site: +1. Add a new entry to `SNOTEL_SITES` in `helpers.py` (coords, elev_m, state, + region, triplet). Triplets can be looked up at + https://www.nrcs.usda.gov/wps/portal/wcc/home/. +2. Set `site_name` and `water_year` in Notebook 1. +3. The notebooks pick up the rest automatically. -- **Medlyn stomatal conductance** (`StomataScheme = "Medlyn"`): VPD-based stomata replacing Ball-Berry. PFT-dependent g1 parameters in `drv_vegp.dat`. -- **CLM5 tanh canopy interception** (`InterceptionScheme = "CLM5Tanh"`): `fpi = tanh(LAI+SAI)`, replacing the CLM3 `0.25*(1-exp(-0.5*LAI))` formula. -- **Canopy clumping** (He et al. 2012): PFT-dependent clumping indices in `drv_vegp.dat`, reducing effective LAI for radiation transfer in forest canopies. -- **CLM4.5 optical/structural parameters**: Updated leaf/stem reflectance and transmittance. -- **PFT-dependent Vcmax**: Custom photosynthesis with PFT-specific `vcmx25` values. +## Choosing a Site and Water Year -CLM input files use headers from the `feature/clm_driver_cleanup` branch. +**For ET:** -## Choosing a Site and Water Year +- **Data coverage**: use the Bokeh preview in the ET Notebook 1 to check + Ameriflux latent-heat coverage. +- **IGBP type**: CRO sites tend to perform best (KGE 0.55-0.68). Forest sites + (DBF, ENF, MF) do well with canopy clumping. WSA/OSH are more WTD-sensitive. +- **Water year**: choose a year with complete CW3E coverage (v1.0 covers + WY2003-present). + +**For snow:** + +- **Pair elevation with air temperature**: shallow or intermittent-snow sites + (Pacific NW low-elevation, arid SW) are harder to fit than continental + high-elevation sites. +- **Forcing bias check**: Notebook 3 plots cumulative CW3E vs. SNOTEL precip. + If CW3E/SNOTEL > 1.15 or < 0.85, expect model SWE to be systematically off. +- **Peak SWE year**: the showcase years were chosen from the low-bias set. + Run them first before experimenting with more difficult site-years. + +## Water Table Depth (ET only) -Any Ameriflux site within the CONUS2 domain can be used. When selecting a site and water year: - -- **Data coverage**: Check that Ameriflux has latent heat data for your chosen water year - using the Bokeh preview in Notebook 1. Gaps in tower data reduce the value of the comparison. -- **IGBP type**: Cropland (CRO) sites tend to perform best (KGE 0.55-0.68). Forest sites - (DBF, ENF, MF) perform well with canopy clumping. Arid/semi-arid sites (WSA, OSH) are - more sensitive to WTD and soil parameters. -- **Water year selection**: Choose a year with complete forcing coverage (CW3E v1.0 covers - WY2003-present). Avoid years with known instrument issues at the tower site. - -## Water Table Depth - -Water table depth (WTD) is the most sensitive parameter for ET in single-column simulations. -It controls how much groundwater is available to sustain transpiration during dry periods. - -**Sources available in Notebook 1:** -- **CONUS2 baseline**: Monthly mean WTD from the CONUS2.1 simulation. Not available for - all water years. -- **Ma et al. (2025) 30m**: Static high-resolution WTD estimate. Available everywhere but - represents a long-term mean, not a specific year. - -**Guidance:** -- Shallow WTD sites (coastal plains, wetlands, river valleys): use CONUS2 baseline if - available, or literature values. US-Slt (Pine Barrens, NJ) has WTD ~0.2m. -- Deep WTD sites (mountains, arid regions): WTD > 8m means the water table is below the - model domain. ET will be entirely rainfall-dependent with no groundwater contribution. - This is physically correct for many western US sites. -- The Ma 2025 product can overestimate WTD in areas with shallow water tables. Cross-check - with site literature or USGS well data when possible. -- WTD can be overridden directly in Notebook 2 without re-running Notebook 1. +WTD is the most sensitive parameter for ET in single-column simulations. +The snow notebook defaults to a deep WTD (-10 m, below the 8m domain) since +SNOTEL ridges are snowmelt-driven, not groundwater-fed. See the ET Notebook 1 +for full WTD sourcing options. ## Troubleshooting **Solver failure (run stops before completing the water year):** -- Check the ParFlow log (`pfclm_sc.out.log`) for timestep cutting. If dt drops below - ~0.01 and the run stalls, the solver cannot converge. -- Common cause: sharp property contrast at the soil/geology interface combined with - a wetting front or freeze/thaw event. Try adjusting WTD or increasing +- Check the ParFlow log (`pfclm_sc.out.log`) for timestep cutting. If dt + drops below ~0.01 and the run stalls, the solver cannot converge. +- Common cause: sharp property contrast at the soil/geology interface combined + with a wetting front or freeze/thaw event. Try adjusting WTD or increasing `Solver.Nonlinear.MaxIter`. -- Verify geology VG parameters are reasonable (n should be >= 2.0). VG n < 1.0 produces - unphysical saturation values. **HydroData errors:** -- `register_api_pin` only needs to run once per machine. If it fails, check your - credentials at https://hydrogen.princeton.edu/pin. -- Some Ameriflux variables (e.g., VPD, wind) are not available at all sites. The forcing - comparison cell skips unavailable variables automatically. -- If `get_gridded_data` fails for CONUS2 baseline WTD, the Ma 2025 product is used as - fallback. Both can be unavailable for some site/year combinations. +- `register_api_pin` only needs to run once per machine. If it fails, check + your credentials at https://hydrogen.princeton.edu/pin. +- Some Ameriflux variables (e.g., VPD, wind) are not available at all sites. +- SNOTEL triplets: if a fetch returns empty, verify the triplet at + https://www.nrcs.usda.gov/wps/portal/wcc/home/ and update `helpers.SNOTEL_SITES`. **Stale kernel after file edits:** -- If `helpers.py` is updated but the notebook throws `ImportError`, restart the kernel. - Jupyter caches imported modules and won't pick up on-disk changes without a restart. - -## File Structure - -``` -clm_inputs/ - drv_vegp.dat # PFT parameters (CLM4.5 + clumping, best-practice) - drv_clmin_template.dat # CLM driver template (30m forcing heights, dewmx=0.2) -helpers.py # vegm builder, clmin patcher, IGBP/soil/geology tables, - # site config read/write -locate_plot_station_get_forcing.ipynb -Single_Column_PFCLM_netcdf.ipynb -pfclm_ameriflux_compare.ipynb -``` +- If `helpers.py` is updated but a notebook throws `ImportError`, restart + the kernel. Jupyter caches imports and won't pick up on-disk changes + without a restart. ## References -- Maxwell, R.M. & Miller, N.L. (2005). Development of a coupled land surface and groundwater model. J. Hydrometeorol. -- Medlyn, B.E. et al. (2011). Reconciling the optimal and empirical approaches to modelling stomatal conductance. Global Change Biology. -- He, L. et al. (2012). Global clumping index map derived from the MODIS BRDF product. Remote Sensing of Environment. +- Maxwell, R.M. & Miller, N.L. (2005). Development of a coupled land surface + and groundwater model. J. Hydrometeorol. +- Medlyn, B.E. et al. (2011). Reconciling the optimal and empirical approaches + to modelling stomatal conductance. Global Change Biology. +- He, L. et al. (2012). Global clumping index map derived from the MODIS BRDF + product. Remote Sensing of Environment. +- Abolafia-Rosenzweig, R. et al. (2022). Evaluation of snow albedo + parameterizations in the CLM. diff --git a/clm_inputs/README.md b/clm_inputs/README.md new file mode 100644 index 0000000..ad78600 --- /dev/null +++ b/clm_inputs/README.md @@ -0,0 +1,97 @@ +# CLM Input Files + +Shipped CLM driver files for both the ET and snow examples. These track the +"best-practice" files published in +[parflow/parflow PR #714](https://github.com/parflow/parflow/pull/714) +(`feature/clm_driver_cleanup` → master). + +## Current Layout + +``` +clm_inputs/ +├── drv_clmin_template.dat # shared (ET + snow) — exact copy of +│ # test/tcl/clm/drv_clmin_bestpractice.dat +├── et/ +│ └── drv_vegp.dat # copy of test/tcl/clm/drv_vegp_bestpractice.dat +└── snow/ + └── drv_vegp.dat # currently identical to et/drv_vegp.dat +``` + +Each notebook copies the topic-specific `drv_vegp.dat` into its run directory +and patches `drv_clmin_template.dat`'s date fields for the selected water year. + +## Why the two `drv_vegp.dat` files are identical + +The PR #714 best-practice `drv_vegp.dat` is topic-agnostic — it carries the +full IGBP PFT table with all improvements needed by either example: + +- CLM4.5 optical/structural corrections (grass, savanna, crop) +- PFT-specific photosynthesis (`vcmx25` + C3/C4 fixes) +- Medlyn stomatal conductance (`g1_medlyn`) +- Canopy clumping index (`clump`, He et al. 2012) +- Foliage nitrogen (`folnmx`) + +None of those hurt a snow run — the snow-specific tuning lives outside this +file, in the ParFlow run script (`run.Solver.CLM.*` keys) and in the +site-specific `drv_vegm.dat` (IGBP=7 for the shrubland snow default). + +So today, `et/drv_vegp.dat` and `snow/drv_vegp.dat` are byte-identical copies +of the same upstream file. The duplication is intentional — it keeps the +door open for topic-specific tuning without a restructure. + +## Planned Consolidation + +When we're confident neither topic needs per-topic vegp overrides, collapse to: + +``` +clm_inputs/ +├── drv_clmin_template.dat # unchanged +└── drv_vegp.dat # single shared best-practice file +``` + +And update the two topic notebooks: + +| File | Change | +|------|--------| +| `et/locate_plot_station_get_forcing.ipynb` | `Path("../clm_inputs/et")` → `Path("../clm_inputs")` | +| `snow/locate_plot_station_get_forcing.ipynb` | `Path("../clm_inputs/snow")` → `Path("../clm_inputs")` | + +That's a three-line change (two notebook edits + delete the two subdirs). + +### Triggers that would *prevent* consolidation + +Reasons we might end up keeping the split: + +1. **Snow-specific PFT tuning**. If future work finds that snow-site canopy + radiation/transpiration benefits from different `z0m`, `displa`, `rhol_*`, + `taul_*` values at high-elevation shrubland sites, those would live in + `snow/drv_vegp.dat` and diverge from the ET file. +2. **Different `vcmx25` regimes for snow vs. ET sites**. Unlikely but + possible if PFT photosynthesis is retuned for cold/dormant-season behavior. +3. **Experimental physics branches**. If a future branch ships a new + `drv_vegp.dat` format with extra columns for one topic only. + +None of these apply today. + +## Upstream Tracking + +When `parflow/master` updates either best-practice file, refresh both copies +here. The fastest check is a `diff` against the upstream raw files: + +```bash +for f in drv_clmin_bestpractice drv_vegp_bestpractice; do + curl -sL "https://raw.githubusercontent.com/parflow/parflow/master/test/tcl/clm/${f}.dat" | \ + diff - clm_inputs/${f%_bestpractice}_template.dat # or et/drv_vegp.dat, snow/drv_vegp.dat +done +``` + +(Adjust target paths — `drv_clmin_bestpractice` lands in +`clm_inputs/drv_clmin_template.dat`, `drv_vegp_bestpractice` lands in both +`clm_inputs/et/drv_vegp.dat` and `clm_inputs/snow/drv_vegp.dat`.) + +## Background + +- PR #714: "Document CLM driver files and remove dead parameters" + (feature/clm_driver_cleanup, merged March 2026). +- Related PRs: #712 (CLM ET improvements), #695/#698/#701/#709 (snow). +- See the top-level `README.md` for the complete merged-PR table. diff --git a/clm_inputs/drv_vegp.dat b/clm_inputs/et/drv_vegp.dat similarity index 100% rename from clm_inputs/drv_vegp.dat rename to clm_inputs/et/drv_vegp.dat diff --git a/clm_inputs/snow/drv_vegp.dat b/clm_inputs/snow/drv_vegp.dat new file mode 100644 index 0000000..4952766 --- /dev/null +++ b/clm_inputs/snow/drv_vegp.dat @@ -0,0 +1,176 @@ +!========================================================================= +! ParFlow-CLM vegetation parameter file (IGBP classification) +! +! Part of ParFlow v3.14 (March 2026) +! https://parflow.org | https://github.com/parflow/parflow +! Originally based on CLM 1.0 (Dai et al. 2003) +! +! References: +! Maxwell & Miller (2005) J. Hydrometeorol. — ParFlow-CLM coupling +! Kollet & Maxwell (2008) Water Resour. Res. — Integrated watershed model +! Maxwell & Condon (2016) Science, SI — Coupling approach, parameter tables +! Kuffour et al. (2020) Geosci. Model Dev. — ParFlow v3.5 description +!========================================================================= +! drv_vegp_bestpractice.dat: +! +! DESCRIPTION: +! Best-practice vegetation parameters with CLM4.5 optical/structural +! corrections, PFT-dependent photosynthesis, and canopy clumping. +! +! Changes from default drv_vegp.dat: +! - CLM4.5 optical corrections for grass/savanna/crop (IGBP 8-10, 12): +! taus_vis, taus_nir → 0.001 (CLM4.5 Table 3.1) +! rhos_vis → 0.16, rhos_nir → 0.39 +! rhol_nir: IGBP 9,10,12: 0.58 → 0.35 +! - CLM4.5 structural: IGBP 10 sai 4.0→0.5, roota 1.0→11.0 (Zeng 2001) +! - C3/C4 photosynthesis parameter fixes +! - Medlyn stomatal conductance (g1_medlyn) +! - Foliage nitrogen (folnmx) +! - Canopy clumping index (clump) for radiation calculations +! +! Presence of vcmx25 keyword activates PFT-specific photosynthesis. +! +! Parameter sources: +! Oleson et al. (2013) CLM4.5 Tech Note — Optical/structural corrections +! Jefferson et al. (2017) J. Hydrometeorol. — Photosynthesis sensitivity +! Medlyn et al. (2011), Lin et al. (2015) — Stomatal conductance +! He et al. (2012) / Lawrence et al. (2019) — Canopy clumping index +! +! REVISION HISTORY: +! 1998-1999: Yongjiu Dai, Xubin Zeng; Original CLM code (NCAR) +! 2004-2009: Reed Maxwell, Stefan Kollet; ParFlow-CLM coupling +! 2009: Ian Ferguson; Irrigation schemes +! 2014: Nick Engdahl; Stomatal resistance multiplier, output controls +! 2016: Reed Maxwell, Bryant Reyes; Under-canopy drag (csoilc) correction +! 2017: Jennifer Jefferson; ET and transpiration improvements +! 2018-2019: Lindsay Bearup, Anna Ryken; Snow parameterizations +! 2022: Reed Maxwell; CLM update +! Jan-Mar 2026: Reed Maxwell; Snow parameterizations, ET improvements, +! PFT photosynthesis, Medlyn stomata, canopy clumping, +! CLM4.5 optical/structural corrections, driver cleanup +!========================================================================= +!IGBP Land Cover Types +! 1 evergreen needleleaf forests +! 2 evergreen broadleaf forests +! 3 deciduous needleleaf forests +! 4 deciduous broadleaf forests +! 5 mixed forests +! 6 closed shrublands +! 7 open shrublands +! 8 woody savannas +! 9 savannas +! 10 grasslands +! 11 permanent wetlands +! 12 croplands +! 13 urban and built-up lands +! 14 cropland / natural vegetation mosaics +! 15 snow and ice +! 16 barren or sparsely vegetated +! 17 water bodies +! 18 bare soil +!========================================================================= +! +itypwat (1-soil, 2-land ice, 3-deep lake, 4-shallow lake, 5-wetland: swamp, marsh) +1 1 1 1 1 1 1 1 1 1 5 1 1 1 2 1 3 1 +! +lai Maximum leaf area index [-] +6.00 6.00 6.00 6.00 6.00 6.00 6.00 6.00 6.00 2.00 6.00 6.00 5.00 6.00 0.00 6.00 0.00 0.00 +! +lai0 Minimum leaf area index [-] +5.00 5.00 1.00 1.00 3.00 2.00 1.00 2.00 1.00 0.50 0.50 0.50 1.00 2.00 0.00 0.50 0.00 0.00 +! +sai Stem area index [-] +2.00 2.00 2.00 2.00 2.00 2.00 2.00 2.00 2.00 0.50 2.00 0.50 2.00 2.00 2.00 2.00 2.00 0.00 +! CLM4.5: IGBP 10 sai 4.0→0.5 +! +z0m Aerodynamic roughness length [m] +1.00 2.20 1.00 0.80 0.80 0.10 0.10 0.70 0.10 0.03 0.03 0.06 0.50 0.06 0.01 0.05 0.002 0.01 +! +displa Displacement height [m] +11.0 23.00 11.0 13.0 13.0 0.30 0.30 6.50 0.70 0.30 0.30 0.30 3.00 0.30 0.00 0.10 0.00 0.00 +! +dleaf Leaf dimension [m] +0.04 0.04 0.04 0.04 0.04 0.04 0.04 0.04 0.04 0.04 0.04 0.04 0.04 0.04 0.04 0.04 0.04 0.00 +! +roota Fitted numerical index of rooting distribution +7.00 7.00 7.00 6.00 5.00 6.00 5.00 6.00 5.00 11.00 6.00 6.00 5.00 5.00 0.00 5.00 0.00 0.00 +! CLM4.5: IGBP 10 roota 1.0→11.0 (Zeng 2001) +! +rootb Fitted numerical index of rooting distribution +2.00 1.00 2.00 2.00 1.50 1.50 2.50 2.50 1.00 2.50 2.00 2.50 2.00 2.00 0.00 2.00 0.00 0.00 +! +rhol_vis !leaf reflectance vis +0.07 0.10 0.07 0.10 0.08 0.08 0.08 0.09 0.11 0.11 0.11 0.11 0.09 0.09 -99. 0.09 -99. -99. +! +rhol_nir !leaf reflectance nir +0.35 0.45 0.35 0.45 0.40 0.40 0.40 0.49 0.35 0.35 0.35 0.35 0.47 0.47 -99. 0.47 -99. -99. +! CLM4.5: IGBP 9,10,12 rhol_nir 0.58→0.35 +! +rhos_vis !stem reflectance vis +0.16 0.16 0.16 0.16 0.16 0.16 0.16 0.16 0.16 0.16 0.36 0.16 0.24 0.24 -99. 0.24 -99. -99. +! CLM4.5: IGBP 8,9,10,12 rhos_vis →0.16 +! +rhos_nir !stem reflectance nir +0.39 0.39 0.39 0.39 0.39 0.39 0.39 0.39 0.39 0.39 0.39 0.39 0.47 0.47 -99. 0.47 -99. -99. +! CLM4.5: IGBP 8,9,10,12 rhos_nir →0.39 +! +taul_vis !leaf transmittance vis +0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.06 0.07 0.07 0.07 0.07 0.06 0.06 -99. 0.06 -99. -99. +! +taul_nir !leaf transmittance nir +0.10 0.25 0.10 0.25 0.17 0.17 0.17 0.21 0.25 0.25 0.10 0.25 0.20 0.20 -99. 0.20 -99. -99. +! +taus_vis !stem transmittance vis +0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.22 0.001 0.09 0.09 -99. 0.09 -99. -99. +! CLM4.5: IGBP 8,9,10,12 taus_vis →0.001 (Table 3.1) +! +taus_nir !stem transmittance nir +0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.15 0.15 -99. 0.15 -99. -99. +! CLM4.5: IGBP 8,9,10,12 taus_nir →0.001 (Table 3.1) +! +xl !leaf/stem orientation index +0.01 0.10 0.01 0.25 0.13 0.13 0.13 -0.08 -0.3 -0.3 -0.3 -0.3 -0.07 -0.07 -99. -0.07 -99. -99. +! +vw !btran exponent:[(h2osoi_vol-watdry)/(watopt-watdry)]**vw +1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. -99. 1. -99. -99. +! +irrig !(irrig=0 -> no irrigation, irrig=1 -> irrigate) +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +! +! --- PFT-Specific Photosynthesis Parameters --- +! Presence of vcmx25 activates PFT-specific photosynthesis in CLM. +! Without these, CLM uses hardcoded defaults from clm_varcon.F90. +! +vcmx25 Maximum rate of carboxylation at 25C (umol CO2/m2/s) +51. 62. 39. 57. 54. 17. 17. 40. 40. 24. 52. 50. 50. 24. 0. 17. 0. 0. +! C3/C4 fix: IGBP 10 (GRA) 52→24 (C4), IGBP 14 50→24 (C4) +! +c3psn Photosynthetic pathway (1=C3, 0=C4) +1. 1. 1. 1. 1. 1. 1. 1. 1. 0. 1. 0. 1. 0. 0. 0. 0. 0. +! IGBP 10 (GRA)=C4, IGBP 12 (CRO)=C4 (temperate default) +! +mp Ball-Berry slope parameter +9. 9. 9. 9. 9. 9. 9. 9. 9. 4. 9. 4. 9. 4. 0. 9. 0. 0. +! C4 types use mp=4 (Collatz et al. 1992) +! +bp Minimum leaf conductance (umol/m2/s) +2000. 2000. 2000. 2000. 2000. 2000. 2000. 2000. 2000. 40000. 2000. 40000. 2000. 40000. 2000. 2000. 2000. 2000. +! C4 types use bp=40000 (Collatz et al. 1992) +! +qe25 Quantum efficiency at 25C (umol CO2/umol photon) +0.06 0.06 0.06 0.06 0.06 0.06 0.06 0.06 0.06 0.05 0.06 0.04 0.06 0.05 0. 0.06 0. 0. +! C4 fix: IGBP 10 0.04→0.05, IGBP 14 0.04→0.05 +! +folnmx Foliage nitrogen concentration when f(N)=1 (%) +1.5 1.5 1.5 1.5 1.5 1.5 1.5 1.5 1.5 1.5 1.5 1.5 1.5 1.5 1.5 1.5 1.5 1.5 +! +g1_medlyn Medlyn stomatal slope parameter (kPa^0.5) +2.35 4.12 2.35 4.45 3.50 4.70 4.70 4.45 4.45 1.62 5.25 1.79 4.00 1.62 0. 4.70 0. 0. +! Medlyn et al. (2011); values from Lin et al. (2015) +! C3/C4 fix: IGBP 9 1.62→4.45 (C3 savanna), IGBP 10 5.25→1.62 (C4), +! IGBP 12 5.79→1.79 (C4), IGBP 14 5.25→1.62 (C4) +! +clump Vegetation clumping index (0-1, 1=no clumping) +0.62 0.72 0.58 0.76 0.69 0.84 0.84 0.78 0.78 0.95 0.80 0.95 0.95 0.95 1.0 0.95 1.0 1.0 +! He et al. (2012); values by IGBP class from global dataset +! diff --git a/Single_Column_PFCLM_netcdf.ipynb b/et/Single_Column_PFCLM_netcdf.ipynb similarity index 77% rename from Single_Column_PFCLM_netcdf.ipynb rename to et/Single_Column_PFCLM_netcdf.ipynb index 4953a43..02c346b 100644 --- a/Single_Column_PFCLM_netcdf.ipynb +++ b/et/Single_Column_PFCLM_netcdf.ipynb @@ -4,19 +4,7 @@ "cell_type": "markdown", "id": "ban9v429qkf", "metadata": {}, - "source": [ - "# Single-Column ParFlow-CLM Simulation\n", - "\n", - "Run a single-column PF-CLM simulation with improved CLM physics:\n", - "- **Medlyn stomatal conductance** \u2014 VPD-based stomata (replaces Ball-Berry)\n", - "- **CLM5 tanh canopy interception** \u2014 physically-based interception scheme\n", - "- **Canopy clumping** \u2014 He et al. 2012 PFT-dependent clumping indices\n", - "\n", - "This notebook uses the forcing and CLM inputs prepared in Notebook 1.\n", - "Output is written as NetCDF for analysis in Notebook 3.\n", - "\n", - "**Requires:** ParFlow v3.14.1+ with CLM ET improvements (`feature/clm_et` branch or later)." - ], + "source": "# Single-Column ParFlow-CLM Simulation (ET example)\n\nRun a single-column PF-CLM simulation with the `clm_et_improved` physics defaults:\n- **Medlyn stomatal conductance** — VPD-based stomata (replaces Ball-Berry)\n- **CLM5 tanh canopy interception** — physically-based interception scheme\n- **Canopy clumping** — He et al. 2012 PFT-dependent clumping indices\n\nThis notebook uses the forcing and CLM inputs prepared in Notebook 1.\nOutput is written as NetCDF for analysis in Notebook 3.\n\n**Requires:** ParFlow master (commit containing PR #712 \"Add CLM ET formulation improvements: PFT photosynthesis, canopy interception, Medlyn stomata\", merged 2026-03-19) or later.\n", "outputs": [] }, { @@ -25,25 +13,7 @@ "id": "0ize3fa2df9", "metadata": {}, "outputs": [], - "source": [ - "import os\n", - "import numpy as np\n", - "import xarray as xr\n", - "import pandas as pd\n", - "import matplotlib.pyplot as plt\n", - "import matplotlib.dates as mdates\n", - "from pathlib import Path\n", - "from datetime import datetime\n", - "\n", - "from helpers import read_site_config\n", - "\n", - "# Set PARFLOW_DIR before importing parflow\n", - "os.environ[\"PARFLOW_DIR\"] = os.environ.get(\n", - " \"PARFLOW_DIR\", str(Path.home() / \"parflow\" / \"parflow_21_Jan_26\"))\n", - "print(f\"PARFLOW_DIR: {os.environ['PARFLOW_DIR']}\")\n", - "\n", - "from parflow import Run" - ] + "source": "import os\nimport sys\nimport numpy as np\nimport xarray as xr\nimport pandas as pd\nimport matplotlib.pyplot as plt\nimport matplotlib.dates as mdates\nfrom pathlib import Path\nfrom datetime import datetime\n\n# Make helpers.py (at repo root) importable from this subdir\nsys.path.insert(0, os.path.abspath(\"..\"))\nfrom helpers import read_site_config\n\n# Set PARFLOW_DIR before importing parflow\nos.environ[\"PARFLOW_DIR\"] = os.environ.get(\n \"PARFLOW_DIR\", str(Path.home() / \"parflow\" / \"parflow_21_Jan_26\"))\nprint(f\"PARFLOW_DIR: {os.environ['PARFLOW_DIR']}\")\n\nfrom parflow import Run\n" }, { "cell_type": "code", @@ -51,7 +21,7 @@ "id": "8qvxm6qie9x", "metadata": {}, "outputs": [], - "source": "# ============================================================\n# Load Site Configuration from Notebook 1\n# ============================================================\nstatic_write_dir = \"runs/US_Slt_WY2012\"\n\nrun_dir = Path(static_write_dir)\ncfg = read_site_config(run_dir / \"site_config.txt\")\n\n# Parse config values\nsite_id = cfg[\"site_id\"]\nwater_year = int(cfg[\"water_year\"])\nstart = cfg[\"start\"]\nend = cfg[\"end\"]\nlat = float(cfg[\"lat\"])\nlon = float(cfg[\"lon\"])\n\n# Soil (top 2m)\nsoil_ksat = float(cfg[\"soil_ksat\"])\nsoil_porosity = float(cfg[\"soil_porosity\"])\nsoil_vg_alpha = float(cfg[\"soil_vg_alpha\"])\nsoil_vg_n = float(cfg[\"soil_vg_n\"])\nsoil_sres = float(cfg[\"soil_sres\"])\n\n# Geology (bottom 6m)\ngeol_ksat = float(cfg[\"geol_ksat\"])\ngeol_porosity = float(cfg[\"geol_porosity\"])\ngeol_vg_alpha = float(cfg[\"geol_vg_alpha\"])\ngeol_vg_n = float(cfg[\"geol_vg_n\"])\ngeol_sres = float(cfg[\"geol_sres\"])\n\n# WTD and forcing\nwtd_m = float(cfg[\"wtd_m\"])\nforcing_filename = cfg[\"forcing_file\"]\n\n# ---- Override any config values here if needed ----\n# wtd_m = -0.2\n\n# Hours\ndt_start = datetime.strptime(start, \"%Y-%m-%d\")\ndt_end = datetime.strptime(end, \"%Y-%m-%d\")\nn_hours = int((dt_end - dt_start).total_seconds() / 3600)\n\nprint(f\"Site: {site_id} WY{water_year}, {n_hours} hours\")\nprint(f\"\\n Soil: Ksat={soil_ksat}, porosity={soil_porosity}, \"\n f\"alpha={soil_vg_alpha}, n={soil_vg_n}\")\nprint(f\" Geology: Ksat={geol_ksat}, porosity={geol_porosity}, \"\n f\"alpha={geol_vg_alpha}, n={geol_vg_n}\")\nprint(f\" WTD BC: {wtd_m} m\")\nprint(f\"\\nEdit values above or uncomment overrides, then run the next cell.\")" + "source": "# ============================================================\n# Load Site Configuration from Notebook 1\n# ============================================================\nstatic_write_dir = \"../runs/et/US_Slt_WY2012\"\n\nrun_dir = Path(static_write_dir)\ncfg = read_site_config(run_dir / \"site_config.txt\")\n\n# Parse config values\nsite_id = cfg[\"site_id\"]\nwater_year = int(cfg[\"water_year\"])\nstart = cfg[\"start\"]\nend = cfg[\"end\"]\nlat = float(cfg[\"lat\"])\nlon = float(cfg[\"lon\"])\n\n# Soil (top 2m)\nsoil_ksat = float(cfg[\"soil_ksat\"])\nsoil_porosity = float(cfg[\"soil_porosity\"])\nsoil_vg_alpha = float(cfg[\"soil_vg_alpha\"])\nsoil_vg_n = float(cfg[\"soil_vg_n\"])\nsoil_sres = float(cfg[\"soil_sres\"])\n\n# Geology (bottom 6m)\ngeol_ksat = float(cfg[\"geol_ksat\"])\ngeol_porosity = float(cfg[\"geol_porosity\"])\ngeol_vg_alpha = float(cfg[\"geol_vg_alpha\"])\ngeol_vg_n = float(cfg[\"geol_vg_n\"])\ngeol_sres = float(cfg[\"geol_sres\"])\n\n# WTD and forcing\nwtd_m = float(cfg[\"wtd_m\"])\nforcing_filename = cfg[\"forcing_file\"]\n\n# ---- Override any config values here if needed ----\n# wtd_m = -0.2\n\n# Hours\ndt_start = datetime.strptime(start, \"%Y-%m-%d\")\ndt_end = datetime.strptime(end, \"%Y-%m-%d\")\nn_hours = int((dt_end - dt_start).total_seconds() / 3600)\n\nprint(f\"Site: {site_id} WY{water_year}, {n_hours} hours\")\nprint(f\"\\n Soil: Ksat={soil_ksat}, porosity={soil_porosity}, \"\n f\"alpha={soil_vg_alpha}, n={soil_vg_n}\")\nprint(f\" Geology: Ksat={geol_ksat}, porosity={geol_porosity}, \"\n f\"alpha={geol_vg_alpha}, n={geol_vg_n}\")\nprint(f\" WTD BC: {wtd_m} m\")\nprint(f\"\\nEdit values above or uncomment overrides, then run the next cell.\")\n" }, { "cell_type": "code", @@ -67,7 +37,7 @@ "# Total physical depth: 8m (10 computational layers with variable dz)\n", "#\n", "# Computational grid (Z=0 bottom, Z=10 top):\n", - "# Layers 0-5 (Z=0-6): geology, 6 \u00d7 1.0m = 6.0m physical\n", + "# Layers 0-5 (Z=0-6): geology, 6 × 1.0m = 6.0m physical\n", "# Layer 6 (Z=6-7): soil, 1.0m physical\n", "# Layer 7 (Z=7-8): soil, 0.6m physical\n", "# Layer 8 (Z=8-9): soil, 0.3m physical\n", @@ -148,10 +118,10 @@ "run.Cell._3.dzScale.Value = 1.0 # geology\n", "run.Cell._4.dzScale.Value = 1.0 # geology\n", "run.Cell._5.dzScale.Value = 1.0 # geology\n", - "run.Cell._6.dzScale.Value = 1.0 # soil \u2014 1.0m\n", - "run.Cell._7.dzScale.Value = 0.6 # soil \u2014 0.6m\n", - "run.Cell._8.dzScale.Value = 0.3 # soil \u2014 0.3m\n", - "run.Cell._9.dzScale.Value = 0.1 # soil \u2014 0.1m (surface)\n", + "run.Cell._6.dzScale.Value = 1.0 # soil — 1.0m\n", + "run.Cell._7.dzScale.Value = 0.6 # soil — 0.6m\n", + "run.Cell._8.dzScale.Value = 0.3 # soil — 0.3m\n", + "run.Cell._9.dzScale.Value = 0.1 # soil — 0.1m (surface)\n", "\n", "# ----- HYDRAULIC PROPERTIES (two-layer) -----\n", "run.Geom.Perm.Names = \"soil geology\"\n", @@ -391,7 +361,7 @@ "id": "stmsqs0crob", "metadata": {}, "outputs": [], - "source": "# Plot ET components and soil moisture time series\ndates = pd.date_range(start, periods=n_steps, freq=\"h\")\n\nfig, axes = plt.subplots(3, 1, figsize=(14, 10), sharex=True)\n\n# Panel 1: Latent heat flux\naxes[0].plot(dates, clm_ds[\"eflx_lh_tot\"].values.flatten(), lw=0.3, color=\"green\")\naxes[0].set_ylabel(\"LH Flux [W/m2]\")\naxes[0].set_title(f\"{site_id} WY{water_year} \u2014 ParFlow-CLM Output ({n_steps}/{n_hours} hrs)\")\n\n# Panel 2: ET components (mm/day)\nfor var, label, color in [\n (\"qflx_tran_veg\", \"Transpiration\", \"green\"),\n (\"qflx_evap_soi\", \"Soil Evaporation\", \"brown\"),\n (\"qflx_evap_veg\", \"Canopy Evaporation\", \"blue\"),\n]:\n if var in clm_ds:\n vals = clm_ds[var].values.flatten() * 3600 * 24 # mm/s -> mm/day\n axes[1].plot(dates, vals, lw=0.3, color=color, label=label)\naxes[1].set_ylabel(\"Flux [mm/day]\")\naxes[1].legend(loc=\"upper right\", fontsize=8)\n\n# Panel 3: Soil moisture (top 2 layers, cell-centered depths)\nif pf_ds is not None and \"saturation\" in pf_ds:\n sat = pf_ds[\"saturation\"]\n pf_dates = pd.date_range(start, periods=sat.shape[0], freq=\"h\")\n axes[2].plot(pf_dates, sat[:, 9, 0, 0].values * soil_porosity * 100,\n lw=0.5, color=\"steelblue\", label=\"5 cm (cell center)\")\n axes[2].plot(pf_dates, sat[:, 8, 0, 0].values * soil_porosity * 100,\n lw=0.5, color=\"navy\", label=\"25 cm (cell center)\")\n axes[2].set_ylabel(\"Vol. Water Content [%]\")\n axes[2].legend(loc=\"upper right\", fontsize=8)\nelse:\n axes[2].text(0.5, 0.5, \"No saturation output\", transform=axes[2].transAxes,\n ha=\"center\", va=\"center\")\n\naxes[2].set_xlabel(\"Date\")\naxes[2].xaxis.set_major_formatter(mdates.DateFormatter(\"%b %Y\"))\nplt.tight_layout()\nplt.show()\n\nclm_ds.close()\nif pf_ds is not None:\n pf_ds.close()" + "source": "# Plot ET components and soil moisture time series\ndates = pd.date_range(start, periods=n_steps, freq=\"h\")\n\nfig, axes = plt.subplots(3, 1, figsize=(14, 10), sharex=True)\n\n# Panel 1: Latent heat flux\naxes[0].plot(dates, clm_ds[\"eflx_lh_tot\"].values.flatten(), lw=0.3, color=\"green\")\naxes[0].set_ylabel(\"LH Flux [W/m2]\")\naxes[0].set_title(f\"{site_id} WY{water_year} — ParFlow-CLM Output ({n_steps}/{n_hours} hrs)\")\n\n# Panel 2: ET components (mm/day)\nfor var, label, color in [\n (\"qflx_tran_veg\", \"Transpiration\", \"green\"),\n (\"qflx_evap_soi\", \"Soil Evaporation\", \"brown\"),\n (\"qflx_evap_veg\", \"Canopy Evaporation\", \"blue\"),\n]:\n if var in clm_ds:\n vals = clm_ds[var].values.flatten() * 3600 * 24 # mm/s -> mm/day\n axes[1].plot(dates, vals, lw=0.3, color=color, label=label)\naxes[1].set_ylabel(\"Flux [mm/day]\")\naxes[1].legend(loc=\"upper right\", fontsize=8)\n\n# Panel 3: Soil moisture (top 2 layers, cell-centered depths)\nif pf_ds is not None and \"saturation\" in pf_ds:\n sat = pf_ds[\"saturation\"]\n pf_dates = pd.date_range(start, periods=sat.shape[0], freq=\"h\")\n axes[2].plot(pf_dates, sat[:, 9, 0, 0].values * soil_porosity * 100,\n lw=0.5, color=\"steelblue\", label=\"5 cm (cell center)\")\n axes[2].plot(pf_dates, sat[:, 8, 0, 0].values * soil_porosity * 100,\n lw=0.5, color=\"navy\", label=\"25 cm (cell center)\")\n axes[2].set_ylabel(\"Vol. Water Content [%]\")\n axes[2].legend(loc=\"upper right\", fontsize=8)\nelse:\n axes[2].text(0.5, 0.5, \"No saturation output\", transform=axes[2].transAxes,\n ha=\"center\", va=\"center\")\n\naxes[2].set_xlabel(\"Date\")\naxes[2].xaxis.set_major_formatter(mdates.DateFormatter(\"%b %Y\"))\nplt.tight_layout()\nplt.show()\n\nclm_ds.close()\nif pf_ds is not None:\n pf_ds.close()" } ], "metadata": { diff --git a/locate_plot_station_get_forcing.ipynb b/et/locate_plot_station_get_forcing.ipynb similarity index 84% rename from locate_plot_station_get_forcing.ipynb rename to et/locate_plot_station_get_forcing.ipynb index dd7010d..ae141e6 100644 --- a/locate_plot_station_get_forcing.ipynb +++ b/et/locate_plot_station_get_forcing.ipynb @@ -13,7 +13,7 @@ "3. Prepare CLM input files (vegm, clmin)\n", "4. Download CW3E meteorological forcing\n", "\n", - "**Prerequisites:** HydroData account (free) \u2014 see README.md for signup instructions." + "**Prerequisites:** HydroData account (free) — see README.md for signup instructions." ], "outputs": [] }, @@ -23,22 +23,7 @@ "id": "0ckol3tfby8", "metadata": {}, "outputs": [], - "source": [ - "# Imports\n", - "import numpy as np\n", - "import pandas as pd\n", - "import matplotlib.pyplot as plt\n", - "import matplotlib.dates as mdates\n", - "import subsettools as st\n", - "import hf_hydrodata as hf\n", - "import shutil\n", - "from pathlib import Path\n", - "from datetime import datetime\n", - "\n", - "from helpers import (IGBP_NAMES, IGBP_FULL_NAMES, IGBP_STR_TO_INT,\n", - " CONUS2_SOILS, CONUS2_GEOLOGY, CONUS2_SUBSURFACE,\n", - " build_vegm, patch_clmin_dates, write_site_config)" - ] + "source": "# Imports\nimport os\nimport sys\nimport numpy as np\nimport pandas as pd\nimport matplotlib.pyplot as plt\nimport matplotlib.dates as mdates\nimport subsettools as st\nimport hf_hydrodata as hf\nimport shutil\nfrom pathlib import Path\nfrom datetime import datetime\n\n# Make helpers.py (at repo root) importable from this subdir\nsys.path.insert(0, os.path.abspath(\"..\"))\nfrom helpers import (IGBP_NAMES, IGBP_FULL_NAMES, IGBP_STR_TO_INT,\n CONUS2_SOILS, CONUS2_GEOLOGY, CONUS2_SUBSURFACE,\n build_vegm, patch_clmin_dates, write_site_config)\n" }, { "cell_type": "code", @@ -46,7 +31,7 @@ "id": "xjgitxx73di", "metadata": {}, "outputs": [], - "source": "# ============================================================\n# User Configuration\n# ============================================================\n# Choose a site and water year. These are used throughout all 3 notebooks.\n\nsite_id = \"US-Slt\" # Ameriflux site ID\nstart = \"2011-10-01\" # Water year start (Oct 1)\nend = \"2012-10-01\" # Water year end (Oct 1)\nwater_year = 2012\n\n# Directory where forcing and CLM input files will be written\n# This is also where you will run ParFlow in Notebook 2\nstatic_write_dir = f\"runs/{site_id.replace('-','_')}_WY{water_year}\"\n\n# Forcing filename (constructed from site_id and dates)\nforcing_filename = f\"forcing1D.{site_id}.{start}-{end}.txt\"\n\n# HydroData PIN (only needed once per machine)\n# Sign up: https://hydrogen.princeton.edu/signup\n# Get PIN: https://hydrogen.princeton.edu/pin\nhf.register_api_pin(\"your_email@example.com\", \"your_pin\")\n\n# -----------------------------------------------------------\n# Tested sites from the ET sensitivity study for reference:\n# US-Slt (DBF, NJ, WY2012) US-xUK (DBF, KS, WY2024)\n# US-xBL (DBF, VA, WY2024) US-GLE (ENF, WY, WY2020)\n# US-Ro5 (CRO, MN, WY2024) US-Ro1 (CRO, MN, WY2016)\n# US-Ro2 (CRO, MN, WY2016) US-Fwf (GRA, AZ, WY2008)\n# US-Mpj (WSA, NM, WY2024) US-Ton (WSA, CA, WY2024)\n# US-xSR (OSH, AZ, WY2024) US-xJR (OSH, NM, WY2024)\n# US-xDL (MF, AL, WY2024)\n# -----------------------------------------------------------" + "source": "# ============================================================\n# User Configuration\n# ============================================================\n# Choose a site and water year. These are used throughout all 3 notebooks.\n\nsite_id = \"US-Slt\" # Ameriflux site ID\nstart = \"2011-10-01\" # Water year start (Oct 1)\nend = \"2012-10-01\" # Water year end (Oct 1)\nwater_year = 2012\n\n# Directory where forcing and CLM input files will be written\n# (path is relative to this notebook's directory, i.e. et/)\n# This is also where you will run ParFlow in Notebook 2\nstatic_write_dir = f\"../runs/et/{site_id.replace('-','_')}_WY{water_year}\"\n\n# Forcing filename (constructed from site_id and dates)\nforcing_filename = f\"forcing1D.{site_id}.{start}-{end}.txt\"\n\n# HydroData PIN (only needed once per machine)\n# Sign up: https://hydrogen.princeton.edu/signup\n# Get PIN: https://hydrogen.princeton.edu/pin\nhf.register_api_pin(\"your_email@example.com\", \"your_pin\")\n\n# -----------------------------------------------------------\n# Tested sites from the ET sensitivity study for reference:\n# US-Slt (DBF, NJ, WY2012) US-xUK (DBF, KS, WY2024)\n# US-xBL (DBF, VA, WY2024) US-GLE (ENF, WY, WY2020)\n# US-Ro5 (CRO, MN, WY2024) US-Ro1 (CRO, MN, WY2016)\n# US-Ro2 (CRO, MN, WY2016) US-Fwf (GRA, AZ, WY2008)\n# US-Mpj (WSA, NM, WY2024) US-Ton (WSA, CA, WY2024)\n# US-xSR (OSH, AZ, WY2024) US-xJR (OSH, NM, WY2024)\n# US-xDL (MF, AL, WY2024)\n# -----------------------------------------------------------\n" }, { "cell_type": "code", @@ -206,7 +191,7 @@ " \"dataset\": \"conus2_domain\", \"variable\": \"pf_indicator\",\n", " \"grid_bounds\": ij_bounds,\n", "})\n", - "print(f\"\\nCONUS2 subsurface profile (bottom \u2192 top):\")\n", + "print(f\"\\nCONUS2 subsurface profile (bottom → top):\")\n", "for i in range(ind.shape[0]):\n", " t = int(ind[i, 0, 0])\n", " info = CONUS2_SUBSURFACE.get(t, (\"Unknown\",))\n", @@ -275,7 +260,7 @@ "outputs": [], "source": [ "# ============================================================\n", - "# Vegm parameters \u2014 edit these if needed\n", + "# Vegm parameters — edit these if needed\n", "# ============================================================\n", "# Sand/clay: try to get from subsettools config_clm, fall back to defaults\n", "import tempfile, os\n", @@ -304,7 +289,7 @@ "\n", "igbp_full = IGBP_FULL_NAMES.get(igbp_type, igbp_str)\n", "print(f\"\\n--- Vegetation / Soil Texture ---\")\n", - "print(f\" IGBP type: {igbp_type} \u2014 {igbp_full} ({IGBP_NAMES.get(igbp_type, '?')})\")\n", + "print(f\" IGBP type: {igbp_type} — {igbp_full} ({IGBP_NAMES.get(igbp_type, '?')})\")\n", "print(f\" Sand: {sand_frac:.0%}\")\n", "print(f\" Clay: {clay_frac:.0%}\")\n", "print(f\" Soil color: {color_index}\")" @@ -327,28 +312,7 @@ "id": "7qmlife8dzd", "metadata": {}, "outputs": [], - "source": [ - "# Create run directory and write CLM input files\n", - "run_dir = Path(static_write_dir)\n", - "run_dir.mkdir(parents=True, exist_ok=True)\n", - "\n", - "# 1. Write drv_vegm.dat\n", - "build_vegm(run_dir / \"drv_vegm.dat\", lat, lon, sand_frac, clay_frac, color_index, igbp_type)\n", - "print(f\"Wrote drv_vegm.dat (IGBP={igbp_type})\")\n", - "\n", - "# 2. Patch drv_clmin.dat dates from template\n", - "clm_inputs = Path(\"clm_inputs\")\n", - "patch_clmin_dates(clm_inputs / \"drv_clmin_template.dat\",\n", - " run_dir / \"drv_clmin.dat\", water_year)\n", - "print(f\"Wrote drv_clmin.dat (WY{water_year})\")\n", - "\n", - "# 3. Copy drv_vegp.dat (shipped with improved PFT parameters + clumping)\n", - "shutil.copy2(clm_inputs / \"drv_vegp.dat\", run_dir / \"drv_vegp.dat\")\n", - "print(f\"Copied drv_vegp.dat (CLM4.5 + canopy clumping)\")\n", - "\n", - "print(f\"\\nRun directory: {run_dir}\")\n", - "print(\"Contents:\", [f.name for f in sorted(run_dir.iterdir())])" - ] + "source": "# Create run directory and write CLM input files\nrun_dir = Path(static_write_dir)\nrun_dir.mkdir(parents=True, exist_ok=True)\n\n# 1. Write drv_vegm.dat\nbuild_vegm(run_dir / \"drv_vegm.dat\", lat, lon, sand_frac, clay_frac, color_index, igbp_type)\nprint(f\"Wrote drv_vegm.dat (IGBP={igbp_type})\")\n\n# 2. Patch drv_clmin.dat dates from shared template\nclm_inputs_shared = Path(\"../clm_inputs\")\npatch_clmin_dates(clm_inputs_shared / \"drv_clmin_template.dat\",\n run_dir / \"drv_clmin.dat\", water_year)\nprint(f\"Wrote drv_clmin.dat (WY{water_year})\")\n\n# 3. Copy drv_vegp.dat (ET-specific: Medlyn + clumping + CLM4.5 optics)\nclm_inputs_et = Path(\"../clm_inputs/et\")\nshutil.copy2(clm_inputs_et / \"drv_vegp.dat\", run_dir / \"drv_vegp.dat\")\nprint(f\"Copied drv_vegp.dat (clm_et_improved config)\")\n\nprint(f\"\\nRun directory: {run_dir}\")\nprint(\"Contents:\", [f.name for f in sorted(run_dir.iterdir())])\n" }, { "cell_type": "markdown", @@ -463,7 +427,7 @@ "metadata": {}, "outputs": [], "source": [ - "# Write site_config.txt \u2014 Notebook 2 reads this\n", + "# Write site_config.txt — Notebook 2 reads this\n", "config_path = run_dir / \"site_config.txt\"\n", "\n", "# Build params dict\n", @@ -514,7 +478,7 @@ "id": "zzdlaftndxe", "metadata": {}, "source": [ - "Proceed to **Single_Column_PFCLM_netcdf.ipynb** \u2014 it will load `site_config.txt` and display the parameters for you to review before running." + "Proceed to **Single_Column_PFCLM_netcdf.ipynb** — it will load `site_config.txt` and display the parameters for you to review before running." ], "outputs": [] } diff --git a/pfclm_ameriflux_compare.ipynb b/et/pfclm_ameriflux_compare.ipynb similarity index 54% rename from pfclm_ameriflux_compare.ipynb rename to et/pfclm_ameriflux_compare.ipynb index c3ed325..32832bb 100644 --- a/pfclm_ameriflux_compare.ipynb +++ b/et/pfclm_ameriflux_compare.ipynb @@ -23,7 +23,7 @@ "id": "z5wcywg9pu", "metadata": {}, "outputs": [], - "source": "import numpy as np\nimport pandas as pd\nimport xarray as xr\nimport matplotlib.pyplot as plt\nimport matplotlib.dates as mdates\nimport hf_hydrodata as hf\nfrom pathlib import Path\nfrom datetime import datetime\n\nfrom helpers import read_site_config\n\n# HydroData PIN (persists from Notebook 1, but register again if needed)\nhf.register_api_pin(\"your_email@example.com\", \"your_pin\")" + "source": "import os\nimport sys\nimport numpy as np\nimport pandas as pd\nimport xarray as xr\nimport matplotlib.pyplot as plt\nimport matplotlib.dates as mdates\nimport hf_hydrodata as hf\nfrom pathlib import Path\nfrom datetime import datetime\n\n# Make helpers.py (at repo root) importable from this subdir\nsys.path.insert(0, os.path.abspath(\"..\"))\nfrom helpers import read_site_config\n\n# HydroData PIN (persists from Notebook 1, but register again if needed)\nhf.register_api_pin(\"your_email@example.com\", \"your_pin\")\n" }, { "cell_type": "code", @@ -31,7 +31,7 @@ "id": "n5foxl9nkrs", "metadata": {}, "outputs": [], - "source": "# ============================================================\n# Load Site Configuration\n# ============================================================\nstatic_write_dir = \"runs/US_Slt_WY2012\"\n\nrun_dir = Path(static_write_dir)\ncfg = read_site_config(run_dir / \"site_config.txt\")\n\nsite_id = cfg[\"site_id\"]\nwater_year = int(cfg[\"water_year\"])\nstart = cfg[\"start\"]\nend = cfg[\"end\"]\nlat = float(cfg[\"lat\"])\nlon = float(cfg[\"lon\"])\nutc_offset = int(cfg[\"utc_offset\"])\nporosity = float(cfg[\"soil_porosity\"])\n\nprint(f\"Site: {site_id}, WY{water_year}, UTC offset: {utc_offset}\")" + "source": "# ============================================================\n# Load Site Configuration\n# ============================================================\nstatic_write_dir = \"../runs/et/US_Slt_WY2012\"\n\nrun_dir = Path(static_write_dir)\ncfg = read_site_config(run_dir / \"site_config.txt\")\n\nsite_id = cfg[\"site_id\"]\nwater_year = int(cfg[\"water_year\"])\nstart = cfg[\"start\"]\nend = cfg[\"end\"]\nlat = float(cfg[\"lat\"])\nlon = float(cfg[\"lon\"])\nutc_offset = int(cfg[\"utc_offset\"])\nporosity = float(cfg[\"soil_porosity\"])\n\nprint(f\"Site: {site_id}, WY{water_year}, UTC offset: {utc_offset}\")\n" }, { "cell_type": "code", @@ -162,7 +162,7 @@ "mod_et = float((mod_c / lv * 86400).sum())\n", "\n", "print(f\"\\\\n{'='*50}\")\n", - "print(f\" {site_id} WY{water_year} \u2014 Daily LH Metrics\")\n", + "print(f\" {site_id} WY{water_year} — Daily LH Metrics\")\n", "print(f\"{'='*50}\")\n", "print(f\" N days: {len(common)}\")\n", "print(f\" Bias: {bias:+.1f} W/m2\")\n", @@ -183,7 +183,7 @@ "outputs": [], "source": [ "# ============================================================\n", - "# Time Series Plot \u2014 Daily Latent Heat\n", + "# Time Series Plot — Daily Latent Heat\n", "# ============================================================\n", "fig, ax = plt.subplots(figsize=(14, 5))\n", "\n", @@ -191,7 +191,7 @@ "ax.plot(mod_c.index, mod_c.values, color=\"steelblue\", lw=1.0, label=\"PF-CLM Model\", alpha=0.8)\n", "\n", "ax.set_ylabel(\"Latent Heat [W/m2]\")\n", - "ax.set_title(f\"{site_id} WY{water_year} \u2014 Daily Mean LH | \"\n", + "ax.set_title(f\"{site_id} WY{water_year} — Daily Mean LH | \"\n", " f\"KGE={kge:.2f}, Bias={bias:+.1f} W/m2, RMSE={rmse:.1f}\")\n", "ax.legend(loc=\"upper right\")\n", "ax.xaxis.set_major_formatter(mdates.DateFormatter(\"%b %Y\"))\n", @@ -208,7 +208,7 @@ "outputs": [], "source": [ "# ============================================================\n", - "# Scatter Plot \u2014 Observed vs Modeled\n", + "# Scatter Plot — Observed vs Modeled\n", "# ============================================================\n", "fig, ax = plt.subplots(figsize=(7, 7))\n", "\n", @@ -237,7 +237,7 @@ "outputs": [], "source": [ "# ============================================================\n", - "# ET Components \u2014 Stacked Area + Obs Overlay\n", + "# ET Components — Stacked Area + Obs Overlay\n", "# ============================================================\n", "et_vars = {\n", " \"qflx_tran_veg\": (\"Transpiration\", \"green\"),\n", @@ -264,7 +264,7 @@ "\n", "ax.plot(obs_c.index, obs_c.values, \"k-\", lw=1.0, label=\"Ameriflux Obs\")\n", "ax.set_ylabel(\"Latent Heat [W/m2]\")\n", - "ax.set_title(f\"{site_id} WY{water_year} \u2014 ET Components\")\n", + "ax.set_title(f\"{site_id} WY{water_year} — ET Components\")\n", "ax.legend(loc=\"upper right\", fontsize=8)\n", "ax.xaxis.set_major_formatter(mdates.DateFormatter(\"%b %Y\"))\n", "plt.tight_layout()\n", @@ -277,7 +277,7 @@ "id": "d11n1t1hn4i", "metadata": {}, "outputs": [], - "source": "# ============================================================\n# Soil Moisture \u2014 Top 2 Layers from ParFlow\n# ============================================================\nif pf_files:\n pf_ds = xr.open_dataset(pf_files[0])\n if \"saturation\" in pf_ds:\n sat = pf_ds[\"saturation\"]\n pf_dates = pd.date_range(start, periods=sat.shape[0], freq=\"h\")\n\n fig, ax = plt.subplots(figsize=(14, 4))\n ax.plot(pf_dates, sat[:, 9, 0, 0].values * porosity * 100,\n lw=0.5, color=\"steelblue\", label=\"5 cm (cell center)\")\n ax.plot(pf_dates, sat[:, 8, 0, 0].values * porosity * 100,\n lw=0.5, color=\"navy\", label=\"25 cm (cell center)\")\n ax.set_ylabel(\"Vol. Water Content [%]\")\n ax.set_title(f\"{site_id} WY{water_year} \u2014 Soil Moisture (porosity={porosity})\")\n ax.legend()\n ax.xaxis.set_major_formatter(mdates.DateFormatter(\"%b %Y\"))\n plt.tight_layout()\n plt.show()\n pf_ds.close()\nelse:\n print(\"No ParFlow pressure/saturation NetCDF output found.\")" + "source": "# ============================================================\n# Soil Moisture — Top 2 Layers from ParFlow\n# ============================================================\nif pf_files:\n pf_ds = xr.open_dataset(pf_files[0])\n if \"saturation\" in pf_ds:\n sat = pf_ds[\"saturation\"]\n pf_dates = pd.date_range(start, periods=sat.shape[0], freq=\"h\")\n\n fig, ax = plt.subplots(figsize=(14, 4))\n ax.plot(pf_dates, sat[:, 9, 0, 0].values * porosity * 100,\n lw=0.5, color=\"steelblue\", label=\"5 cm (cell center)\")\n ax.plot(pf_dates, sat[:, 8, 0, 0].values * porosity * 100,\n lw=0.5, color=\"navy\", label=\"25 cm (cell center)\")\n ax.set_ylabel(\"Vol. Water Content [%]\")\n ax.set_title(f\"{site_id} WY{water_year} — Soil Moisture (porosity={porosity})\")\n ax.legend()\n ax.xaxis.set_major_formatter(mdates.DateFormatter(\"%b %Y\"))\n plt.tight_layout()\n plt.show()\n pf_ds.close()\nelse:\n print(\"No ParFlow pressure/saturation NetCDF output found.\")" }, { "cell_type": "code", @@ -285,7 +285,7 @@ "id": "aguq5vz6n2b", "metadata": {}, "outputs": [], - "source": "# ============================================================\n# Forcing Quality \u2014 CW3E vs Ameriflux (5 variables)\n# ============================================================\n# Compare: SWdown, LWdown, Tair, Wind speed, VPD\n# Not all variables are available at all sites \u2014 skips missing ones.\n\n# Load CW3E forcing\nforcing_path = run_dir / cfg[\"forcing_file\"]\nforcing_raw = np.loadtxt(forcing_path)\nforc_dates = pd.date_range(start, periods=len(forcing_raw), freq=\"h\")\n# Shift to local time for diurnal/daily comparisons\nforc_dates_local = forc_dates + pd.Timedelta(hours=utc_offset)\n\nforcing_df = pd.DataFrame(forcing_raw, index=forc_dates_local,\n columns=[\"dswr\", \"dlwr\", \"precip\", \"temp_k\", \"ugrd\", \"vgrd\", \"press\", \"spfh\"])\nforcing_df[\"temp_c\"] = forcing_df[\"temp_k\"] - 273.15\nforcing_df[\"wind_speed\"] = np.sqrt(forcing_df[\"ugrd\"]**2 + forcing_df[\"vgrd\"]**2)\n\n# VPD: e = q*P/(0.622+0.378*q), esat = 611.2*exp(17.67*Tc/(Tc+243.5)), VPD = (esat-e)/100\nq, P, Tc = forcing_df[\"spfh\"], forcing_df[\"press\"], forcing_df[\"temp_c\"]\ne = q * P / (0.622 + 0.378 * q)\nesat = 611.2 * np.exp(17.67 * Tc / (Tc + 243.5))\nforcing_df[\"vpd_hpa\"] = ((esat - e) / 100.0).clip(lower=0)\n\n# Variables to compare: (display_name, ameriflux_var, aggregation, forcing_col, unit)\nvar_map = [\n (\"SWdown\", \"downward_shortwave\", \"mean\", \"dswr\", \"W/m2\"),\n (\"LWdown\", \"downward_longwave\", \"mean\", \"dlwr\", \"W/m2\"),\n (\"Tair\", \"air_temp\", \"mean\", \"temp_c\", \"C\"),\n (\"Wind\", \"wind_speed\", \"mean\", \"wind_speed\", \"m/s\"),\n (\"VPD\", \"vapor_pressure_deficit\", \"mean\", \"vpd_hpa\", \"hPa\"),\n]\n\n# Fetch obs and compute metrics for each variable\nresults = {}\nobs_data = {}\nprint(f\"Forcing quality check: {site_id} WY{water_year}\\n\")\n\nfor disp, obs_var, agg, forc_col, unit in var_map:\n try:\n df = hf.get_point_data(\n dataset=\"ameriflux\", variable=obs_var,\n temporal_resolution=\"hourly\", aggregation=agg,\n date_start=start, date_end=end, site_ids=site_id,\n )\n if df is None or df.empty or site_id not in df.columns:\n print(f\" {disp:<8}: no data available\")\n continue\n obs = pd.Series(df[site_id].values.astype(float),\n index=pd.to_datetime(df[\"date\"])).replace(-9999, np.nan).dropna()\n # Shift to local time\n obs.index = obs.index + pd.Timedelta(hours=utc_offset)\n # Temperature: Ameriflux is in C, forcing temp_c is also C\n # But if obs_var is air_temp, Ameriflux may be in K \u2014 check range\n if obs_var == \"air_temp\" and obs.median() > 200:\n obs = obs - 273.15 # K -> C\n\n obs_data[disp] = obs\n forc = forcing_df[forc_col]\n\n common = forc.index.intersection(obs.index)\n f_c = forc.loc[common].dropna()\n o_c = obs.loc[common].dropna()\n common = f_c.index.intersection(o_c.index)\n if len(common) < 100:\n print(f\" {disp:<8}: insufficient overlap ({len(common)} hours)\")\n continue\n f_c, o_c = f_c.loc[common], o_c.loc[common]\n\n bias = float((f_c - o_c).mean())\n rmse = float(np.sqrt(((f_c - o_c)**2).mean()))\n r = float(np.corrcoef(f_c, o_c)[0, 1])\n results[disp] = {\"bias\": bias, \"rmse\": rmse, \"r\": r, \"unit\": unit, \"n\": len(common)}\n print(f\" {disp:<8}: bias={bias:+.2f} {unit}, RMSE={rmse:.2f}, r={r:.2f} ({len(common)} hrs)\")\n except Exception as e:\n print(f\" {disp:<8}: skipped ({e})\")\n\nif not results:\n print(\"\\nNo forcing comparison data available for this site.\")\nelse:\n # Multi-panel plot: daily time series, diurnal (May-Sep), scatter\n vars_ok = [v for v in [\"SWdown\", \"LWdown\", \"Tair\", \"Wind\", \"VPD\"] if v in results]\n n_vars = len(vars_ok)\n\n fig, axes = plt.subplots(n_vars, 3, figsize=(18, 4 * n_vars),\n gridspec_kw={\"width_ratios\": [3, 1.5, 1.5]})\n if n_vars == 1:\n axes = axes.reshape(1, -1)\n\n forc_col_map = {\"SWdown\": \"dswr\", \"LWdown\": \"dlwr\", \"Tair\": \"temp_c\",\n \"Wind\": \"wind_speed\", \"VPD\": \"vpd_hpa\"}\n\n for row, var in enumerate(vars_ok):\n m = results[var]\n forc = forcing_df[forc_col_map[var]]\n obs = obs_data[var]\n\n common = forc.index.intersection(obs.dropna().index)\n f_c = forc.loc[common]\n o_c = obs.loc[common].dropna()\n common = f_c.index.intersection(o_c.index)\n f_c, o_c = f_c.loc[common], o_c.loc[common]\n\n # Panel 1: Daily time series\n ax = axes[row, 0]\n ax.plot(o_c.resample(\"D\").mean(), \"k-\", lw=1.5, label=\"Ameriflux\")\n ax.plot(f_c.resample(\"D\").mean(), \"r-\", lw=1, alpha=0.8, label=\"CW3E\")\n ax.set_ylabel(f\"{var} ({m['unit']})\")\n ax.set_title(f\"{var}: bias={m['bias']:+.2f}, RMSE={m['rmse']:.2f}, r={m['r']:.2f}\")\n ax.legend(fontsize=8)\n ax.xaxis.set_major_formatter(mdates.DateFormatter(\"%b\"))\n ax.grid(True, alpha=0.3)\n\n # Panel 2: Diurnal cycle (May-Sep)\n ax = axes[row, 1]\n gs = (common.month >= 5) & (common.month <= 9)\n if gs.sum() > 100:\n ax.plot(o_c.loc[common[gs]].groupby(o_c.loc[common[gs]].index.hour).mean(),\n \"k-\", lw=2, label=\"Obs\")\n ax.plot(f_c.loc[common[gs]].groupby(f_c.loc[common[gs]].index.hour).mean(),\n \"r-\", lw=1.5, alpha=0.8, label=\"CW3E\")\n ax.set_xlabel(\"Local Hour\")\n ax.set_title(\"Diurnal (May-Sep)\")\n ax.set_xticks(range(0, 24, 6))\n ax.legend(fontsize=8)\n ax.grid(True, alpha=0.3)\n\n # Panel 3: Scatter\n ax = axes[row, 2]\n step = max(1, len(f_c) // 2000)\n ax.scatter(o_c.values[::step], f_c.values[::step], s=3, alpha=0.3, color=\"steelblue\")\n lims = [min(o_c.min(), f_c.min()), max(o_c.max(), f_c.max())]\n ax.plot(lims, lims, \"k--\", alpha=0.5)\n ax.set_xlabel(f\"Obs ({m['unit']})\")\n ax.set_ylabel(f\"CW3E ({m['unit']})\")\n ax.set_title(f\"r={m['r']:.2f}\")\n ax.set_aspect(\"equal\")\n ax.grid(True, alpha=0.3)\n\n fig.suptitle(f\"{site_id} WY{water_year} \u2014 CW3E Forcing vs Ameriflux Observations\",\n fontsize=14, fontweight=\"bold\")\n plt.tight_layout()\n plt.show()\n\nclm_ds.close()" + "source": "# ============================================================\n# Forcing Quality — CW3E vs Ameriflux (5 variables)\n# ============================================================\n# Compare: SWdown, LWdown, Tair, Wind speed, VPD\n# Not all variables are available at all sites — skips missing ones.\n\n# Load CW3E forcing\nforcing_path = run_dir / cfg[\"forcing_file\"]\nforcing_raw = np.loadtxt(forcing_path)\nforc_dates = pd.date_range(start, periods=len(forcing_raw), freq=\"h\")\n# Shift to local time for diurnal/daily comparisons\nforc_dates_local = forc_dates + pd.Timedelta(hours=utc_offset)\n\nforcing_df = pd.DataFrame(forcing_raw, index=forc_dates_local,\n columns=[\"dswr\", \"dlwr\", \"precip\", \"temp_k\", \"ugrd\", \"vgrd\", \"press\", \"spfh\"])\nforcing_df[\"temp_c\"] = forcing_df[\"temp_k\"] - 273.15\nforcing_df[\"wind_speed\"] = np.sqrt(forcing_df[\"ugrd\"]**2 + forcing_df[\"vgrd\"]**2)\n\n# VPD: e = q*P/(0.622+0.378*q), esat = 611.2*exp(17.67*Tc/(Tc+243.5)), VPD = (esat-e)/100\nq, P, Tc = forcing_df[\"spfh\"], forcing_df[\"press\"], forcing_df[\"temp_c\"]\ne = q * P / (0.622 + 0.378 * q)\nesat = 611.2 * np.exp(17.67 * Tc / (Tc + 243.5))\nforcing_df[\"vpd_hpa\"] = ((esat - e) / 100.0).clip(lower=0)\n\n# Variables to compare: (display_name, ameriflux_var, aggregation, forcing_col, unit)\nvar_map = [\n (\"SWdown\", \"downward_shortwave\", \"mean\", \"dswr\", \"W/m2\"),\n (\"LWdown\", \"downward_longwave\", \"mean\", \"dlwr\", \"W/m2\"),\n (\"Tair\", \"air_temp\", \"mean\", \"temp_c\", \"C\"),\n (\"Wind\", \"wind_speed\", \"mean\", \"wind_speed\", \"m/s\"),\n (\"VPD\", \"vapor_pressure_deficit\", \"mean\", \"vpd_hpa\", \"hPa\"),\n]\n\n# Fetch obs and compute metrics for each variable\nresults = {}\nobs_data = {}\nprint(f\"Forcing quality check: {site_id} WY{water_year}\\n\")\n\nfor disp, obs_var, agg, forc_col, unit in var_map:\n try:\n df = hf.get_point_data(\n dataset=\"ameriflux\", variable=obs_var,\n temporal_resolution=\"hourly\", aggregation=agg,\n date_start=start, date_end=end, site_ids=site_id,\n )\n if df is None or df.empty or site_id not in df.columns:\n print(f\" {disp:<8}: no data available\")\n continue\n obs = pd.Series(df[site_id].values.astype(float),\n index=pd.to_datetime(df[\"date\"])).replace(-9999, np.nan).dropna()\n # Shift to local time\n obs.index = obs.index + pd.Timedelta(hours=utc_offset)\n # Temperature: Ameriflux is in C, forcing temp_c is also C\n # But if obs_var is air_temp, Ameriflux may be in K — check range\n if obs_var == \"air_temp\" and obs.median() > 200:\n obs = obs - 273.15 # K -> C\n\n obs_data[disp] = obs\n forc = forcing_df[forc_col]\n\n common = forc.index.intersection(obs.index)\n f_c = forc.loc[common].dropna()\n o_c = obs.loc[common].dropna()\n common = f_c.index.intersection(o_c.index)\n if len(common) < 100:\n print(f\" {disp:<8}: insufficient overlap ({len(common)} hours)\")\n continue\n f_c, o_c = f_c.loc[common], o_c.loc[common]\n\n bias = float((f_c - o_c).mean())\n rmse = float(np.sqrt(((f_c - o_c)**2).mean()))\n r = float(np.corrcoef(f_c, o_c)[0, 1])\n results[disp] = {\"bias\": bias, \"rmse\": rmse, \"r\": r, \"unit\": unit, \"n\": len(common)}\n print(f\" {disp:<8}: bias={bias:+.2f} {unit}, RMSE={rmse:.2f}, r={r:.2f} ({len(common)} hrs)\")\n except Exception as e:\n print(f\" {disp:<8}: skipped ({e})\")\n\nif not results:\n print(\"\\nNo forcing comparison data available for this site.\")\nelse:\n # Multi-panel plot: daily time series, diurnal (May-Sep), scatter\n vars_ok = [v for v in [\"SWdown\", \"LWdown\", \"Tair\", \"Wind\", \"VPD\"] if v in results]\n n_vars = len(vars_ok)\n\n fig, axes = plt.subplots(n_vars, 3, figsize=(18, 4 * n_vars),\n gridspec_kw={\"width_ratios\": [3, 1.5, 1.5]})\n if n_vars == 1:\n axes = axes.reshape(1, -1)\n\n forc_col_map = {\"SWdown\": \"dswr\", \"LWdown\": \"dlwr\", \"Tair\": \"temp_c\",\n \"Wind\": \"wind_speed\", \"VPD\": \"vpd_hpa\"}\n\n for row, var in enumerate(vars_ok):\n m = results[var]\n forc = forcing_df[forc_col_map[var]]\n obs = obs_data[var]\n\n common = forc.index.intersection(obs.dropna().index)\n f_c = forc.loc[common]\n o_c = obs.loc[common].dropna()\n common = f_c.index.intersection(o_c.index)\n f_c, o_c = f_c.loc[common], o_c.loc[common]\n\n # Panel 1: Daily time series\n ax = axes[row, 0]\n ax.plot(o_c.resample(\"D\").mean(), \"k-\", lw=1.5, label=\"Ameriflux\")\n ax.plot(f_c.resample(\"D\").mean(), \"r-\", lw=1, alpha=0.8, label=\"CW3E\")\n ax.set_ylabel(f\"{var} ({m['unit']})\")\n ax.set_title(f\"{var}: bias={m['bias']:+.2f}, RMSE={m['rmse']:.2f}, r={m['r']:.2f}\")\n ax.legend(fontsize=8)\n ax.xaxis.set_major_formatter(mdates.DateFormatter(\"%b\"))\n ax.grid(True, alpha=0.3)\n\n # Panel 2: Diurnal cycle (May-Sep)\n ax = axes[row, 1]\n gs = (common.month >= 5) & (common.month <= 9)\n if gs.sum() > 100:\n ax.plot(o_c.loc[common[gs]].groupby(o_c.loc[common[gs]].index.hour).mean(),\n \"k-\", lw=2, label=\"Obs\")\n ax.plot(f_c.loc[common[gs]].groupby(f_c.loc[common[gs]].index.hour).mean(),\n \"r-\", lw=1.5, alpha=0.8, label=\"CW3E\")\n ax.set_xlabel(\"Local Hour\")\n ax.set_title(\"Diurnal (May-Sep)\")\n ax.set_xticks(range(0, 24, 6))\n ax.legend(fontsize=8)\n ax.grid(True, alpha=0.3)\n\n # Panel 3: Scatter\n ax = axes[row, 2]\n step = max(1, len(f_c) // 2000)\n ax.scatter(o_c.values[::step], f_c.values[::step], s=3, alpha=0.3, color=\"steelblue\")\n lims = [min(o_c.min(), f_c.min()), max(o_c.max(), f_c.max())]\n ax.plot(lims, lims, \"k--\", alpha=0.5)\n ax.set_xlabel(f\"Obs ({m['unit']})\")\n ax.set_ylabel(f\"CW3E ({m['unit']})\")\n ax.set_title(f\"r={m['r']:.2f}\")\n ax.set_aspect(\"equal\")\n ax.grid(True, alpha=0.3)\n\n fig.suptitle(f\"{site_id} WY{water_year} — CW3E Forcing vs Ameriflux Observations\",\n fontsize=14, fontweight=\"bold\")\n plt.tight_layout()\n plt.show()\n\nclm_ds.close()" } ], "metadata": { diff --git a/helpers.py b/helpers.py index 190d2ff..9e8e3bf 100644 --- a/helpers.py +++ b/helpers.py @@ -77,6 +77,56 @@ # Combined lookup for any pf_indicator value CONUS2_SUBSURFACE = {**CONUS2_SOILS, **CONUS2_GEOLOGY} +# ============================================================================= +# SNOTEL Site Catalog (for snow/ notebooks) +# Subset of low-bias site-years from the snow_model_sensitivity study. +# triplet is the NRCS/HydroData site identifier (e.g. "679:WA:SNTL"). +# ============================================================================= +SNOTEL_SITES = { + # Cascades (maritime, heavy snowfall) + "Paradise": {"coords": (46.79, -121.74), "elev_m": 1564, "state": "WA", + "region": "Cascades", "triplet": "679:WA:SNTL"}, + "Stevens_Pass": {"coords": (47.74, -121.09), "elev_m": 1241, "state": "WA", + "region": "Cascades", "triplet": "791:WA:SNTL"}, + # Sierra Nevada + "CSS_Lab": {"coords": (39.31, -120.37), "elev_m": 2103, "state": "CA", + "region": "Sierra Nevada", "triplet": "428:CA:SNTL"}, + "Leavitt_Lake": {"coords": (38.28, -119.56), "elev_m": 2987, "state": "CA", + "region": "Sierra Nevada", "triplet": "518:CA:SNTL"}, + "Donner_Summit": {"coords": (39.32, -120.33), "elev_m": 2103, "state": "CA", + "region": "Sierra Nevada", "triplet": "428:CA:SNTL"}, + # Wyoming (continental) + "Togwotee_Pass": {"coords": (43.76, -110.07), "elev_m": 2936, "state": "WY", + "region": "Wyoming", "triplet": "822:WY:SNTL"}, + "Canyon": {"coords": (44.73, -110.50), "elev_m": 2438, "state": "WY", + "region": "Wyoming", "triplet": "384:WY:SNTL"}, + # Colorado Rockies + "Berthoud_Summit": {"coords": (39.80, -105.78), "elev_m": 3536, "state": "CO", + "region": "Colorado Rockies", "triplet": "335:CO:SNTL"}, + "Schofield_Pass": {"coords": (39.02, -107.05), "elev_m": 3261, "state": "CO", + "region": "Colorado Rockies", "triplet": "737:CO:SNTL"}, + # Wasatch / Utah + "Brighton": {"coords": (40.60, -111.58), "elev_m": 2667, "state": "UT", + "region": "Wasatch", "triplet": "366:UT:SNTL"}, + # Southern Rockies / arid + "Quemazon": {"coords": (35.93, -106.40), "elev_m": 2926, "state": "NM", + "region": "Southern Rockies", "triplet": "708:NM:SNTL"}, + # Northern Rockies + "Northeast_Entrance":{"coords": (45.00, -109.93), "elev_m": 2286, "state": "MT", + "region": "Northern Rockies", "triplet": "656:MT:SNTL"}, +} + +# Showcase subset — geographically diverse, recent CW3E-covered water years. +# Default WY choices come from the snow_model_sensitivity low-bias pool. +SNOTEL_SHOWCASE = [ + ("Paradise", 2022), + ("CSS_Lab", 2022), + ("Togwotee_Pass", 2020), + ("Berthoud_Summit", 2024), + ("Brighton", 2024), + ("Quemazon", 2023), +] + def write_site_config(filepath, params): """Write site parameters to a text file for passing between notebooks. diff --git a/snow/Single_Column_PFCLM_netcdf.ipynb b/snow/Single_Column_PFCLM_netcdf.ipynb new file mode 100644 index 0000000..ee36813 --- /dev/null +++ b/snow/Single_Column_PFCLM_netcdf.ipynb @@ -0,0 +1,501 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "snow_nb2_title", + "metadata": {}, + "source": [ + "# Single-Column ParFlow-CLM Simulation (snow example)\n", + "\n", + "Run a single-column PF-CLM simulation with the `clm_snow_improved` physics defaults:\n", + "- **Fractional snow cover**: `FracSnoScheme=\"CLM\"` with `FracSnoRoughness=1e-8`\n", + "- **Albedo decay**: slow aging (`AlbedoDecayVis=0.3`, `AlbedoDecayNir=0.1`)\n", + "- **Fresh snow albedo**: tuned high (`AlbedoVisNew=0.97`, `AlbedoNirNew=0.67`)\n", + "- **Shrubland land cover** (IGBP 7): snow-adapted canopy radiation per `snow_model_sensitivity`\n", + "\n", + "This notebook uses the forcing and CLM inputs prepared in Notebook 1.\n", + "Output is written as NetCDF for analysis in Notebook 3.\n", + "\n", + "**Requires:** ParFlow master (commit containing PR #695 \"Feature: CLM snow parameterization options\" merged 2026-01-28 and PR #709 \"SZA-modulated fSCA\" merged 2026-03-11) or later.\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "snow_nb2_imports", + "metadata": {}, + "outputs": [], + "source": [ + "import os\n", + "import sys\n", + "import numpy as np\n", + "import xarray as xr\n", + "import pandas as pd\n", + "import matplotlib.pyplot as plt\n", + "import matplotlib.dates as mdates\n", + "from pathlib import Path\n", + "from datetime import datetime\n", + "\n", + "sys.path.insert(0, os.path.abspath(\"..\"))\n", + "from helpers import read_site_config\n", + "\n", + "os.environ[\"PARFLOW_DIR\"] = os.environ.get(\n", + " \"PARFLOW_DIR\", str(Path.home() / \"parflow\" / \"parflow_21_Jan_26\"))\n", + "print(f\"PARFLOW_DIR: {os.environ['PARFLOW_DIR']}\")\n", + "\n", + "from parflow import Run\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "snow_nb2_loadcfg", + "metadata": {}, + "outputs": [], + "source": [ + "# ============================================================\n", + "# Load Site Configuration from Notebook 1\n", + "# ============================================================\n", + "static_write_dir = \"../runs/snow/Berthoud_Summit_WY2024\"\n", + "\n", + "run_dir = Path(static_write_dir)\n", + "cfg = read_site_config(run_dir / \"site_config.txt\")\n", + "\n", + "site_name = cfg[\"site_name\"]\n", + "water_year = int(cfg[\"water_year\"])\n", + "start = cfg[\"start\"]\n", + "end = cfg[\"end\"]\n", + "lat = float(cfg[\"lat\"])\n", + "lon = float(cfg[\"lon\"])\n", + "\n", + "# Soil (top 2m)\n", + "soil_ksat = float(cfg[\"soil_ksat\"])\n", + "soil_porosity = float(cfg[\"soil_porosity\"])\n", + "soil_vg_alpha = float(cfg[\"soil_vg_alpha\"])\n", + "soil_vg_n = float(cfg[\"soil_vg_n\"])\n", + "soil_sres = float(cfg[\"soil_sres\"])\n", + "\n", + "# Geology (bottom 6m)\n", + "geol_ksat = float(cfg[\"geol_ksat\"])\n", + "geol_porosity = float(cfg[\"geol_porosity\"])\n", + "geol_vg_alpha = float(cfg[\"geol_vg_alpha\"])\n", + "geol_vg_n = float(cfg[\"geol_vg_n\"])\n", + "geol_sres = float(cfg[\"geol_sres\"])\n", + "\n", + "wtd_m = float(cfg[\"wtd_m\"])\n", + "forcing_filename = cfg[\"forcing_file\"]\n", + "\n", + "# ---- Override any config values here if needed ----\n", + "# wtd_m = -8.0\n", + "\n", + "dt_start = datetime.strptime(start, \"%Y-%m-%d\")\n", + "dt_end = datetime.strptime(end, \"%Y-%m-%d\")\n", + "n_hours = int((dt_end - dt_start).total_seconds() / 3600)\n", + "\n", + "print(f\"Site: {site_name} WY{water_year}, {n_hours} hours\")\n", + "print(f\" Soil: Ksat={soil_ksat}, porosity={soil_porosity}\")\n", + "print(f\" Geology: Ksat={geol_ksat}, porosity={geol_porosity}\")\n", + "print(f\" WTD BC: {wtd_m} m\")\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "snow_nb2_setup", + "metadata": {}, + "outputs": [], + "source": [ + "# ============================================================\n", + "# ParFlow-CLM Run Setup (snow example)\n", + "# ============================================================\n", + "# Subsurface geometry is the same 8m / 10-layer column as the ET example.\n", + "# Snow-specific CLM keys are in the ----- CLM SNOW PHYSICS ----- section.\n", + "\n", + "run_name = \"pfclm_sc\"\n", + "run = Run(run_name)\n", + "run.FileVersion = 4\n", + "\n", + "# ----- TIMING -----\n", + "run.TimingInfo.BaseUnit = 1.0\n", + "run.TimingInfo.StartCount = 0\n", + "run.TimingInfo.StartTime = 0.0\n", + "run.TimingInfo.StopTime = n_hours\n", + "run.TimingInfo.DumpInterval = 1.0\n", + "run.TimeStep.Type = \"Constant\"\n", + "run.TimeStep.Value = 1.0\n", + "\n", + "# ----- GRID (single column, 10 layers, 8m deep) -----\n", + "run.Process.Topology.P = 1\n", + "run.Process.Topology.Q = 1\n", + "run.Process.Topology.R = 1\n", + "\n", + "run.ComputationalGrid.Lower.X = 0.0\n", + "run.ComputationalGrid.Lower.Y = 0.0\n", + "run.ComputationalGrid.Lower.Z = 0.0\n", + "run.ComputationalGrid.DX = 2.0\n", + "run.ComputationalGrid.DY = 2.0\n", + "run.ComputationalGrid.DZ = 1.0\n", + "run.ComputationalGrid.NX = 1\n", + "run.ComputationalGrid.NY = 1\n", + "run.ComputationalGrid.NZ = 10\n", + "\n", + "# ----- DOMAIN + TWO-LAYER GEOMETRY -----\n", + "run.GeomInput.Names = \"domain_input soil_input geol_input\"\n", + "\n", + "run.GeomInput.domain_input.InputType = \"Box\"\n", + "run.GeomInput.domain_input.GeomName = \"domain\"\n", + "run.Geom.domain.Lower.X = 0.0\n", + "run.Geom.domain.Lower.Y = 0.0\n", + "run.Geom.domain.Lower.Z = 0.0\n", + "run.Geom.domain.Upper.X = 2.0\n", + "run.Geom.domain.Upper.Y = 2.0\n", + "run.Geom.domain.Upper.Z = 10.0\n", + "run.Geom.domain.Patches = \"x_lower x_upper y_lower y_upper z_lower z_upper\"\n", + "\n", + "run.GeomInput.soil_input.InputType = \"Box\"\n", + "run.GeomInput.soil_input.GeomName = \"soil\"\n", + "run.Geom.soil.Lower.X = 0.0\n", + "run.Geom.soil.Lower.Y = 0.0\n", + "run.Geom.soil.Lower.Z = 6.0\n", + "run.Geom.soil.Upper.X = 2.0\n", + "run.Geom.soil.Upper.Y = 2.0\n", + "run.Geom.soil.Upper.Z = 10.0\n", + "\n", + "run.GeomInput.geol_input.InputType = \"Box\"\n", + "run.GeomInput.geol_input.GeomName = \"geology\"\n", + "run.Geom.geology.Lower.X = 0.0\n", + "run.Geom.geology.Lower.Y = 0.0\n", + "run.Geom.geology.Lower.Z = 0.0\n", + "run.Geom.geology.Upper.X = 2.0\n", + "run.Geom.geology.Upper.Y = 2.0\n", + "run.Geom.geology.Upper.Z = 6.0\n", + "\n", + "# ----- VARIABLE DZ -----\n", + "run.Solver.Nonlinear.VariableDz = True\n", + "run.dzScale.GeomNames = \"domain\"\n", + "run.dzScale.Type = \"nzList\"\n", + "run.dzScale.nzListNumber = 10\n", + "run.Cell._0.dzScale.Value = 1.0\n", + "run.Cell._1.dzScale.Value = 1.0\n", + "run.Cell._2.dzScale.Value = 1.0\n", + "run.Cell._3.dzScale.Value = 1.0\n", + "run.Cell._4.dzScale.Value = 1.0\n", + "run.Cell._5.dzScale.Value = 1.0\n", + "run.Cell._6.dzScale.Value = 1.0\n", + "run.Cell._7.dzScale.Value = 0.6\n", + "run.Cell._8.dzScale.Value = 0.3\n", + "run.Cell._9.dzScale.Value = 0.1\n", + "\n", + "# ----- HYDRAULIC PROPERTIES (two-layer) -----\n", + "run.Geom.Perm.Names = \"soil geology\"\n", + "run.Geom.soil.Perm.Type = \"Constant\"\n", + "run.Geom.soil.Perm.Value = soil_ksat\n", + "run.Geom.geology.Perm.Type = \"Constant\"\n", + "run.Geom.geology.Perm.Value = geol_ksat\n", + "\n", + "run.Perm.TensorType = \"TensorByGeom\"\n", + "run.Geom.Perm.TensorByGeom.Names = \"domain\"\n", + "run.Geom.domain.Perm.TensorValX = 1.0\n", + "run.Geom.domain.Perm.TensorValY = 1.0\n", + "run.Geom.domain.Perm.TensorValZ = 1.0\n", + "\n", + "run.SpecificStorage.Type = \"Constant\"\n", + "run.SpecificStorage.GeomNames = \"domain\"\n", + "run.Geom.domain.SpecificStorage.Value = 1.0e-4\n", + "\n", + "run.Geom.Porosity.GeomNames = \"soil geology\"\n", + "run.Geom.soil.Porosity.Type = \"Constant\"\n", + "run.Geom.soil.Porosity.Value = soil_porosity\n", + "run.Geom.geology.Porosity.Type = \"Constant\"\n", + "run.Geom.geology.Porosity.Value = geol_porosity\n", + "\n", + "# ----- PHASES -----\n", + "run.Phase.Names = \"water\"\n", + "run.Phase.water.Density.Type = \"Constant\"\n", + "run.Phase.water.Density.Value = 1.0\n", + "run.Phase.water.Viscosity.Type = \"Constant\"\n", + "run.Phase.water.Viscosity.Value = 1.0\n", + "run.Phase.water.Mobility.Type = \"Constant\"\n", + "run.Phase.water.Mobility.Value = 1.0\n", + "\n", + "run.Phase.RelPerm.Type = \"VanGenuchten\"\n", + "run.Phase.RelPerm.GeomNames = \"soil geology\"\n", + "run.Geom.soil.RelPerm.Alpha = soil_vg_alpha\n", + "run.Geom.soil.RelPerm.N = soil_vg_n\n", + "run.Geom.geology.RelPerm.Alpha = geol_vg_alpha\n", + "run.Geom.geology.RelPerm.N = geol_vg_n\n", + "\n", + "run.Phase.Saturation.Type = \"VanGenuchten\"\n", + "run.Phase.Saturation.GeomNames = \"soil geology\"\n", + "run.Geom.soil.Saturation.Alpha = soil_vg_alpha\n", + "run.Geom.soil.Saturation.N = soil_vg_n\n", + "run.Geom.soil.Saturation.SRes = soil_sres\n", + "run.Geom.soil.Saturation.SSat = 1.0\n", + "run.Geom.geology.Saturation.Alpha = geol_vg_alpha\n", + "run.Geom.geology.Saturation.N = geol_vg_n\n", + "run.Geom.geology.Saturation.SRes = geol_sres\n", + "run.Geom.geology.Saturation.SSat = 1.0\n", + "\n", + "# ----- MISC -----\n", + "run.Contaminants.Names = \"\"\n", + "run.Gravity = 1.0\n", + "run.Domain.GeomName = \"domain\"\n", + "run.Wells.Names = \"\"\n", + "run.KnownSolution = \"NoKnownSolution\"\n", + "\n", + "# ----- TIME CYCLES -----\n", + "run.Cycle.Names = \"constant\"\n", + "run.Cycle.constant.Names = \"alltime\"\n", + "run.Cycle.constant.alltime.Length = 1\n", + "run.Cycle.constant.Repeat = -1\n", + "\n", + "# ----- BOUNDARY CONDITIONS -----\n", + "run.BCPressure.PatchNames = \"x_lower x_upper y_lower y_upper z_lower z_upper\"\n", + "\n", + "run.Patch.x_lower.BCPressure.Type = \"FluxConst\"\n", + "run.Patch.x_lower.BCPressure.Cycle = \"constant\"\n", + "run.Patch.x_lower.BCPressure.alltime.Value = 0.0\n", + "run.Patch.y_lower.BCPressure.Type = \"FluxConst\"\n", + "run.Patch.y_lower.BCPressure.Cycle = \"constant\"\n", + "run.Patch.y_lower.BCPressure.alltime.Value = 0.0\n", + "\n", + "run.Patch.z_lower.BCPressure.Type = \"DirEquilRefPatch\"\n", + "run.Patch.z_lower.BCPressure.RefGeom = \"domain\"\n", + "run.Patch.z_lower.BCPressure.RefPatch = \"z_upper\"\n", + "run.Patch.z_lower.BCPressure.Cycle = \"constant\"\n", + "run.Patch.z_lower.BCPressure.alltime.Value = wtd_m\n", + "\n", + "run.Patch.x_upper.BCPressure.Type = \"FluxConst\"\n", + "run.Patch.x_upper.BCPressure.Cycle = \"constant\"\n", + "run.Patch.x_upper.BCPressure.alltime.Value = 0.0\n", + "run.Patch.y_upper.BCPressure.Type = \"FluxConst\"\n", + "run.Patch.y_upper.BCPressure.Cycle = \"constant\"\n", + "run.Patch.y_upper.BCPressure.alltime.Value = 0.0\n", + "run.Patch.z_upper.BCPressure.Type = \"OverlandFlow\"\n", + "run.Patch.z_upper.BCPressure.Cycle = \"constant\"\n", + "run.Patch.z_upper.BCPressure.alltime.Value = 0.0\n", + "\n", + "# ----- SLOPES -----\n", + "run.TopoSlopesX.Type = \"Constant\"\n", + "run.TopoSlopesX.GeomNames = \"domain\"\n", + "run.TopoSlopesX.Geom.domain.Value = 0.1\n", + "run.TopoSlopesY.Type = \"Constant\"\n", + "run.TopoSlopesY.GeomNames = \"domain\"\n", + "run.TopoSlopesY.Geom.domain.Value = 0.0\n", + "run.Mannings.Type = \"Constant\"\n", + "run.Mannings.GeomNames = \"domain\"\n", + "run.Mannings.Geom.domain.Value = 2.e-6\n", + "run.PhaseSources.water.Type = \"Constant\"\n", + "run.PhaseSources.water.GeomNames = \"domain\"\n", + "run.PhaseSources.water.Geom.domain.Value = 0.0\n", + "\n", + "# ----- SOLVER -----\n", + "run.Solver = \"Richards\"\n", + "run.Solver.MaxIter = 9000\n", + "run.Solver.Nonlinear.MaxIter = 100\n", + "run.Solver.Nonlinear.ResidualTol = 1e-7\n", + "run.Solver.Nonlinear.EtaChoice = \"EtaConstant\"\n", + "run.Solver.Nonlinear.EtaValue = 0.01\n", + "run.Solver.Nonlinear.UseJacobian = True\n", + "run.Solver.Nonlinear.DerivativeEpsilon = 1e-12\n", + "run.Solver.Nonlinear.StepTol = 1e-15\n", + "run.Solver.Nonlinear.Globalization = \"LineSearch\"\n", + "run.Solver.Linear.KrylovDimension = 100\n", + "run.Solver.Linear.MaxRestarts = 5\n", + "run.Solver.Linear.Preconditioner = \"PFMG\"\n", + "run.Solver.PrintSubsurf = False\n", + "run.Solver.Drop = 1E-20\n", + "run.Solver.AbsTol = 1E-9\n", + "\n", + "# ----- OUTPUT: PFB (disabled) -----\n", + "run.Solver.PrintSubsurfData = False\n", + "run.Solver.PrintPressure = False\n", + "run.Solver.PrintSaturation = False\n", + "run.Solver.PrintCLM = False\n", + "run.Solver.PrintMask = False\n", + "run.Solver.PrintSpecificStorage = False\n", + "run.Solver.PrintEvapTrans = False\n", + "run.Solver.WriteSiloMannings = False\n", + "run.Solver.WriteSiloMask = False\n", + "run.Solver.WriteSiloSlopes = False\n", + "run.Solver.WriteSiloSaturation = False\n", + "\n", + "# ----- OUTPUT: NetCDF -----\n", + "run.NetCDF.NumStepsPerFile = n_hours\n", + "run.NetCDF.WritePressure = True\n", + "run.NetCDF.WriteSubsurface = False\n", + "run.NetCDF.WriteSaturation = True\n", + "run.NetCDF.WriteCLM = True\n", + "run.NetCDF.CLMNumStepsPerFile = n_hours\n", + "\n", + "# ----- CLM CORE CONFIGURATION -----\n", + "run.Solver.LSM = \"CLM\"\n", + "run.Solver.CLM.MetForcing = \"1D\"\n", + "run.Solver.CLM.MetFileName = forcing_filename\n", + "run.Solver.CLM.MetFilePath = \".\"\n", + "\n", + "run.Solver.CLM.EvapBeta = \"Linear\"\n", + "run.Solver.CLM.VegWaterStress = \"Saturation\"\n", + "run.Solver.CLM.ResSat = 0.2\n", + "run.Solver.CLM.WiltingPoint = 0.2\n", + "run.Solver.CLM.FieldCapacity = 1.0\n", + "run.Solver.CLM.IrrigationType = \"none\"\n", + "run.Solver.CLM.RZWaterStress = 1\n", + "run.Solver.CLM.RootZoneNZ = 5\n", + "run.Solver.CLM.SoiLayer = 4\n", + "\n", + "# ----- CLM SNOW PHYSICS (clm_snow_improved) -----\n", + "# Fractional snow cover: CLM default scheme with minimal roughness\n", + "run.Solver.CLM.FracSnoScheme = \"CLM\"\n", + "run.Solver.CLM.FracSnoRoughness = 1e-8\n", + "\n", + "# Albedo \u2014 slow aging + tuned fresh-snow values (best low-bias result:\n", + "# |bias|=13.7%, melt-out +7d, NSE=0.87 in the shrub+combined study)\n", + "run.Solver.CLM.AlbedoDecayVis = 0.3\n", + "run.Solver.CLM.AlbedoDecayNir = 0.1\n", + "run.Solver.CLM.AlbedoVisNew = 0.97\n", + "run.Solver.CLM.AlbedoNirNew = 0.67\n", + "\n", + "# --- Optional: switch to SZA-modulated fSCA (PR #709) ---\n", + "# Uncomment to replace FracSnoScheme=CLM with the solar-zenith-angle\n", + "# scheme. This is the \"interpolating\" formulation that transitions\n", + "# between bare-soil and snow roughness based on sun angle.\n", + "# run.Solver.CLM.FracSnoScheme = \"SZA\"\n", + "# run.Solver.CLM.FracSnoRoughnessMin = 1e-8\n", + "# run.Solver.CLM.FracSnoRoughnessMax = 0.2\n", + "# run.Solver.CLM.FracSnoGammaSZA = 4.0\n", + "# run.Solver.CLM.FracSnoAvgWindow = 72.0\n", + "\n", + "# ----- CLM OUTPUT -----\n", + "run.Solver.PrintLSMSink = False\n", + "run.Solver.CLM.CLMDumpInterval = 1\n", + "run.Solver.CLM.CLMFileDir = \"output/\"\n", + "run.Solver.CLM.BinaryOutDir = False\n", + "run.Solver.CLM.IstepStart = 1\n", + "run.Solver.WriteCLMBinary = False\n", + "run.Solver.WriteSiloCLM = False\n", + "run.Solver.CLM.WriteLogs = False\n", + "run.Solver.CLM.WriteLastRST = True\n", + "run.Solver.CLM.DailyRST = False\n", + "run.Solver.CLM.SingleFile = True\n", + "\n", + "# ----- INITIAL CONDITIONS -----\n", + "run.ICPressure.Type = \"HydroStaticPatch\"\n", + "run.ICPressure.GeomNames = \"domain\"\n", + "run.Geom.domain.ICPressure.Value = wtd_m\n", + "run.Geom.domain.ICPressure.RefGeom = \"domain\"\n", + "run.Geom.domain.ICPressure.RefPatch = \"z_lower\"\n", + "\n", + "print(\"ParFlow run configured (snow):\")\n", + "print(f\" Subsurface: soil Ksat={soil_ksat}, geol Ksat={geol_ksat}, WTD={wtd_m} m\")\n", + "print(f\" Snow: FracSno=CLM (z0=1e-8), Albedo slow-decay + high-fresh\")\n", + "print(f\" Land cover: IGBP 7 (Open Shrubland)\")\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "snow_nb2_run", + "metadata": {}, + "outputs": [], + "source": [ + "# Run ParFlow\n", + "run.run(working_directory=static_write_dir, skip_validation=True)\n" + ] + }, + { + "cell_type": "markdown", + "id": "snow_nb2_verify_md", + "metadata": {}, + "source": [ + "## Quick-Look Output Verification\n", + "\n", + "Load the NetCDF output and plot SWE + snow depth to verify the run completed successfully.\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "snow_nb2_load", + "metadata": {}, + "outputs": [], + "source": [ + "# Load CLM output\n", + "run_dir = Path(static_write_dir)\n", + "clm_files = sorted(run_dir.glob(\"*.out.CLM.*.nc\"))\n", + "\n", + "if not clm_files:\n", + " raise FileNotFoundError(f\"No CLM NetCDF output found in {run_dir}\")\n", + "\n", + "clm_ds = xr.open_dataset(clm_files[0])\n", + "n_steps = clm_ds.sizes.get(\"time\", 0)\n", + "\n", + "if n_steps < n_hours:\n", + " print(f\"WARNING: Run incomplete! Got {n_steps}/{n_hours} timesteps \"\n", + " f\"({n_steps/n_hours*100:.0f}%). Check pfclm_sc.out.log for solver failures.\")\n", + "else:\n", + " print(f\"Run complete: {n_steps} timesteps\")\n", + "\n", + "print(f\"CLM variables: {list(clm_ds.data_vars)}\")\n", + "if \"swe_out\" in clm_ds:\n", + " swe = clm_ds[\"swe_out\"].values.flatten()\n", + " print(f\" SWE range: {swe.min():.1f} to {swe.max():.1f} mm, peak={swe.max():.1f} mm\")\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "snow_nb2_plot", + "metadata": {}, + "outputs": [], + "source": [ + "# Plot SWE and snow depth\n", + "dates = pd.date_range(start, periods=n_steps, freq=\"h\")\n", + "\n", + "fig, axes = plt.subplots(2, 1, figsize=(14, 7), sharex=True)\n", + "\n", + "if \"swe_out\" in clm_ds:\n", + " axes[0].plot(dates, clm_ds[\"swe_out\"].values.flatten(),\n", + " lw=0.5, color=\"steelblue\")\n", + " axes[0].set_ylabel(\"SWE [mm]\")\n", + " axes[0].set_title(f\"{site_name} WY{water_year} \u2014 PF-CLM Snow Output \"\n", + " f\"({n_steps}/{n_hours} hrs)\")\n", + " axes[0].grid(alpha=0.3)\n", + "\n", + "# Snow depth (snowdp) if available\n", + "snowdp_var = None\n", + "for v in [\"snowdp\", \"SnowDepth\", \"snow_depth\"]:\n", + " if v in clm_ds:\n", + " snowdp_var = v\n", + " break\n", + "if snowdp_var is not None:\n", + " axes[1].plot(dates, clm_ds[snowdp_var].values.flatten(),\n", + " lw=0.5, color=\"navy\")\n", + " axes[1].set_ylabel(f\"Snow depth [m] ({snowdp_var})\")\n", + "else:\n", + " axes[1].text(0.5, 0.5, \"No snow depth variable in output\",\n", + " transform=axes[1].transAxes, ha=\"center\", va=\"center\")\n", + "\n", + "axes[1].set_xlabel(\"Date\")\n", + "axes[1].xaxis.set_major_formatter(mdates.DateFormatter(\"%b %Y\"))\n", + "axes[1].grid(alpha=0.3)\n", + "plt.tight_layout()\n", + "plt.show()\n", + "\n", + "clm_ds.close()\n" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + }, + "language_info": { + "name": "python" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} \ No newline at end of file diff --git a/snow/locate_plot_station_get_forcing.ipynb b/snow/locate_plot_station_get_forcing.ipynb new file mode 100644 index 0000000..4fcaa17 --- /dev/null +++ b/snow/locate_plot_station_get_forcing.ipynb @@ -0,0 +1,467 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "snow_nb1_title", + "metadata": {}, + "source": [ + "# Locate SNOTEL Station & Download Forcing (snow example)\n", + "\n", + "This notebook helps you:\n", + "1. Pick a SNOTEL snow-pillow site from the showcase catalog (or enter your own)\n", + "2. Query CONUS2 subsurface parameters at the site location\n", + "3. Prepare CLM input files (vegm, clmin)\n", + "4. Download CW3E meteorological forcing\n", + "5. Preview the SNOTEL SWE observations the comparison notebook will use\n", + "\n", + "**Prerequisites:** HydroData account (free) \u2014 see README.md for signup instructions.\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "snow_nb1_imports", + "metadata": {}, + "outputs": [], + "source": [ + "# Imports\n", + "import os\n", + "import sys\n", + "import numpy as np\n", + "import pandas as pd\n", + "import matplotlib.pyplot as plt\n", + "import matplotlib.dates as mdates\n", + "import subsettools as st\n", + "import hf_hydrodata as hf\n", + "import shutil\n", + "from pathlib import Path\n", + "from datetime import datetime\n", + "\n", + "# Make helpers.py (at repo root) importable from this subdir\n", + "sys.path.insert(0, os.path.abspath(\"..\"))\n", + "from helpers import (IGBP_NAMES, IGBP_FULL_NAMES,\n", + " CONUS2_SOILS, CONUS2_GEOLOGY, CONUS2_SUBSURFACE,\n", + " SNOTEL_SITES, SNOTEL_SHOWCASE,\n", + " build_vegm, patch_clmin_dates, write_site_config)\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "snow_nb1_userconfig", + "metadata": {}, + "outputs": [], + "source": [ + "# ============================================================\n", + "# User Configuration\n", + "# ============================================================\n", + "# Pick a SNOTEL site + water year. Defaults to a showcase entry.\n", + "\n", + "site_name = \"Berthoud_Summit\" # key into SNOTEL_SITES (see below)\n", + "water_year = 2024 # Oct 1 (WY-1) to Oct 1 (WY)\n", + "\n", + "start = f\"{water_year - 1}-10-01\"\n", + "end = f\"{water_year}-10-01\"\n", + "\n", + "# Directory where forcing and CLM input files will be written\n", + "# (path is relative to this notebook's directory, i.e. snow/)\n", + "static_write_dir = f\"../runs/snow/{site_name}_WY{water_year}\"\n", + "\n", + "# Forcing filename\n", + "forcing_filename = f\"forcing1D.{site_name}.{start}-{end}.txt\"\n", + "\n", + "# HydroData PIN (only needed once per machine)\n", + "# Sign up: https://hydrogen.princeton.edu/signup\n", + "# Get PIN: https://hydrogen.princeton.edu/pin\n", + "hf.register_api_pin(\"your_email@example.com\", \"your_pin\")\n", + "\n", + "# -----------------------------------------------------------\n", + "# Showcase sites (geographically diverse, recent CW3E-covered years)\n", + "# -----------------------------------------------------------\n", + "print(\"Showcase site-years:\")\n", + "for s, wy in SNOTEL_SHOWCASE:\n", + " info = SNOTEL_SITES[s]\n", + " print(f\" {s:<18s} WY{wy} {info['region']:<18s} {info['elev_m']} m triplet={info['triplet']}\")\n", + "\n", + "# -----------------------------------------------------------\n", + "# Full SNOTEL catalog: see helpers.SNOTEL_SITES for the 12 curated sites.\n", + "# To extend, add entries following the same schema \u2014 coords, elev_m, state,\n", + "# region, triplet. The low-bias site-years from snow_model_sensitivity\n", + "# (|T bias|<1C, |P bias|<10%) are documented in the repo README.\n", + "# -----------------------------------------------------------\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "snow_nb1_resolve", + "metadata": {}, + "outputs": [], + "source": [ + "# Resolve site metadata\n", + "if site_name not in SNOTEL_SITES:\n", + " raise KeyError(f\"{site_name} not in SNOTEL_SITES \u2014 add it to helpers.SNOTEL_SITES\")\n", + "\n", + "site = SNOTEL_SITES[site_name]\n", + "lat, lon = site[\"coords\"]\n", + "elev_m = site[\"elev_m\"]\n", + "region = site[\"region\"]\n", + "triplet = site[\"triplet\"]\n", + "\n", + "print(f\"Site: {site_name} ({region}, {site['state']})\")\n", + "print(f\" lat={lat:.4f}, lon={lon:.4f}, elev={elev_m} m\")\n", + "print(f\" SNOTEL triplet: {triplet}\")\n" + ] + }, + { + "cell_type": "markdown", + "id": "snow_nb1_subsurface_md", + "metadata": {}, + "source": [ + "## CONUS2 Grid & Subsurface Parameters\n", + "\n", + "Snow ET sensitivity is not especially sensitive to subsurface \u2014 but we still\n", + "need sensible soil/geology for the ParFlow domain. Query CONUS2 to populate\n", + "them; water table depth defaults to deep (well below the 8m domain) since\n", + "SNOTEL sites are mountain ridges where deep water tables are the norm.\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "snow_nb1_cg_subsurf", + "metadata": {}, + "outputs": [], + "source": [ + "# Define single-column CONUS2 domain (single point = single CONUS2 cell)\n", + "latlon_bounds = [[lat, lon], [lat, lon]]\n", + "ij_bounds, mask = st.define_latlon_domain(latlon_bounds=latlon_bounds, grid=\"conus2\")\n", + "print(f\"CONUS2 grid bounds (i_start, j_start, i_end, j_end): {ij_bounds}\")\n", + "\n", + "# Query full pf_indicator profile (10 layers, bottom to top)\n", + "ind = hf.get_gridded_data({\n", + " \"dataset\": \"conus2_domain\", \"variable\": \"pf_indicator\",\n", + " \"grid_bounds\": ij_bounds,\n", + "})\n", + "print(f\"\\nCONUS2 subsurface profile (bottom -> top):\")\n", + "for i in range(ind.shape[0]):\n", + " t = int(ind[i, 0, 0])\n", + " info = CONUS2_SUBSURFACE.get(t, (\"Unknown\",))\n", + " print(f\" Layer {i}: type {t} ({info[0]})\")\n", + "\n", + "# Soil = top layer, Geology = layer 4 (from bottom)\n", + "soil_type = int(ind[-1, 0, 0])\n", + "geol_type = int(ind[4, 0, 0])\n", + "\n", + "soil_info = CONUS2_SOILS.get(soil_type)\n", + "geol_info = CONUS2_GEOLOGY.get(geol_type)\n", + "\n", + "print(f\"\\n--- Soil (top layer, type {soil_type}: {soil_info[0]}) ---\")\n", + "print(f\" Ksat={soil_info[1]} m/h, porosity={soil_info[2]}, \"\n", + " f\"VG alpha={soil_info[3]}, VG n={soil_info[4]}, Sres={soil_info[5]}\")\n", + "if geol_info:\n", + " print(f\"\\n--- Geology (layer 4, type {geol_type}: {geol_info[0]}) ---\")\n", + " print(f\" Ksat={geol_info[1]} m/h, porosity={geol_info[2]}, \"\n", + " f\"VG alpha={geol_info[3]}, VG n={geol_info[4]}, Sres={geol_info[5]}\")\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "snow_nb1_wtd", + "metadata": {}, + "outputs": [], + "source": [ + "## Water table depth \u2014 default deep for mountain snow sites.\n", + "## Override here if you have site-specific WTD (negative = below surface).\n", + "wtd_m = -10.0 # m \u2014 below the 8m domain; snow is rainfall-dominated at depth\n", + "print(f\"WTD BC: {wtd_m} m\")\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "snow_nb1_vegm", + "metadata": {}, + "outputs": [], + "source": [ + "# ============================================================\n", + "# Vegm parameters \u2014 forced to IGBP 7 (shrubland) for snow default config\n", + "# ============================================================\n", + "# The snow sensitivity study found shrubland IGBP produces the best SWE\n", + "# timing across sites. Override if a site has a known different land cover.\n", + "igbp_type = 7 # Open Shrubland \u2014 best snow config default\n", + "\n", + "# Sand/clay: try CONUS2 config_clm, fall back to defaults\n", + "import tempfile\n", + "_tmp_dir = tempfile.mkdtemp(prefix=\"clm_tmp_\")\n", + "try:\n", + " _ = st.config_clm(ij_bounds, start=start, end=end,\n", + " dataset=\"conus2_domain\", write_dir=_tmp_dir)\n", + " with open(os.path.join(_tmp_dir, \"drv_vegm.dat\")) as _f:\n", + " _lines = _f.readlines()\n", + " for _line in _lines:\n", + " _parts = _line.split()\n", + " if len(_parts) >= 25:\n", + " sand_frac = float(_parts[4])\n", + " clay_frac = float(_parts[5])\n", + " color_index = int(_parts[6])\n", + " break\n", + " print(\"Sand/clay from CONUS2 domain\")\n", + "except Exception as e:\n", + " print(f\"config_clm failed ({e}), using defaults\")\n", + " sand_frac, clay_frac, color_index = 0.40, 0.20, 4\n", + "\n", + "print(f\"\\n--- Vegetation / Soil Texture ---\")\n", + "print(f\" IGBP: {igbp_type} \u2014 {IGBP_FULL_NAMES[igbp_type]} ({IGBP_NAMES[igbp_type]})\")\n", + "print(f\" Sand: {sand_frac:.0%}\")\n", + "print(f\" Clay: {clay_frac:.0%}\")\n", + "print(f\" Soil color: {color_index}\")\n" + ] + }, + { + "cell_type": "markdown", + "id": "snow_nb1_prep_md", + "metadata": {}, + "source": [ + "## Prepare Run Directory\n", + "\n", + "Write CLM input files (drv_vegm.dat, drv_clmin.dat, drv_vegp.dat) to the run directory.\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "snow_nb1_prep", + "metadata": {}, + "outputs": [], + "source": [ + "# Create run directory and write CLM input files\n", + "run_dir = Path(static_write_dir)\n", + "run_dir.mkdir(parents=True, exist_ok=True)\n", + "\n", + "# 1. drv_vegm.dat (shrubland = IGBP 7)\n", + "build_vegm(run_dir / \"drv_vegm.dat\", lat, lon, sand_frac, clay_frac, color_index, igbp_type)\n", + "print(f\"Wrote drv_vegm.dat (IGBP={igbp_type}, shrubland)\")\n", + "\n", + "# 2. drv_clmin.dat (shared template, patched dates)\n", + "clm_inputs_shared = Path(\"../clm_inputs\")\n", + "patch_clmin_dates(clm_inputs_shared / \"drv_clmin_template.dat\",\n", + " run_dir / \"drv_clmin.dat\", water_year)\n", + "print(f\"Wrote drv_clmin.dat (WY{water_year})\")\n", + "\n", + "# 3. drv_vegp.dat (snow-specific: copied from snow_model_sensitivity reference)\n", + "clm_inputs_snow = Path(\"../clm_inputs/snow\")\n", + "shutil.copy2(clm_inputs_snow / \"drv_vegp.dat\", run_dir / \"drv_vegp.dat\")\n", + "print(f\"Copied drv_vegp.dat (clm_snow_improved config)\")\n", + "\n", + "print(f\"\\nRun directory: {run_dir}\")\n", + "print(\"Contents:\", [f.name for f in sorted(run_dir.iterdir())])\n" + ] + }, + { + "cell_type": "markdown", + "id": "snow_nb1_forcing_md", + "metadata": {}, + "source": [ + "## Download CW3E Forcing\n", + "\n", + "Download 8 meteorological forcing variables from CW3E via HydroData.\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "snow_nb1_download", + "metadata": {}, + "outputs": [], + "source": [ + "# CW3E forcing in ParFlow-CLM 1D order:\n", + "# DSWR, DLWR, APCP [mm/s], Temp [K], UGRD, VGRD, Press [Pa], SPFH [kg/kg]\n", + "forcing_vars = [\n", + " \"downward_shortwave\", \"downward_longwave\", \"precipitation\",\n", + " \"air_temp\", \"east_windspeed\", \"north_windspeed\",\n", + " \"atmospheric_pressure\", \"specific_humidity\",\n", + "]\n", + "\n", + "dt_start = datetime.strptime(start, \"%Y-%m-%d\")\n", + "dt_end = datetime.strptime(end, \"%Y-%m-%d\")\n", + "n_hours = int((dt_end - dt_start).total_seconds() / 3600)\n", + "print(f\"Downloading {n_hours} hours for {site_name}...\")\n", + "\n", + "forcing_data = np.zeros((n_hours, 8))\n", + "for i, var in enumerate(forcing_vars):\n", + " print(f\" [{i+1}/8] {var}...\", end=\" \", flush=True)\n", + " options = {\n", + " \"dataset\": \"CW3E\", \"period\": \"hourly\", \"variable\": var,\n", + " \"start_time\": start, \"end_time\": end,\n", + " \"grid_bounds\": ij_bounds, \"dataset_version\": \"1.0\",\n", + " }\n", + " data = hf.get_gridded_data(options)\n", + " col = data[:, 0, 0] if data.ndim == 3 else data.flatten()\n", + " if len(col) != n_hours:\n", + " col = col[:n_hours] if len(col) > n_hours else np.pad(col, (0, n_hours - len(col)))\n", + " forcing_data[:, i] = col\n", + " print(f\"OK (range: {col.min():.3g} to {col.max():.3g})\")\n", + "\n", + "forcing_path = run_dir / forcing_filename\n", + "np.savetxt(forcing_path, forcing_data, fmt=\"%.6e\", delimiter=\" \")\n", + "print(f\"\\nSaved: {forcing_path} ({forcing_data.shape[0]} hours, \"\n", + " f\"{forcing_path.stat().st_size/1024:.0f} KB)\")\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "snow_nb1_plot_forcing", + "metadata": {}, + "outputs": [], + "source": [ + "# Plot forcing: precipitation and temperature\n", + "forcing_data = np.loadtxt(run_dir / forcing_filename)\n", + "n_hours = forcing_data.shape[0]\n", + "dates = pd.date_range(start, periods=n_hours, freq=\"h\")\n", + "\n", + "fig, axes = plt.subplots(2, 1, figsize=(14, 6), sharex=True)\n", + "axes[0].plot(dates, forcing_data[:, 2] * 3600, color=\"steelblue\", lw=0.5)\n", + "axes[0].set_ylabel(\"Precip [mm/h]\")\n", + "axes[0].set_title(f\"{site_name} CW3E Forcing (WY{water_year})\")\n", + "\n", + "axes[1].plot(dates, forcing_data[:, 3] - 273.15, color=\"firebrick\", lw=0.3)\n", + "axes[1].axhline(0, color=\"gray\", lw=0.5, ls=\"--\")\n", + "axes[1].set_ylabel(\"Air Temp [C]\")\n", + "axes[1].set_xlabel(\"Date\")\n", + "axes[1].xaxis.set_major_formatter(mdates.DateFormatter(\"%b %Y\"))\n", + "plt.tight_layout()\n", + "plt.show()\n" + ] + }, + { + "cell_type": "markdown", + "id": "snow_nb1_swe_md", + "metadata": {}, + "source": [ + "## Preview SNOTEL SWE Observations\n", + "\n", + "Fetch the daily SWE time series the comparison notebook will use.\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "snow_nb1_swe_plot", + "metadata": {}, + "outputs": [], + "source": [ + "# Fetch SNOTEL SWE for preview\n", + "try:\n", + " df = hf.get_point_data(\n", + " dataset=\"snotel\", variable=\"swe\",\n", + " temporal_resolution=\"daily\", aggregation=\"sod\",\n", + " site_ids=[triplet],\n", + " date_start=start, date_end=end,\n", + " )\n", + " if df is None or df.empty or triplet not in df.columns:\n", + " print(f\"No SNOTEL SWE data returned for triplet {triplet}\")\n", + " obs_swe = None\n", + " else:\n", + " obs_dates = pd.to_datetime(df[\"date\"])\n", + " obs_swe = pd.Series(df[triplet].values.astype(float),\n", + " index=pd.DatetimeIndex(obs_dates), name=\"swe_mm\").dropna()\n", + " print(f\"SNOTEL SWE: {len(obs_swe)} daily values, \"\n", + " f\"peak={obs_swe.max():.0f} mm, mean={obs_swe.mean():.0f} mm\")\n", + "except Exception as e:\n", + " print(f\"SNOTEL fetch failed: {e}\")\n", + " obs_swe = None\n", + "\n", + "if obs_swe is not None and len(obs_swe) > 0:\n", + " fig, ax = plt.subplots(figsize=(14, 4))\n", + " ax.plot(obs_swe.index, obs_swe.values, color=\"steelblue\", lw=1)\n", + " ax.set_ylabel(\"SWE [mm]\")\n", + " ax.set_title(f\"{site_name} WY{water_year} \u2014 SNOTEL SWE observations\")\n", + " ax.xaxis.set_major_formatter(mdates.DateFormatter(\"%b %Y\"))\n", + " ax.grid(alpha=0.3)\n", + " plt.tight_layout()\n", + " plt.show()\n" + ] + }, + { + "cell_type": "markdown", + "id": "snow_nb1_config_md", + "metadata": {}, + "source": [ + "## Write Site Configuration\n", + "\n", + "Save all site parameters to a config file that Notebooks 2 and 3 will load.\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "snow_nb1_writecfg", + "metadata": {}, + "outputs": [], + "source": [ + "# Write site_config.txt\n", + "config_path = run_dir / \"site_config.txt\"\n", + "s = CONUS2_SOILS[soil_type]\n", + "params = {\n", + " \"site_name\": site_name,\n", + " \"water_year\": water_year,\n", + " \"start\": start,\n", + " \"end\": end,\n", + " \"lat\": f\"{lat:.4f}\",\n", + " \"lon\": f\"{lon:.4f}\",\n", + " \"elev_m\": elev_m,\n", + " \"region\": region,\n", + " \"triplet\": triplet,\n", + " \"igbp\": igbp_type,\n", + " \"utc_offset\": round(lon / 15),\n", + " \"_comment_soil\": f\"Soil (CONUS2 top layer, type {soil_type}: {s[0]})\",\n", + " \"soil_ksat\": s[1],\n", + " \"soil_porosity\": s[2],\n", + " \"soil_vg_alpha\": s[3],\n", + " \"soil_vg_n\": s[4],\n", + " \"soil_sres\": s[5],\n", + "}\n", + "if geol_info:\n", + " g = geol_info\n", + " params[\"_comment_geol\"] = f\"Geology (CONUS2 layer 4, type {geol_type}: {g[0]})\"\n", + " params[\"geol_ksat\"] = g[1]\n", + " params[\"geol_porosity\"] = g[2]\n", + " params[\"geol_vg_alpha\"] = g[3]\n", + " params[\"geol_vg_n\"] = g[4]\n", + " params[\"geol_sres\"] = g[5]\n", + "\n", + "params[\"_comment_wtd\"] = \"Water table depth (negative = below surface)\"\n", + "params[\"wtd_m\"] = wtd_m\n", + "params[\"_comment_forcing\"] = \"Forcing\"\n", + "params[\"forcing_file\"] = forcing_filename\n", + "\n", + "write_site_config(config_path, params)\n", + "print(f\"Wrote: {config_path}\\n\")\n", + "print(open(config_path).read())\n" + ] + }, + { + "cell_type": "markdown", + "id": "snow_nb1_done", + "metadata": {}, + "source": [ + "Proceed to **Single_Column_PFCLM_netcdf.ipynb** \u2014 it will load `site_config.txt` and display the parameters for you to review before running.\n" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + }, + "language_info": { + "name": "python" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} \ No newline at end of file diff --git a/snow/pfclm_snotel_compare.ipynb b/snow/pfclm_snotel_compare.ipynb new file mode 100644 index 0000000..bee31e0 --- /dev/null +++ b/snow/pfclm_snotel_compare.ipynb @@ -0,0 +1,326 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "snow_nb3_title", + "metadata": {}, + "source": [ + "# ParFlow-CLM vs SNOTEL Comparison\n", + "\n", + "Load model output from Notebook 2, fetch SNOTEL observations, and evaluate:\n", + "- Daily SWE time series (obs vs. model)\n", + "- Scatter plot with 1:1 line\n", + "- Metrics: NSE, RMSE, |bias|%, peak SWE bias, melt-out day delta\n", + "- CW3E precipitation vs. SNOTEL precipitation (forcing quality check)\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "snow_nb3_imports", + "metadata": {}, + "outputs": [], + "source": [ + "import os\n", + "import sys\n", + "import numpy as np\n", + "import pandas as pd\n", + "import xarray as xr\n", + "import matplotlib.pyplot as plt\n", + "import matplotlib.dates as mdates\n", + "import hf_hydrodata as hf\n", + "from pathlib import Path\n", + "from datetime import datetime\n", + "\n", + "sys.path.insert(0, os.path.abspath(\"..\"))\n", + "from helpers import read_site_config\n", + "\n", + "# HydroData PIN (persists from Notebook 1)\n", + "hf.register_api_pin(\"your_email@example.com\", \"your_pin\")\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "snow_nb3_loadcfg", + "metadata": {}, + "outputs": [], + "source": [ + "# ============================================================\n", + "# Load Site Configuration\n", + "# ============================================================\n", + "static_write_dir = \"../runs/snow/Berthoud_Summit_WY2024\"\n", + "\n", + "run_dir = Path(static_write_dir)\n", + "cfg = read_site_config(run_dir / \"site_config.txt\")\n", + "\n", + "site_name = cfg[\"site_name\"]\n", + "water_year = int(cfg[\"water_year\"])\n", + "start = cfg[\"start\"]\n", + "end = cfg[\"end\"]\n", + "triplet = cfg[\"triplet\"]\n", + "region = cfg[\"region\"]\n", + "\n", + "print(f\"Site: {site_name} ({region}), WY{water_year}, triplet={triplet}\")\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "snow_nb3_loadmodel", + "metadata": {}, + "outputs": [], + "source": [ + "# ============================================================\n", + "# Load Model SWE from CLM NetCDF\n", + "# ============================================================\n", + "clm_files = sorted(run_dir.glob(\"*.out.CLM.*.nc\"))\n", + "if not clm_files:\n", + " raise FileNotFoundError(f\"No CLM output in {run_dir}. Run Notebook 2 first.\")\n", + "\n", + "clm_ds = xr.open_dataset(clm_files[0])\n", + "if \"swe_out\" not in clm_ds:\n", + " raise KeyError(f\"swe_out not in CLM output. Variables: {list(clm_ds.data_vars)}\")\n", + "\n", + "swe_h = clm_ds[\"swe_out\"].values.flatten()\n", + "dt_start = datetime.strptime(start, \"%Y-%m-%d\")\n", + "hourly_dates = pd.date_range(dt_start, periods=len(swe_h), freq=\"h\")\n", + "\n", + "model_swe_hourly = pd.Series(swe_h, index=hourly_dates, name=\"model_swe_mm\")\n", + "model_swe_daily = model_swe_hourly.resample(\"D\").mean()\n", + "\n", + "print(f\"Loaded {len(swe_h)} hourly SWE values ({len(model_swe_daily)} days)\")\n", + "print(f\" Peak model SWE: {model_swe_daily.max():.0f} mm\")\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "snow_nb3_loadobs", + "metadata": {}, + "outputs": [], + "source": [ + "# ============================================================\n", + "# Fetch SNOTEL SWE observations\n", + "# ============================================================\n", + "df = hf.get_point_data(\n", + " dataset=\"snotel\", variable=\"swe\",\n", + " temporal_resolution=\"daily\", aggregation=\"sod\",\n", + " site_ids=[triplet],\n", + " date_start=start, date_end=end,\n", + ")\n", + "\n", + "if df is None or df.empty or triplet not in df.columns:\n", + " raise RuntimeError(f\"No SNOTEL SWE data for triplet {triplet}\")\n", + "\n", + "obs_dates = pd.to_datetime(df[\"date\"])\n", + "obs_swe_daily = pd.Series(df[triplet].values.astype(float),\n", + " index=pd.DatetimeIndex(obs_dates),\n", + " name=\"obs_swe_mm\").dropna()\n", + "# SNOTEL SWE is already in mm \u2014 do NOT apply any unit conversion\n", + "print(f\"SNOTEL SWE: {len(obs_swe_daily)} daily values, peak={obs_swe_daily.max():.0f} mm\")\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "snow_nb3_metrics", + "metadata": {}, + "outputs": [], + "source": [ + "# ============================================================\n", + "# Compute metrics on common dates\n", + "# ============================================================\n", + "common = obs_swe_daily.index.intersection(model_swe_daily.index)\n", + "obs_c = obs_swe_daily.loc[common]\n", + "mod_c = model_swe_daily.loc[common]\n", + "mask = (~obs_c.isna()) & (~mod_c.isna())\n", + "obs_c = obs_c[mask]\n", + "mod_c = mod_c[mask]\n", + "\n", + "# Basic error stats\n", + "bias = float((mod_c - obs_c).mean())\n", + "rmse = float(np.sqrt(((mod_c - obs_c) ** 2).mean()))\n", + "peak_obs = float(obs_c.max())\n", + "peak_mod = float(mod_c.max())\n", + "bias_pct = (peak_mod - peak_obs) / peak_obs * 100.0 if peak_obs > 0 else np.nan\n", + "\n", + "# Nash-Sutcliffe\n", + "mu_o = obs_c.mean()\n", + "ss_res = ((mod_c - obs_c) ** 2).sum()\n", + "ss_tot = ((obs_c - mu_o) ** 2).sum()\n", + "nse = float(1 - ss_res / ss_tot) if ss_tot > 0 else np.nan\n", + "\n", + "# Peak SWE dates\n", + "peak_obs_date = obs_c.idxmax()\n", + "peak_mod_date = mod_c.idxmax()\n", + "peak_date_delta = (peak_mod_date - peak_obs_date).days\n", + "\n", + "# Melt-out date (first day SWE < 10 mm after peak)\n", + "def melt_out(s, threshold=10.0):\n", + " spring = s[s.index >= pd.Timestamp(f\"{water_year}-02-01\")]\n", + " if spring.empty:\n", + " return None\n", + " peak_idx = spring.idxmax()\n", + " post = spring[spring.index >= peak_idx]\n", + " below = post[post < threshold]\n", + " return below.index[0] if not below.empty else None\n", + "\n", + "mo_obs = melt_out(obs_c)\n", + "mo_mod = melt_out(mod_c)\n", + "mo_delta = (mo_mod - mo_obs).days if (mo_obs and mo_mod) else None\n", + "\n", + "print(f\"\\n{'='*55}\")\n", + "print(f\" {site_name} WY{water_year} \u2014 SWE Comparison\")\n", + "print(f\"{'='*55}\")\n", + "print(f\" N days: {len(obs_c)}\")\n", + "print(f\" Mean bias: {bias:+.1f} mm\")\n", + "print(f\" RMSE: {rmse:.1f} mm\")\n", + "print(f\" NSE: {nse:.3f}\")\n", + "print(f\" Peak SWE obs: {peak_obs:.0f} mm on {peak_obs_date.date()}\")\n", + "print(f\" Peak SWE model: {peak_mod:.0f} mm on {peak_mod_date.date()}\")\n", + "print(f\" Peak |bias|%: {abs(bias_pct):.1f}%\")\n", + "print(f\" Peak date delta: {peak_date_delta:+d} days\")\n", + "if mo_delta is not None:\n", + " print(f\" Obs melt-out: {mo_obs.date()}\")\n", + " print(f\" Model melt-out: {mo_mod.date()}\")\n", + " print(f\" Melt-out delta: {mo_delta:+d} days\")\n", + "else:\n", + " print(f\" Melt-out: not reached in this period\")\n", + "print(f\"{'='*55}\")\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "snow_nb3_ts", + "metadata": {}, + "outputs": [], + "source": [ + "# ============================================================\n", + "# Time Series Plot \u2014 Daily SWE\n", + "# ============================================================\n", + "fig, ax = plt.subplots(figsize=(14, 5))\n", + "ax.plot(obs_c.index, obs_c.values, color=\"black\", lw=1.2, label=\"SNOTEL Obs\")\n", + "ax.plot(mod_c.index, mod_c.values, color=\"steelblue\", lw=1.2, label=\"PF-CLM Model\")\n", + "\n", + "if mo_obs is not None:\n", + " ax.axvline(mo_obs, color=\"black\", ls=\"--\", lw=0.8, alpha=0.5)\n", + "if mo_mod is not None:\n", + " ax.axvline(mo_mod, color=\"steelblue\", ls=\"--\", lw=0.8, alpha=0.5)\n", + "\n", + "ax.set_ylabel(\"SWE [mm]\")\n", + "ax.set_title(f\"{site_name} WY{water_year} \u2014 Daily SWE | \"\n", + " f\"NSE={nse:.2f}, peak|bias|={abs(bias_pct):.0f}%, \"\n", + " f\"melt-out \u0394={mo_delta if mo_delta is not None else 'NA'}d\")\n", + "ax.legend(loc=\"upper right\")\n", + "ax.xaxis.set_major_formatter(mdates.DateFormatter(\"%b %Y\"))\n", + "ax.grid(alpha=0.3)\n", + "plt.tight_layout()\n", + "plt.show()\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "snow_nb3_scatter", + "metadata": {}, + "outputs": [], + "source": [ + "# ============================================================\n", + "# Scatter Plot \u2014 Observed vs Modeled\n", + "# ============================================================\n", + "fig, ax = plt.subplots(figsize=(7, 7))\n", + "ax.scatter(obs_c.values, mod_c.values, s=12, alpha=0.5, color=\"steelblue\", edgecolors=\"none\")\n", + "\n", + "lim_max = max(obs_c.max(), mod_c.max()) * 1.1\n", + "ax.plot([0, lim_max], [0, lim_max], \"k--\", lw=1, label=\"1:1\")\n", + "ax.set_xlim(0, lim_max)\n", + "ax.set_ylim(0, lim_max)\n", + "ax.set_xlabel(\"SNOTEL SWE [mm]\")\n", + "ax.set_ylabel(\"PF-CLM SWE [mm]\")\n", + "ax.set_title(f\"{site_name} WY{water_year} | NSE={nse:.2f}, RMSE={rmse:.0f} mm\")\n", + "ax.set_aspect(\"equal\")\n", + "ax.legend()\n", + "ax.grid(alpha=0.3)\n", + "plt.tight_layout()\n", + "plt.show()\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "snow_nb3_precip", + "metadata": {}, + "outputs": [], + "source": [ + "# ============================================================\n", + "# Forcing Quality \u2014 CW3E vs SNOTEL Precipitation\n", + "# ============================================================\n", + "# SNOTEL precip (daily totals, already mm) vs CW3E precip (cumulative).\n", + "# CW3E is typically high-biased in the mountains \u2014 check the scale.\n", + "\n", + "forcing_path = run_dir / cfg[\"forcing_file\"]\n", + "forcing_raw = np.loadtxt(forcing_path)\n", + "forc_dates = pd.date_range(start, periods=len(forcing_raw), freq=\"h\")\n", + "\n", + "# CW3E precip is mm/s -> convert to mm/hr then daily sum\n", + "cw3e_precip_hourly = pd.Series(forcing_raw[:, 2] * 3600, index=forc_dates, name=\"cw3e_precip_mm\")\n", + "cw3e_precip_daily = cw3e_precip_hourly.resample(\"D\").sum()\n", + "cw3e_precip_cum = cw3e_precip_daily.cumsum()\n", + "\n", + "# SNOTEL precip\n", + "try:\n", + " df = hf.get_point_data(\n", + " dataset=\"snotel\", variable=\"precipitation\",\n", + " temporal_resolution=\"daily\", aggregation=\"sum\",\n", + " site_ids=[triplet],\n", + " date_start=start, date_end=end,\n", + " )\n", + " if df is not None and not df.empty and triplet in df.columns:\n", + " snotel_precip_daily = pd.Series(df[triplet].values.astype(float),\n", + " index=pd.to_datetime(df[\"date\"]),\n", + " name=\"snotel_precip_mm\").dropna()\n", + " # SNOTEL precip is already in mm\n", + " snotel_precip_cum = snotel_precip_daily.cumsum()\n", + " else:\n", + " snotel_precip_cum = None\n", + "except Exception as e:\n", + " print(f\"SNOTEL precip fetch failed: {e}\")\n", + " snotel_precip_cum = None\n", + "\n", + "fig, ax = plt.subplots(figsize=(14, 5))\n", + "ax.plot(cw3e_precip_cum.index, cw3e_precip_cum.values, color=\"red\", lw=1.5, label=\"CW3E\")\n", + "if snotel_precip_cum is not None:\n", + " ax.plot(snotel_precip_cum.index, snotel_precip_cum.values, color=\"black\",\n", + " lw=1.5, label=\"SNOTEL\")\n", + " ratio = cw3e_precip_cum.iloc[-1] / snotel_precip_cum.iloc[-1] if snotel_precip_cum.iloc[-1] > 0 else np.nan\n", + " ax.set_title(f\"{site_name} WY{water_year} \u2014 Cumulative Precip | \"\n", + " f\"CW3E/SNOTEL = {ratio:.2f}\")\n", + "else:\n", + " ax.set_title(f\"{site_name} WY{water_year} \u2014 CW3E Cumulative Precip (no SNOTEL precip)\")\n", + "\n", + "ax.set_ylabel(\"Cumulative precip [mm]\")\n", + "ax.xaxis.set_major_formatter(mdates.DateFormatter(\"%b %Y\"))\n", + "ax.legend(loc=\"upper left\")\n", + "ax.grid(alpha=0.3)\n", + "plt.tight_layout()\n", + "plt.show()\n", + "\n", + "clm_ds.close()\n" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + }, + "language_info": { + "name": "python" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} \ No newline at end of file