diff --git a/.claude/skills/audit-pr/SKILL.md b/.claude/skills/audit-pr/SKILL.md index c56a73006..4fe66aad4 100644 --- a/.claude/skills/audit-pr/SKILL.md +++ b/.claude/skills/audit-pr/SKILL.md @@ -18,7 +18,7 @@ gh pr diff {pr} --name-only - **Step 1 - Context**: Gather files, classify changes - **Step 2 - Code Style**: lint-code -- **Step 3 - Scientific**: Physics validation (if applicable) +- **Step 3 - Scientific**: Physics validation + description rigour (if applicable) - **Step 4 - Testing**: Coverage, FIRST principles - **Step 5 - Docs**: CHANGELOG, PR description - **Step 6 - Governance**: Check feature status tags (see below) @@ -27,6 +27,41 @@ gh pr diff {pr} --name-only - **Step 9 - Approval**: Wait for human confirmation - **Step 10 - Post**: Only after approval +## Scientific PR Description Rigour (Step 3) + +For PRs that introduce or materially change a scientific component +(physics, numerics, calibration, benchmark), the PR description + thread +must expose the science clearly so reviewers can judge it in-place. This +is a description-quality check, not a separate artefact requirement — +evidence lives in the PR, not in a sidecar dashboard. + +**A PR is "scientific" if it touches any of**: + +- `src/suews/src/suews_phys_*.f95` +- `src/suews/src/suews_ctrl_daily_state.f95` +- `src/supy/data_model/` (new physics-facing fields or enums) +- New YAML parameters that affect model output +- A CHANGELOG entry tagged `[feature]` or `[change]` that is not pure + tooling / build / docs + +**Required in the PR description (or posted as comments in the thread)**: + +- [ ] **Methodology** — what was changed and how it was developed/validated + (governing equation, numerical approach, calibration procedure) +- [ ] **Scientific decisions** — key choices made and rationale + (why this parameterisation, why this threshold, what alternatives + were considered) +- [ ] **Results** — quantitative evidence that the change works, posted + in the PR (numbers, plots, before/after comparisons, test outputs). + A bare "tests pass" is not sufficient for scientific changes. + +**If missing or thin**: flag as a blocking review comment and ask the +author to expand the PR body or post a follow-up comment with the +missing pieces. Do not approve a scientific PR whose description reads +as purely mechanical ("refactored X", "added Y") without scientific +narrative. If scope is ambiguous, ask the author to classify rather +than guess. + ## Governance Check (Step 6) For PRs that add new features or breaking changes: @@ -73,6 +108,7 @@ Details: `references/style-checks.md` |----------|--------| | Code Style | PASS/FAIL | | Scientific | PASS/FAIL/N/A | +| Description rigour | PASS/FAIL/N/A | | Testing | PASS/FAIL | | Docs | PASS/FAIL | diff --git a/.claude/skills/audit-pr/references/review-checklist.md b/.claude/skills/audit-pr/references/review-checklist.md index 0545ee1d9..e000e2328 100644 --- a/.claude/skills/audit-pr/references/review-checklist.md +++ b/.claude/skills/audit-pr/references/review-checklist.md @@ -66,6 +66,56 @@ Detailed checklist for comprehensive PR review. --- +## Scientific PR Description Rigour Checklist + +Applies to any PR that introduces or materially changes a scientific +component. The evidence lives in the PR itself — description and +thread comments — not in a sidecar artefact. Missing pieces are a +**blocking** finding: ask the author to expand the body or post +follow-up comments before approving. + +### Triggers — a PR is "scientific" if it touches + +- [ ] Any file matching `src/suews/src/suews_phys_*.f95` +- [ ] `src/suews/src/suews_ctrl_daily_state.f95` +- [ ] `src/supy/data_model/` (new physics-facing fields, enums, validators) +- [ ] New YAML parameters that alter model output +- [ ] A CHANGELOG entry tagged `[feature]` or `[change]` that is not pure + tooling / build / docs +- [ ] Benchmark additions or regressions against published reference runs + +### Required in the PR description / thread + +- [ ] **Methodology** — what was changed and how it was developed or + validated (governing equation, numerical approach, calibration + procedure, data source) +- [ ] **Scientific decisions** — the key choices made and why: which + parameterisation, which threshold, which alternatives were + considered and rejected +- [ ] **Results** — quantitative evidence that the change works, posted + in the PR as numbers, figures, or comparison tables. Plots can be + embedded via `draft-post`'s figure upload or linked from a + comment. A bare "tests pass" is not sufficient for a scientific + change +- [ ] **Scope** — which issue this resolves, which follow-ups remain + +### How to handle a thin or mechanical PR description + +When a scientific PR arrives with a description that only names files +changed, or says "refactored X / added Y" without scientific narrative: + +1. Flag it as "Description rigour: FAIL" in the review summary +2. Post a blocking review comment listing which of the four items above + are missing, and invite the author to either edit the PR body or + post the missing content as a follow-up comment in the thread +3. Do not approve the PR for merge until the description covers + methodology, decisions, and quantitative results +4. If scope is ambiguous (is this really scientific, or is it tooling?), + ask the author to classify in the PR description and defer to a + maintainer rather than guess + +--- + ## Testing Checklist ### Coverage Requirements @@ -176,6 +226,7 @@ PR can be merged when ALL of: - [ ] CI tests pass - [ ] Code review approved - [ ] Scientific review approved (if physics changes) +- [ ] PR description covers methodology, scientific decisions, and results (if scientific change) - [ ] Documentation updated - [ ] CHANGELOG entry present - [ ] No blocking issues diff --git a/CHANGELOG.md b/CHANGELOG.md index 97747ca4d..50c263e74 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -65,6 +65,7 @@ EXAMPLES: ### 18 Apr 2026 +- [feature][experimental] Added parameter-sweep calibration tool for moisture-aware LAI (#1292). `scripts/verify/moisture_phenology_sweep.py` scans any of the six moisture parameters (`w_wilt`, `w_opt`, `f_shape`, `w_on`, `w_off`, `tau_w`) over a user-supplied or default value list, reruns SUEWS per value, and writes a dose-response summary (mean LAI, RMSE vs `LAIType = 0` baseline, seasonal amplitude, green-up DOY) plus a dual-axis sensitivity plot to `.context/gh1292//sweep_.{json,png}`. The `--dry-start` flag depletes vegetation `soilstore_surf` so the moisture gate engages at otherwise well-watered samples. Includes a pairwise-validator co-adjustment that keeps `w_opt > w_wilt` and `w_on > w_off` satisfied during extreme sweep values. Design note extended with a new calibration-methodology appendix (London sensitivity results, FLUXNET calibration roadmap, status of the optional SDD term). - [maintenance] Register `physics` / `api` pytest markers as an orthogonal axis to the existing tier markers (#1300) - `physics` covers numerical / binary correctness (one canonical Python per `(OS, arch)` suffices post-abi3); `api` covers Python wrapper surface (pandas / numpy / pydantic, CLI, `SUEWSSimulation`) and needs full `(platform x Python)` coverage - Applied module-level `pytestmark` to all 43 non-UMEP test files; the 5 UMEP test files pick up `api` via the existing `pytest_collection_modifyitems` hook @@ -79,6 +80,8 @@ EXAMPLES: ### 17 Apr 2026 +- [feature][experimental] Activated moisture-aware LAI phenology numerics (`LAIType = 2`) with the Design C formulation (#1292). Replaces the PR1 no-op branch with (a) a Jarvis-style water-stress factor on `delta_GDD` parametrised by `w_wilt`, `w_opt`, and `f_shape` on dimensionless relative soil water `w = 1 - smd/smdcap`, and (b) a CLM5-style persistence latch (Lawrence et al. 2019) that freezes thermal accumulation when the `tau_w`-day running mean of `w` drops below `w_off` and reopens when it rises above `w_on` while a thermal companion (Tbar >= BaseT_GDD) is met. `PHENOLOGY_STATE.leaf_on_permitted` now defaults to true so well-watered sites do not spend `tau_w` days ramping up, and `wbar_id` is seeded from the first measured `w` on the first call. Default moisture thresholds (`w_wilt = 0.15`, `w_opt = 0.40`, `w_on = 0.35`, `w_off = 0.20`) degrade gracefully to thermal-only behaviour at well-watered sites -- the bundled London sample is bit-identical to `LAIType = 0`. Aggressive thresholds or depleted soil-store initial conditions produce the expected strictly lower LAI; a dry-start regression test locks in the spring green-up delay. Calibration against FLUXNET2015 dryland / monsoon sites is deferred to a follow-up. Internal design note: `dev-ref/design-notes/gh1292-moisture-phenology.md`. +- [feature][experimental] Added scaffolding for moisture-aware LAI phenology (`LAIType = 2`) behind a new per-vegetation-surface switch (#1292). Extends `LAI_PRM` with six parameter slots (`w_wilt`, `w_opt`, `f_shape`, `w_on`, `w_off`, `tau_w`) and `PHENOLOGY_STATE` with three per-veg-surface state slots (`wbar_id`, `w_id_prev`, `leaf_on_permitted`); bumps the Rust/C bridge schema for `LAI_PRM`, `PHENOLOGY_STATE`, `LC_DECTR_PRM`, `LC_EVETR_PRM`, `LC_GRASS_PRM` to version 2. The `LAIType = 2` branch is a no-op that reproduces `LAIType = 0` bit-identically in this release -- Design C numerics (Jarvis water-stress factor plus CLM5-style persistence trigger) land in a later PR. Internal design note: `dev-ref/design-notes/gh1292-moisture-phenology.md`. - [maintenance] Remove NumPy build-time dependency (post-f2py cleanup) (#1298) - [maintenance] Enforce numpy-style docstrings via ruff `D` rules (#1294) - Removed D100-D105/D107 from the global gradual-adoption ignore list in `.ruff.toml` diff --git a/dev-ref/design-notes/README.md b/dev-ref/design-notes/README.md new file mode 100644 index 000000000..f160dc860 --- /dev/null +++ b/dev-ref/design-notes/README.md @@ -0,0 +1,13 @@ +# Design Notes + +Internal design notes for substantive SUEWS design decisions that need a durable record before implementation. + +## Convention + +- One markdown file per issue-scoped design decision. +- Filename: `gh-.md` (e.g. `gh1292-moisture-phenology.md`). +- Audience: SUEWS core developers. These are internal working documents, not user-facing docs. +- Link upstream issues, papers, and community threads; do not paraphrase public threads verbatim. +- British English, plain ASCII where practical, consistent with `.claude/rules/00-project-essentials.md`. + +Design notes that turn into shipped features should be referenced from the PR body and may later be retired or moved under `docs/` if they become user-relevant. diff --git a/dev-ref/design-notes/gh1292-moisture-phenology.md b/dev-ref/design-notes/gh1292-moisture-phenology.md new file mode 100644 index 000000000..399e62025 --- /dev/null +++ b/dev-ref/design-notes/gh1292-moisture-phenology.md @@ -0,0 +1,365 @@ +# Moisture-aware phenology for rainfall-driven sites + +**Upstream issue:** [GH-1292](https://github.com/UMEP-dev/SUEWS/issues/1292) +**Related:** [GH-1291](https://github.com/UMEP-dev/SUEWS/issues/1291) (observed-LAI forcing silently ignored) +**Status:** Design, pre-implementation. Needs scientific review before any code lands. +**Labels:** `1-feature`, `2-module:dailystate`, `4-needs-science`. +**Scope note:** This is an internal design record for core developers. The public issue on GitHub is sufficient for community visibility -- we do not paraphrase community forum threads here. + +--- + +## 1. Problem statement + +SUEWS phenology currently evolves LAI through growing-degree-day (GDD) and senescence-degree-day (SDD) thresholds computed from daily mean air temperature. There is no precipitation or soil-moisture term anywhere in the governing equations. + +Thermal-only accumulation is implemented in `src/suews/src/suews_phys_dailystate.f95`: + +- `update_GDDLAI` signature at lines 538-560. +- GDD/SDD accumulation loop at lines 582-617. +- LAI power-law update at lines 627-672. + +For sites where the annual LAI cycle is controlled by rainfall rather than temperature (monsoon grasslands, arid/semi-arid systems, pulsed wet-season savannas, double-cropping agricultural systems) the scheme cannot reproduce observed seasonality at any tuning of `BaseT`, `GDDFull`, `SDDFull`, or `LAIPower`. The limitation is acknowledged in Omidvar et al. (2022), *Geosci. Model Dev.* 15, 3041-3078, Appendix D (Figs D1/D2; US-MMS deciduous broadleaf works, US-SRG monsoon grassland fails). The paper explicitly flags hydrological feedback as future work. + +--- + +## 2. Empirical evidence: the thermal scheme's failure mode + +Ting's prior FLUXNET2015 evaluation (the `SUEWS-FLUXNET2015` repo that produced Omidvar et al. 2022) already derived per-site thermal-only LAI parameters for ~100 sites via the `LAIModel` Mathematica package. The pattern at moisture-limited sites is diagnostic: the fit collapses to near-zero seasonality rather than signalling failure. + +Selected rows from the derived parameters (`data/csv-prm/prm-LAI-Albedo.csv` in that repo): + +- **AU-ASM** (Alice Springs, mulga, semi-arid, MAP ~300 mm): LAImin = 0.20, LAImax = 0.30, GDDFull = 93. The fitted seasonal amplitude is essentially zero and GDDFull is tiny -- the thermal driver saturates year-round, so LAI is effectively constant. +- **AU-DaS** (Daly River savanna, tropical monsoon): LAImin = 0.97, LAImax = 1.86, GDDFull = 351. Amplitude under-estimated compared with MODIS. +- **AU-Gin** (Gingin, Mediterranean woody savanna): LAImin = 0.60, LAImax = 1.09, GDDFull = 239. Again compressed. + +Contrast with thermally-controlled sites fitted by the same pipeline: + +- **US-MMS** (Morgan-Monroe, temperate deciduous broadleaf): LAImin = 1.01, LAImax = 5.00, GDDFull = 69. Strong, physically sensible cycle. +- **US-Syv** (Sylvania, temperate mixed forest): LAImin = 0.66, LAImax = 4.81. + +Taken together, the thermal-only fits at dryland and monsoon sites do not merely mistime leaf-on -- they erase the seasonality altogether. This is the empirical evidence that motivated Omidvar Appendix D and is the basis for introducing a moisture control. + +**Data available locally.** The `SUEWS-FLUXNET2015` repo under `/Users/tingsun/Dropbox (Personal)/6.Repos/SUEWS-FLUXNET2015/` contains per-site pre-processed forcings and SUEWS simulation outputs in HDF5, MODIS LAI, SoilGrids soil properties, and the derived-parameters CSVs. The Mathematica parameter-fitting package (`LAIModel.m` from `/Users/tingsun/Dropbox (Personal)/6.Repos/LAI-model/`) is the pipeline to extend for a moisture-aware refit. The `AmeriFluxEvaluation` repo has companion Python notebooks (`GDD_SDD_LAI_OBS.ipynb`, `LAI.ipynb`) that explore GDD/SDD vs observed LAI. + +--- + +## 3. Verified facts about the current code + +Confirmed in the repo at time of writing (branch `sunt05/gh1292-lai-moisture`): + +- **Driver call order.** In `src/suews/src/suews_ctrl_driver.f95`: + - Step A1 `SUEWS_cal_DailyState` at line 320. + - Step A10 `SUEWS_cal_SoilState` at line 3085. + - Consequence: when end-of-day phenology fires inside DailyState, current-timestep SMD is not yet computed. DailyState therefore sees the previous day's end-of-day SMD, which is exactly what a moisture-aware phenology scheme wants. No driver-level reordering is needed. +- **`LAIType` switch.** Declared as `FlexibleRefValue(int)` in `src/supy/data_model/core/site.py` around line 401. Existing values: `0` = original, `1` = new high-latitude variant. Adding values 2+ is the natural extension. +- **`LAIParams`.** Defined in `src/supy/data_model/core/site.py` lines 348-496, nested inside `EVETRParams` / `DECTRParams` / `GRASSParams`. `to_df_state` / `from_df_state` already tolerate missing columns with defaults (pattern at 441-450). +- **Per-surface vegetation SMD.** Computed by `cal_smd_veg()` in `src/suews/src/suews_phys_waterdist.f95` and carried as `modState%hydroState%smd_nsurf`. Capacity information accessible via the soil-store parameters on the same state object. +- **`LAI_obs` forcing path.** `LAI_obs` is read into the `forcing` struct in the dailystate module and used only when `LAICalcYes == 0`. That branch is currently broken -- tracked in GH-1291. The interim workaround recommended to users (prescribe observed LAI) therefore does not work today. + +--- + +## 4. Literature review: how other models couple water to phenology + +Our design decision should not invent a new formulation if a published scheme already fits. The literature divides into three schools. + +### 4.1 Empirical bioclimatic indices + +- **Jolly, Nemani & Running (2005).** *Glob. Change Biol.* 11, 619-632. doi:10.1111/j.1365-2486.2005.00930.x. The **Growing Season Index (GSI)** is a multiplicative product of piecewise-linear indicators on minimum temperature, vapour-pressure deficit and photoperiod. No soil-moisture term; water stress enters only indirectly via VPD. This is close to what SUEWS does now, minus the photoperiod term, and is known to be inadequate for monsoon systems because leaf-on there follows rainfall pulses regardless of VPD. +- **Stockli et al. (2008, 2011).** *JGR Biogeosci.* doi:10.1029/2008JG000781 and doi:10.1029/2010JG001545. Prognostic extension of GSI with a fourth multiplicative factor on soil moisture; MODIS FPAR/LAI assimilated for parameter estimation. The SM factor is Jarvis-style. +- **Caldararu, Purves & Palmer (2014).** *Biogeosciences* 11, 763-778. doi:10.5194/bg-11-763-2014. Prognostic LAI optimising net canopy carbon return under light, temperature and soil-moisture constraints. Grasslands emerge explicitly as water-limited in the Bayesian fit. +- **Knorr et al. (2010).** *JGR Biogeosci.* doi:10.1029/2009JG001119. Generic 3-PFT-class BETHY phenology where leaf dynamics respond prognostically to temperature, water and light. Documents a known oscillation pathology under sustained drought -- a design warning for any NPP-coupled leaf balance. +- **Forkel et al. (2014).** *Glob. Change Biol.* 21, 3414-3435. doi:10.1111/gcb.12950. LPJmL with explicit water-availability factor; finds water is a co-dominant control on start-of-season and peak greenness globally, dominant in drylands. + +### 4.2 Major LSM precedents + +- **CLM5 (Lawrence et al. 2019).** *J. Adv. Model. Earth Syst.* 11, 4245-4287. doi:10.1029/2018MS001583. Defines three phenology types: evergreen, seasonally deciduous (GDD-triggered), and **stress-deciduous**. Stress-deciduous triggers leaf-on when both a root-zone soil-water threshold **and** a soil-temperature threshold persist for ~15 days; leaf-off triggered when either falls below the off-threshold. This is the closest published precedent for what SUEWS needs. +- **Noah-MP (Niu et al. 2011).** *JGR Atmos.* 116, D12109. doi:10.1029/2010JD015139. Dynamic-vegetation option carries prognostic leaf/stem/root/wood C pools; LAI = f(leaf C, SLA). SM enters indirectly through GPP stomatal stress. Known to underperform in arid/semi-arid unless retuned. +- **JULES/TRIFFID (Best et al. 2011, Clark et al. 2011).** *GMD* 4, 677-699 and 701-722. GDD-controlled leaf turnover; no direct SM trigger. Water acts on GPP through the beta factor. +- **ORCHIDEE (Krinner et al. 2005).** *Glob. Biogeochem. Cycles* 19, GB1015. doi:10.1029/2003GB002199. PFT-specific schemes: GDD for temperate trees, **GDD plus moisture trigger for tropical raingreen**, and a **moisture-onset trigger** for C4 grass; water stress accelerates senescence. A direct precedent for per-surface differentiation of the moisture response. +- **LPJ-DGVM / LPJmL (Sitch et al. 2003; Smith et al. 2014).** *Glob. Change Biol.* 9, 161-185; *Biogeosciences* 11, 2027-2054. The **raingreen PFT** holds LAI at maximum while water-stress scalar W >= 0.35 and drops leaves when W < 0.35 -- the simplest published threshold form. +- **ISBA-A-gs / SURFEX (Gibelin et al. 2006).** *JGR Atmos.* 111, D18102. doi:10.1029/2005JD006691. Implicit phenology: LAI = leaf biomass / constant ratio; biomass driven by A-gs photosynthesis so SM enters via stomatal conductance. Reproduces precipitation-controlled interannual LAI variability reasonably without a dedicated SM trigger. + +### 4.3 Cropping / rangeland schemes with explicit water stress + +- **WOFOST / DSSAT / APSIM.** Jarvis-style water-stress factor Ta/Tp (actual over potential transpiration, or TURFAC/SWFAC in DSSAT) multiplies daily leaf area expansion and accelerates senescence. Jones et al. (2003) *Eur. J. Agron.* 18, 235-265. Holzworth et al. (2014) *Environ. Model. Softw.* 62, 327-350. +- **Biome-BGC (White et al. 2000; Thornton et al. 2002).** Daily bucket soil water; GPP-driven leaf allocation. Phenology timing uses climatologic triggers, LAI accumulation scales with GPP. +- **Grassland-specific.** LINGRA (Schapendonk et al. 1998), PaSim (Riedo et al. 1998): leaf expansion * water-stress factor. DayCent/CENTURY (Parton et al. 1998) was historically criticised for thermal-only phenology and has been patched; Liu et al. (2022), *Agric. For. Meteorol.* 327, 109230, doi:10.1016/j.agrformet.2022.109230, propose a nonlinear soil-water senescence function specifically for grassland brown-down. + +### 4.4 Dryland and monsoon studies + +- **Dahlin, Fisher & Lawrence (2015).** *Biogeosciences* 12, 5061-5074. doi:10.5194/bg-12-5061-2015. Tests CLM's drought-deciduous trigger against FLUXNET and MODIS; the double-threshold design works but is sensitive to rooting depth and the 15-day persistence window. +- **Broich et al. (2014).** *Biogeosciences* 11, 5181-5198. doi:10.5194/bg-11-5181-2014. Australian semi-arid LSP; precipitation alone explains ~80 % of interannual variance in greening-season length on the North Australian Tropical Transect. +- **Zhang et al. (2005).** *Int. J. Remote Sens.* 26, 4517-4533. doi:10.1080/01431160500113971. 1-3 month precipitation-to-NDVI lags across Africa. +- **Jolly & Running (2004).** *Glob. Change Biol.* 10, 303-308. doi:10.1111/j.1365-2486.2003.00701.x. Kalahari field evidence that **soil-water potential** outperforms cumulative rainfall as a leaf-out predictor. +- **Hickler et al. (2005).** *Geophys. Res. Lett.* 32, L21415. doi:10.1029/2005GL024370. Precipitation controls on Sahel greening trend. + +### 4.5 Common mathematical forms + +- **Multiplicative Jarvis-style stress factor** f(w) in [0,1]: used by WOFOST, DSSAT, APSIM, Biome-BGC, ORCHIDEE (beta), Stockli 2008, ISBA implicitly. Two parameters (wilting point / field capacity or equivalent thresholds), optional curvature exponent. +- **Threshold / persistence trigger.** Leaf-on when running-mean SM > W_on for tau days; leaf-off when SM < W_off. Used by CLM5 stress-deciduous, LPJ raingreen (single 0.35 threshold). Three parameters. +- **Logistic** f = 1 / (1 + exp(-k(w - w50))). Smooth alternative; Caldararu et al., Liu et al. 2022. Two parameters. +- **Additive senescence stress** dSDD_eff = dSDD + gamma*(1 - f(w)). Used by ORCHIDEE, WOFOST family. One extra parameter. +- **Bucket trigger** growth allowed only when theta/theta_FC > theta_crit. Stockli 2008 fourth factor; raingreen variants. One parameter. + +### 4.6 Observational consensus + +Across the remote-sensing phenology literature (Zhang et al. 2005; Broich et al. 2014, 2015; Wang, Rich & Price 2003 *IJRS* 24, 2345-2364; Jolly & Running 2004) the consistent findings are: + +- Antecedent precipitation over a 1-3-month window is a competitive predictor of onset timing in drylands and monsoon systems. +- Root-zone soil moisture is a better predictor than rainfall itself once a simple bucket model is available. +- In monsoon grass/savanna leaf-on is triggered by rainfall onset regardless of GDD state, so a multiplicative stress factor alone is necessary but not sufficient -- a trigger is needed. + +--- + +## 5. Design space for SUEWS + +Given the literature, four candidate designs remain on the table: + +- **Design A: Multiplicative SM factor only.** Keep the GDD accumulator; multiply dGDD by f(w). Two new parameters. Closest to WOFOST-family schemes. Will **reduce** GDD accumulation during dry periods but cannot **prevent** leaf-on if GDD has already saturated -- a real risk in dryland sites where the thermal driver saturates year-round (cf AU-ASM in section 2 above). +- **Design B: Persistence trigger only.** Gate leaf-on on a 15-day running mean of SM exceeding W_on. Three parameters (W_on, W_off, tau). Equivalent to CLM5 stress-deciduous. Solves the AU-ASM pathology but discards the temperature signal when it matters (e.g. high-latitude continental where cold snaps should still delay leaf-out). +- **Design C: Hybrid (multiplicative + trigger).** Apply both: f(w) multiplies dGDD during the build-up phase, and leaf-on is additionally gated on the persistence trigger. Four parameters but no redundancy -- each addresses a distinct failure mode. This is effectively CLM5's stress-deciduous logic wrapped around the SUEWS GDD accumulator. +- **Design D: GPP-driven prognostic LAI.** Replace the GDD/SDD scheme entirely with a Noah-MP or ISBA-A-gs style leaf-biomass pool that accumulates from net GPP. Cleanest scientifically but a full architectural refactor; requires a leaf-C state and a Q10 or similar respiration term. Out of scope for an incremental PR roadmap. + +Sub-choices orthogonal to A-D: + +- **Water state variable**: cumulative antecedent precipitation (API) vs soil-moisture deficit (SMD) vs volumetric soil water. SUEWS already carries per-surface SMD via `cal_smd_veg`; adopting API would introduce a new diagnostic with no clear advantage (Jolly & Running 2004: SM potential outperforms cumulative rainfall). +- **Applied to growth only, or growth and senescence**: symmetric gating (ORCHIDEE-style) captures drought-driven brown-down. Option to layer the senescence term in a later PR without changing PR1/PR2 plumbing. +- **Per-surface versus shared parameters**: EVETR, DECTR, GRASS differ in rooting depth and drought strategy; per-surface parameters are the ORCHIDEE precedent. +- **Backward compatibility**: must be opt-in via a new `LAIType` value; existing YAMLs produce bit-identical output. + +--- + +## 6. Recommendation + +**Design C (Hybrid)**, introduced as `LAIType = 2`. Rationale: + +- **Addresses both empirical failure modes.** AU-ASM thermal saturation (needs trigger) and monsoon-grass rainfall pulse (needs multiplicative response during build-up) are both captured. A single-mechanism scheme would fix one and expose the other. +- **Directly portable from CLM5 stress-deciduous.** Lawrence et al. 2019 provides the governing equations and default thresholds; we adopt rather than invent. +- **Reuses SUEWS hydrology.** `modState%hydroState%smd_nsurf` from `cal_smd_veg` is the natural input. Previous-day SMD is already available thanks to driver ordering (section 3). +- **Backward-compatible.** `LAIType` default stays `0`; all existing YAMLs and existing fitted parameters in `SUEWS-FLUXNET2015/data/csv-prm/prm-LAI-Albedo.csv` are untouched. +- **Incremental.** Senescence gating (Design C plus option) can land as a later PR without altering the PR1/PR2 scaffolding. + +### 6.1 Governing equations (LAIType = 2) + +Operating on previous-day vegetation SMD per surface iv: + +``` +w(iv) = clamp(1 - xsmd_prev(iv) / smdcap(iv), 0, 1) # relative root-zone water +wbar(iv) = tau-day running mean of w(iv) # persistence window + +# Multiplicative stress on GDD build-up (Jarvis) +f_SMD(iv) = clamp((w(iv) - w_wilt(iv)) / (w_opt(iv) - w_wilt(iv)), + 0, 1) ** shape(iv) +dGDD_eff = dGDD * f_SMD(iv) + +# CLM5-style persistence trigger +if not leaf_on_permitted(iv): + leaf_on_permitted(iv) = (wbar(iv) >= w_on(iv)) AND (Tbar >= Tbar_on) +if leaf_on_permitted(iv) AND (wbar(iv) < w_off(iv)): + leaf_off_trigger(iv) = TRUE # force SDD accumulation +``` + +LAI update then reuses the existing power-law on the gated dGDD_eff and on the (possibly triggered) SDD accumulator. + +### 6.2 Parameter inventory (per vegetated surface) + +- `w_wilt`: dimensionless, lower bound of Jarvis factor. Default 0.15 (wilting-point equivalent). Range [0.0, 0.5]. +- `w_opt`: dimensionless, upper bound. Default 0.4. Range [0.1, 0.9]. Validator: `w_opt > w_wilt`. +- `shape`: dimensionless exponent, default 1.0. Range [0.25, 4.0]. +- `w_on`: dimensionless trigger threshold for leaf-on. Default 0.35 (Sitch et al. 2003 raingreen value). Range [0.05, 0.95]. +- `w_off`: trigger threshold for leaf-off. Default 0.20. Range [0.05, 0.95]. Validator: `w_off < w_on`. +- `tau`: persistence window in days. Default 15 (CLM5). Range [1, 60]. + +Total: six per-surface parameters. Four can be held at sensible defaults; two (`w_wilt`, `w_opt`) may need site calibration against MODIS LAI -- this is what the existing FLUXNET2015 fitting pipeline will estimate for us. + +### 6.3 Defaults when not calibrated + +Defaults are defensible without site calibration: `w_on=0.35` (Sitch 2003), `w_off=0.20`, `tau=15` (CLM5), `w_wilt=0.15`, `w_opt=0.4`, `shape=1.0`. Under these defaults the scheme degrades gracefully: in well-watered regimes `f_SMD -> 1` and the trigger permits leaf-on so behaviour collapses toward the thermal scheme; under sustained drought the trigger prevents spurious leaf-on from GDD saturation. + +--- + +## 7. Implementation sketch (no code in this branch) + +### 7.1 Fortran + +- Extend `update_GDDLAI` signature in `src/suews/src/suews_phys_dailystate.f95` (lines 538-560) with: + - `smd_id_prev(nvegsurf)`, `smdcap_id_prev(nvegsurf)` -- previous-day SMD and capacity, mm. + - `tbar_id_prev(nvegsurf)` already effectively available via `Tmin_id_prev` / `Tmax_id_prev` -- no new arg. + - Persistent state for the `tau`-day running-mean `wbar` and the boolean `leaf_on_permitted`: add to the daily-state derived type in `src/suews/src/suews_ctrl_type.f95` parallel to `Tmin_id`, `Tmax_id`. +- In the accumulation loop (lines 582-617) branch on `LAItype(iv) == 2`: compute `f_SMD`, update `wbar`, evaluate the trigger, gate `dGDD` and `dSDD`. +- Caller `SUEWS_cal_DailyState` in the same file (~lines 100-300): thread `modState%hydroState%smd_nsurf` and soil-store capacity through to `update_GDDLAI`. + +### 7.2 Python data model + +- Extend `LAIParams` in `src/supy/data_model/core/site.py` (lines 348-496) with six new optional `FlexibleRefValue(float)` fields listed in section 6.2, each with default, range, unit and description in line with the existing `LAIParams` style. +- Update the `laitype` description at line 403: `0: original, 1: new high latitude, 2: moisture-aware (hybrid multiplicative + persistence trigger, CLM5-style)`. +- Extend `to_df_state` / `from_df_state` to (de)serialise the six new fields with defaults; legacy YAMLs must round-trip unchanged. +- Pydantic validators: `w_opt > w_wilt`, `w_off < w_on`, all fields within stated ranges; when `laitype == 2` and fields unset, emit an info-level log that defaults apply (do not raise). + +### 7.3 Backward-compatibility contract + +`LAIType` default stays `0`. Output must be bit-identical for existing YAMLs after the eventual implementation. A regression test built on `SUEWSSimulation.from_sample_data()` covering `LAIType in {0, 1}` is a required part of the first scaffolding PR, independent of the numerical changes. + +--- + +## 8. Validation plan + +The existing `SUEWS-FLUXNET2015` pipeline is the validation harness. No new data collection needed. + +### 8.1 Site set + +- **Thermal control sites** (scheme must not regress): US-MMS (DBF), US-Syv, DE-Hai, DK-Sor. +- **Semi-arid / dryland** (scheme should recover seasonality): AU-ASM (mulga), AU-Stp (Sturt Plains), AU-Dry, AU-Cpr, US-SRG (Santa Rita grassland -- needs to be added to the archive if not already), US-Wkg (Walnut Gulch), US-Ton (Mediterranean savanna). +- **Monsoon / tropical** (scheme should retime leaf-on to rainfall pulse): AU-DaS, AU-DaP (Daly River monsoon), AU-Ade (Adelaide River). +- **Mediterranean** (strong summer drought, winter growth): AU-Gin, FR-Pue, IT-Noe. + +### 8.2 Fitting pipeline + +Extend the Mathematica `LAIModel.m` package (at `/Users/tingsun/Dropbox (Personal)/6.Repos/LAI-model/LAIModel.m`, functions `fitTSLAI`, `calLAIGDDSDD`, `simLAI`) to expose the hybrid Design C model as a fittable form. Fit the six per-surface parameters per site against MODIS LAI; retain `BaseTGDD`, `GDDFull` etc. from the existing thermal fits unchanged. + +### 8.3 Success metrics + +- **Regression (thermal sites)**: new fit with `LAIType = 2` recovers the existing thermal fit within a small tolerance (i.e. moisture modifier is near-unity at these sites). Acceptance: RMSE of simulated vs MODIS LAI not worse than the thermal-only baseline. +- **Dryland sites**: LAImin-LAImax seasonal amplitude recovers to within some fraction of the MODIS-observed amplitude -- where thermal-only collapsed to near-flat (see section 2). +- **Monsoon sites**: green-up timing lag (modelled vs MODIS) reduced relative to the thermal-only baseline. +- **Energy-balance consistency**: evaluate albedo, Qh and Qe against FLUXNET eddy-covariance data to ensure the improved LAI does not degrade the surface energy balance. + +### 8.4 Reporting + +A comparison mirroring Omidvar et al. 2022 Appendix D: reproduce the thermal-only failure at US-SRG (and add AU-ASM and AU-Stp as additional dryland stress tests), then show the Design C fix. Package the result as a short note for the next SUEWS release and as supplementary material if the scheme reaches publication. + +--- + +## 9. Incremental PR roadmap + +Follow `dev-ref/FEATURE_DEVELOPMENT_WORKFLOW.md` (bottom-up): + +- **PR1 -- scaffolding (no numerical change).** + - `LAIParams` extended with the six new optional fields; `laitype` description updated. + - `update_GDDLAI` signature extended; caller plumbs SMD, capacity and running-mean state through. + - New daily-state slots (`wbar`, `leaf_on_permitted`) added to the derived type in `suews_ctrl_type.f95` but consumed only under `LAItype == 2`. + - `LAItype == 2` branch is a no-op: `f_SMD = 1`, trigger always permits leaf-on. Thermal outputs bit-identical to `LAItype = 0`. + - Regression test: `SUEWSSimulation.from_sample_data()` on `master` vs branch for `LAItype in {0, 1, 2}`. +- **PR2 -- implementation of Design C.** + - Enable the Jarvis-style `f_SMD` and the CLM5-style trigger. + - Synthetic dry-then-wet unit test: identical temperature, zero precipitation for 20 days then normal for 20 days; assert effective GDD reduction and delayed leaf-on. + - Soft validator warnings when fields unset. +- **PR3 -- calibration and optional senescence gating.** + - Extend `LAIModel.m` to fit Design C; refit FLUXNET2015 dryland and monsoon sites. + - Publish per-IGBP default `w_wilt`, `w_opt`, `w_on`, `w_off` from the fitting results. + - Optional: add additive senescence stress term (ORCHIDEE / WOFOST style) behind the same `LAIType = 2` switch. + +Prerequisite: **GH-1291** (observed-LAI forcing silently ignored) should be resolved before PR3 so users have a viable fallback while calibrated parameters propagate into the sample configs. + +--- + +## 10. Open scientific questions + +- **Per-surface vs shared parameters.** Rooting depth and drought strategy differ between EVETR, DECTR and GRASS; per-surface is the ORCHIDEE precedent but increases parameter count. Can we collapse to IGBP-class defaults with per-site override? +- **Symmetric gating of SDD.** Drought-driven brown-down (ORCHIDEE, Liu et al. 2022) is well-supported observationally. PR2 vs PR3 split is a pragmatic choice; does the monsoon-site fit collapse without it? +- **Irrigated urban grass.** Omidvar et al. 2022 flagged irrigation as a confounder. Should the scheme bypass the moisture response when a site is flagged as irrigated, falling back to thermal-only? Or should urban irrigation be modelled explicitly through the water balance and implicitly influence SMD? +- **Interaction with observed-LAI forcing (post GH-1291).** When the user supplies observed LAI the moisture-aware calculation is bypassed. This is likely correct; worth documenting explicitly. +- **tau window choice.** CLM5 uses 15 days; Dahlin et al. 2015 note sensitivity to this choice. Should it be per-surface and fittable, or class-default? +- **Rooting depth parameterisation.** The SMD capacity SUEWS uses is not the same as CLM5's root-zone soil-water potential. Does the SUEWS bucket depth need re-examining for dryland PFTs, or is the per-surface `SurfaceStoreCap` sufficient? + +--- + +## 11. Prerequisite and related work + +- **GH-1291** (observed-LAI forcing silently ignored, OPEN): the documented interim workaround for rainfall-driven sites does not work until fixed. Not a hard blocker for PR1 but must be resolved before PR3 calibration is meaningful. +- **Omidvar et al. 2022** Appendix D provides the motivating contrast and is the target for PR3 validation. +- **Existing FLUXNET2015 archive** under `/Users/tingsun/Dropbox (Personal)/6.Repos/SUEWS-FLUXNET2015/` and `/Users/tingsun/Dropbox (Personal)/6.Repos/LAI-model/` is the calibration and validation substrate. Read-only; copy required inputs into the active workspace rather than operating in place. + +--- + +## 12. References + +- Omidvar, H., Sun, T., Grimmond, S. et al. (2022). Surface Urban Energy and Water balance scheme (v2020a) in non-urban areas: developments, evaluation and impacts. *Geosci. Model Dev.* 15, 3041-3078. https://doi.org/10.5194/gmd-15-3041-2022 +- Lawrence, D. M. et al. (2019). The Community Land Model Version 5 (CLM5). *J. Adv. Model. Earth Syst.* 11, 4245-4287. https://doi.org/10.1029/2018MS001583 +- Sitch, S. et al. (2003). Evaluation of ecosystem dynamics, plant geography and terrestrial carbon cycling in the LPJ dynamic global vegetation model. *Glob. Change Biol.* 9, 161-185. https://doi.org/10.1046/j.1365-2486.2003.00569.x +- Smith, B. et al. (2014). Implications of incorporating N cycling and N limitations on primary production in an individual-based dynamic vegetation model. *Biogeosciences* 11, 2027-2054. https://doi.org/10.5194/bg-11-2027-2014 +- Krinner, G. et al. (2005). A dynamic global vegetation model for studies of the coupled atmosphere-biosphere system. *Glob. Biogeochem. Cycles* 19, GB1015. https://doi.org/10.1029/2003GB002199 +- Niu, G.-Y. et al. (2011). The community Noah land surface model with multiparameterization options (Noah-MP). *JGR Atmos.* 116, D12109. https://doi.org/10.1029/2010JD015139 +- Best, M. J. et al. (2011). The Joint UK Land Environment Simulator (JULES), model description - part 1. *GMD* 4, 677-699. https://doi.org/10.5194/gmd-4-677-2011 +- Clark, D. B. et al. (2011). JULES - part 2. *GMD* 4, 701-722. https://doi.org/10.5194/gmd-4-701-2011 +- Gibelin, A.-L. et al. (2006). ISBA-A-gs evaluation using interactive vegetation. *JGR Atmos.* 111, D18102. https://doi.org/10.1029/2005JD006691 +- Jolly, W. M., Nemani, R. & Running, S. W. (2005). A generalized, bioclimatic index to predict foliar phenology in response to climate. *Glob. Change Biol.* 11, 619-632. https://doi.org/10.1111/j.1365-2486.2005.00930.x +- Jolly, W. M. & Running, S. W. (2004). Effects of precipitation and soil water potential on drought deciduous phenology in the Kalahari. *Glob. Change Biol.* 10, 303-308. https://doi.org/10.1111/j.1365-2486.2003.00701.x +- Stockli, R. et al. (2008). Use of FLUXNET in the Community Land Model development. *JGR Biogeosci.* https://doi.org/10.1029/2007JG000562. And Stockli et al. (2011). A global reanalysis of vegetation phenology. *JGR Biogeosci.* 116, G03020. https://doi.org/10.1029/2010JG001545 +- Caldararu, S., Purves, D. W. & Palmer, P. I. (2014). Phenology as a strategy for carbon optimality: a global model. *Biogeosciences* 11, 763-778. https://doi.org/10.5194/bg-11-763-2014 +- Knorr, W. et al. (2010). Carbon cycle data assimilation with a generic phenology model. *JGR Biogeosci.* https://doi.org/10.1029/2009JG001119 +- Forkel, M. et al. (2014/2015). Codominant water control on global interannual variability and trends in land surface phenology and greenness. *Glob. Change Biol.* 21, 3414-3435. https://doi.org/10.1111/gcb.12950 +- Dahlin, K. M., Fisher, R. A. & Lawrence, P. J. (2015). Environmental drivers of drought deciduous phenology in the Community Land Model. *Biogeosciences* 12, 5061-5074. https://doi.org/10.5194/bg-12-5061-2015 +- Broich, M. et al. (2014). Land surface phenology performance across the semi-arid Australian continent. *Biogeosciences* 11, 5181-5198. https://doi.org/10.5194/bg-11-5181-2014 +- Zhang, X. et al. (2005). Monitoring the response of vegetation phenology to precipitation in Africa. *Int. J. Remote Sens.* 26, 4517-4533. https://doi.org/10.1080/01431160500113971 +- Liu, Y. et al. (2022). Green-up and brown-down: modelling grassland foliage phenology responses to soil moisture availability. *Agric. For. Meteorol.* 327, 109230. https://doi.org/10.1016/j.agrformet.2022.109230 +- Hickler, T. et al. (2005). Precipitation controls Sahel greening trend. *Geophys. Res. Lett.* 32, L21415. https://doi.org/10.1029/2005GL024370 +- Wang, J., Rich, P. M. & Price, K. P. (2003). Temporal responses of NDVI to precipitation and temperature in the central Great Plains, USA. *Int. J. Remote Sens.* 24, 2345-2364. https://doi.org/10.1080/01431160210154812 + +Internal references: + +- Upstream issue: https://github.com/UMEP-dev/SUEWS/issues/1292 +- Prerequisite: https://github.com/UMEP-dev/SUEWS/issues/1291 +- Developer workflow: `dev-ref/FEATURE_DEVELOPMENT_WORKFLOW.md` +- Repository style: `.claude/rules/00-project-essentials.md` +- Ting's prior FLUXNET2015 work: `/Users/tingsun/Dropbox (Personal)/6.Repos/SUEWS-FLUXNET2015/`, `/Users/tingsun/Dropbox (Personal)/6.Repos/LAI-model/`, `/Users/tingsun/Dropbox (Personal)/6.Repos/AmeriFluxEvaluation/` (read-only; see `ref_dropbox-legacy-repos` memory for handling rules) + +--- + +## Appendix A. Calibration methodology (PR3) + +The moisture-aware scheme exposes six parameters per vegetation surface: `w_wilt`, `w_opt`, `f_shape`, `w_on`, `w_off`, `tau_w`. PR3 ships the calibration **scaffolding**: a parameter-sweep tool that scans any one of these over a list of candidate values, reruns the sample, and reports RMSE vs an `LAIType = 0` baseline plus seasonal amplitude and green-up DOY. PR3 does not ship calibrated per-IGBP defaults -- that requires dedicated FLUXNET site configurations (lat/lon, land cover fractions, soil-store capacity, thermal GDD fits from Omidvar et al. 2022 Table S1), which live outside this repo. See A.3 below. + +### A.1 Sweep tool + +`scripts/verify/moisture_phenology_sweep.py` exposes two invocations: + +``` +# Scan one parameter with user-supplied values: +uv run python scripts/verify/moisture_phenology_sweep.py \ + --param w_opt --values 0.25 0.40 0.55 0.70 0.90 --dry-start + +# Scan all six with the bundled default ranges: +uv run python scripts/verify/moisture_phenology_sweep.py --all --dry-start +``` + +The `--dry-start` flag depletes `soilstore_surf` for the three vegetation surfaces to 10 mm on day 1 (the minimum the data model allows). Without it the bundled London sample is well-watered and the moisture gate never engages, so every sweep value collapses onto the thermal baseline. For a real dryland calibration the flag is unnecessary because the site is dry by construction. + +Outputs land in `.context/gh1292//sweep_.json` (tabular dose-response: mean LAI, RMSE vs baseline, seasonal amplitude, green-up DOY) and `.context/gh1292//sweep_.png` (dual-axis sensitivity plot). A pairwise-validator co-adjustment in the tool ensures `w_opt > w_wilt` and `w_on > w_off` remain satisfied during extreme sweep values. + +### A.2 London sensitivity results (PR3 demo, depleted soil initial state) + +From the `--all --dry-start` scan on the bundled London sample with default `LAIType = 0` baseline mean LAI = 4.546 m2/m2 (note London is a UK suburban site, not a dryland; the depleted soil gives a short "synthetic drought" pulse for ~30 days before winter recharge): + +- `w_wilt` strong effect: raising from 0.15 to 0.70 drops mean LAI from 4.545 to 4.435 (~2.4%) and delays green-up from DOY 104 to DOY 109. The Jarvis lower bound is the most influential parameter under transient drought. +- `w_opt` moderate: raising from 0.40 to 0.90 drops mean LAI from 4.545 to 4.470 (~1.6%) and delays green-up by five days. Tightening `w_opt` while keeping `w_wilt` unchanged narrows the Jarvis transition and can eclipse `w_wilt` for high values. +- `w_off` asymmetric: values below the default 0.20 are silent (the latch stays open); raising to 0.45 drops mean LAI 4.546 -> 4.483 (~1.4%) and lags green-up to DOY 106 by triggering the CLM5 off-latch early after the initial dry pulse. +- `w_on`, `tau_w`, `f_shape`: **no measurable effect** on the London synthetic-drought scenario with all other parameters at default. This is consistent with the latch mechanics: once the initial transient clears and the soil recharges, `wbar_id` sits comfortably above `w_off` and the leaf-on latch is never consulted again, so `w_on` / `tau_w` / `f_shape` stay off the integration. + +The practical implication for calibration: **start with `w_wilt`, `w_opt`, and `w_off`** as the sensitive trio; treat `w_on`, `tau_w`, and `f_shape` as fine-tuning parameters that only matter when the run actually experiences sustained drought transitions. Real dryland sites (AU-ASM, US-SRG, AU-DaS) will almost certainly exercise all six. + +Raw sweep outputs are checked into the `.context/` cache and can be regenerated with the command in A.1. PR3 does not check those outputs into the repo because they depend on the specific supy version and sample configuration. + +### A.3 Roadmap for real FLUXNET2015 calibration (out of scope for PR3) + +Full per-IGBP calibration requires: + +- Per-site SUEWS YAML configurations matching the FLUXNET2015 tower locations -- latitude, land-cover fractions, soil-store capacity, roughness, and the thermal fits already derived in `SUEWS-FLUXNET2015/data/csv-prm/prm-LAI-Albedo.csv`. These do not ship with SUEWS; recreating them from the Omidvar et al. (2022) pipeline is a half-day of work per site. +- MODIS LAI time series aligned to the simulation output at daily aggregation; the `h5-lsm/df_LSM_.h5` files in the Dropbox archive carry the MODIS fields already. +- An optimiser wrapper around the sweep tool. The natural choice is a bounded Nelder-Mead or differential-evolution run in the 6-dimensional parameter space minimising RMSE(MODIS, simulated) subject to the pairwise validators. Default IGBP priors from the literature (Sitch et al. 2003 raingreen `0.35`, CLM5 stress-deciduous `~15 d` persistence window) bracket the starting point. +- Primary tier (dryland + monsoon): AU-ASM, AU-DaS, AU-DaP, US-SRG, US-Wkg, US-Ton, AU-Stp. +- Control tier (non-regression at temperate sites): US-MMS, DE-Hai, US-Syv, DK-Sor. + +Acceptance bars (from section 8.3): dryland amplitude recovery >= 60 % of MODIS, monsoon green-up DOY error <= 21 days, temperate-site RMSE degradation <= 0.15 m2/m2, monthly flux (Qh, Qe) deviation <= 10 W/m2 outside the growing season. + +This calibration is genuinely scientific work, not engineering, and is recommended as a follow-up issue under the same GH-1292 umbrella (or a new issue that cites the fitted parameters). + +### A.4 Optional moisture-stress senescence term (deferred) + +Design C Section 5 option (d) proposes an **additive** SDD drought-stress term `delta_SDD_eff = delta_SDD - mu * max(0, w_off - wbar_id)` that accelerates senescence during sustained drought. In the CLM5-style latch implementation shipped with PR2, drought already shuts off thermal accumulation entirely via the latch, so the additive SDD term is scientifically optional rather than strictly necessary. PR3 does not add it because: + +- The PR2 latch captures the dominant effect (drought halts leaf-on) without new parameters. +- Adding the term would cascade through the bridge (LAI_PRM schema 2 -> 3, LC_*_PRM +1) and warrants its own review. +- Whether the additive term is needed is a calibration outcome: if the latch-only scheme under-predicts brown-down at monsoon sites (Liu et al. 2022 would suggest it might), add the term; if not, the cleaner single-mechanism scheme wins. + +The term is tracked as future work and will be revisited after A.3 calibration delivers empirical evidence. diff --git a/scripts/verify/moisture_phenology_drought_demo.py b/scripts/verify/moisture_phenology_drought_demo.py new file mode 100644 index 000000000..5f687e532 --- /dev/null +++ b/scripts/verify/moisture_phenology_drought_demo.py @@ -0,0 +1,424 @@ +"""GH-1292 PR3 drought demonstration figure (four-panel). + +Produces a single four-panel PNG that shows both control branches of the +Design C LAI phenology -- the moisture branch (rainfall -> SMD) and the +thermal branch (air temperature -> GDD / SDD accumulators) -- converging +on the LAI trajectory. + +Takes the bundled London sample (Ward et al. 2016 KCL benchmark), clones the forcing with rainfall zeroed +for a spring-through-early-summer window (so SUEWS's internal water +balance drains the root-zone store during the growth phase), runs the +model at ``LAIType = 0`` and ``LAIType = 2``, and overlays: + +* Panel A -- MOISTURE INPUT: daily rainfall (baseline vs drought + forcing). Shows the intervention directly. +* Panel B -- MOISTURE STATE: daily vegetation-surface SMD with the + Jarvis threshold markers (w_opt, w_wilt as SMD-equivalents on this + SoilStoreCap). The drought forcing drives SMD up toward capacity + during the zeroed-rainfall window. +* Panel C -- THERMAL STATE: daily mean air temperature on the left + axis together with the vegetation GDD and SDD accumulators on the + right axis. Shows that the thermal driver is identical between the + two LAIType runs -- any LAI divergence must come from the moisture + control. +* Panel D -- OUTCOME: simulated LAI under ``LAIType = 0`` (thermal + only) and ``LAIType = 2`` (moisture-aware). LAIType=2 flatlines at + LAImin during the drought window because the CLM5-style latch + prevents leaf-on while wbar_id sits below w_off. + +Output: ``.context/gh1292/drought-demo/moisture_control.png`` plus a +JSON summary. The figure is what gets embedded in the dashboard. + +Why synthetic rather than a real dryland forcing +------------------------------------------------ +A real FLUXNET dryland site (AU-ASM, US-SRG) would be the scientifically +correct test, but it requires a dedicated SUEWS configuration (site +geometry, land cover, thermal GDD fits) that does not ship with this +repo. The London sample config is UK-urban, so a synthetic drought +imposed through the rainfall forcing is the cleanest way to demonstrate +the moisture control within the shipped model without muddying the +message with a site-config mismatch. +""" + +from __future__ import annotations + +import json +from pathlib import Path + +import numpy as np +import pandas as pd + +OUT_DIR = Path(".context/gh1292/drought-demo") +LAI_TYPE_COLS = [("laitype", f"({i},)") for i in range(3)] +# Drought window placed EARLY (late winter through late spring) so it catches +# the spring leaf-out phase. London's thermal GDD saturates in mid-to-late +# spring; the drought blocks that. Ending the window in May (DOY ~130) rather +# than late June gives all three vegetation surfaces enough time to catch up +# through summer (they have staggered recovery because per-surface soil stores +# drain and refill at slightly different rates) and then senesce normally via +# SDD in autumn. Longer drought windows leave DecTr still in catch-up growth +# through October/November, at which point the cumulative-GDD scheme has no way +# of knowing it is calendar autumn -- the resulting late "growth" is a known +# cumulative-GDD-scheme limitation, not a moisture-control bug. +DROUGHT_DOY_START = 30 # 30 Jan +DROUGHT_DOY_END = 130 # 9 May + +# London's soil-store capacity is 150 mm, so the drought forcing pushes SMD +# to ~80 mm (w ~ 0.47). Default moisture thresholds (w_wilt=0.15, w_opt=0.40) +# would leave the Jarvis factor saturated at 1.0 throughout -- the scheme +# would degrade gracefully to LAIType=0, which is the correct physics but a +# boring demo. To make the moisture control *visible* we tighten the thresholds +# so the Jarvis curve engages in the actual w range the drought produces. +DEMO_MOISTURE_PARAMS = { + "w_wilt": 0.65, # below this -> f_w = 0 (complete shutdown) + "w_opt": 0.95, # above this -> f_w = 1 (no stress) + "f_shape": 1.0, + "w_on": 0.85, # leaf-on persistence threshold + "w_off": 0.60, # drop below -> CLM5 latch flips off + "tau_w": 10.0, +} + + +def _build_drought_forcing(df_forcing: pd.DataFrame) -> pd.DataFrame: + df_drought = df_forcing.copy() + # The forcing DataFrame index typically has datetime as level 1 (grid, datetime); + # support either shape defensively. + if isinstance(df_drought.index, pd.MultiIndex): + dt = df_drought.index.get_level_values("datetime") + else: + dt = df_drought.index + doy = dt.dayofyear + mask = (doy >= DROUGHT_DOY_START) & (doy <= DROUGHT_DOY_END) + if "rain" in df_drought.columns: + df_drought.loc[mask, "rain"] = 0.0 + return df_drought + + +def _run_scenario( + laitype: int, + df_forcing: pd.DataFrame, + moisture_params: dict | None = None, + dry_start: bool = False, +): + """Return (SUEWS, DailyState) output frames for the given laitype.""" + + from supy import SUEWSSimulation + + sim = SUEWSSimulation.from_sample_data() + state = sim.state_init.copy() + for col in LAI_TYPE_COLS: + state.loc[:, col] = laitype + if moisture_params and laitype == 2: + for name, value in moisture_params.items(): + for i in range(3): + state.loc[:, (name, f"({i},)")] = value + if dry_start: + # Pydantic enforces soilstore >= 10 mm; use 10 to start as close to + # full depletion as the data model allows. + for veg_idx in (2, 3, 4): + state.loc[:, ("soilstore_surf", f"({veg_idx},)")] = 10.0 + + scenario = SUEWSSimulation.from_state(state) + scenario.update_forcing(df_forcing) + output = scenario.run() + return output.SUEWS, output.DailyState + + +def _daily(series_index, values) -> pd.Series: + s = pd.Series(values, index=series_index) + if isinstance(s.index, pd.MultiIndex): + s = s.droplevel(0) + s = s[~s.index.duplicated(keep="first")] + return s.resample("1D").mean() + + +def _daily_sum(series_index, values) -> pd.Series: + s = pd.Series(values, index=series_index) + if isinstance(s.index, pd.MultiIndex): + s = s.droplevel(0) + s = s[~s.index.duplicated(keep="first")] + return s.resample("1D").sum() + + +def main() -> int: + from supy import SUEWSSimulation + + OUT_DIR.mkdir(parents=True, exist_ok=True) + + print("[drought-demo] Building baseline and drought forcings...") + sample = SUEWSSimulation.from_sample_data() + df_forcing_baseline = sample.forcing + # supy exposes forcing via a wrapper object with a DataFrame-like API. + # Pull the raw DataFrame so we can mutate the rain column deterministically. + df_baseline_raw = pd.DataFrame( + {col: df_forcing_baseline[col] for col in df_forcing_baseline.columns} + ) + df_drought_raw = _build_drought_forcing(df_baseline_raw) + + print("[drought-demo] Running LAIType=0 under drought forcing (thermal baseline, dry-start)...") + out_l0, daily_l0 = _run_scenario(0, df_drought_raw, dry_start=True) + print( + "[drought-demo] Running LAIType=2 under drought forcing " + f"with tightened demo parameters {DEMO_MOISTURE_PARAMS} + dry-start..." + ) + out_l2, daily_l2 = _run_scenario( + 2, df_drought_raw, moisture_params=DEMO_MOISTURE_PARAMS, dry_start=True + ) + + # Aggregate to daily series for plotting clarity. + dt_idx = out_l0.index + rain_daily_baseline = _daily_sum( + pd.MultiIndex.from_frame( + pd.DataFrame({"grid": [1] * len(df_baseline_raw), "datetime": df_baseline_raw.index}) + ) + if not isinstance(df_baseline_raw.index, pd.MultiIndex) + else df_baseline_raw.index, + df_baseline_raw["rain"].values, + ) + rain_daily_drought = _daily_sum( + pd.MultiIndex.from_frame( + pd.DataFrame({"grid": [1] * len(df_drought_raw), "datetime": df_drought_raw.index}) + ) + if not isinstance(df_drought_raw.index, pd.MultiIndex) + else df_drought_raw.index, + df_drought_raw["rain"].values, + ) + # SMD is identical under the same forcing regardless of LAIType -- use either run. + smd_grass = _daily(dt_idx, out_l0["SMDGrass"].values) + lai_l0 = _daily(dt_idx, out_l0["LAI"].values) + lai_l2 = _daily(dt_idx, out_l2["LAI"].values) + + # Thermal series (same for both LAIType runs since air-temperature forcing is identical). + tair_daily = _daily(dt_idx, out_l0["T2"].values) + # DailyState is emitted once per day at the end-of-day timestep; drop duplicated grid + # level and subset to the grass surface which is the representative vegetation type. + ds = daily_l0.copy() + if isinstance(ds.index, pd.MultiIndex): + ds = ds.droplevel(0) + ds = ds[~ds.index.duplicated(keep="first")] + gdd_grass = ds["GDD_Grass"].dropna() + sdd_grass = ds["SDD_Grass"].dropna() + + # --- Figure --------------------------------------------------------------- + import matplotlib + + matplotlib.use("Agg") + import matplotlib.pyplot as plt + import matplotlib.dates as mdates + + fig, axes = plt.subplots( + 4, 1, figsize=(11, 11.5), sharex=True, + gridspec_kw={"hspace": 0.13}, + ) + + def panel_label(ax, letter, description=None, x=-0.06, fontsize_solo=13, fontsize_titled=11): + """Panel label with left edge aligned to the y-axis label's left edge. + + Two scenarios, one alignment principle (``x`` always fixed): + + * **Solo** (``description=None``) -- the letter alone sits at the + upper-left corner of the panel, top edge on the axes top spine + (``va="top"``, ``y=1.0``). Use when the panel has no descriptive + title (e.g. stacked time-series with side labels). + * **Labelled title** (``description`` given) -- the letter is + concatenated with the description into one bold line placed above + the axes (``va="bottom"``, ``y=1.02``). Use when the panel's title + serves as the panel label. + + In both cases the left edge lines up with the y-axis label's left + edge (default ``x=-0.06`` in axes fraction); adjust per-figure if + the y-label sits at a different x. + """ + if description: + text = f"{letter}) {description}" + y, va, fs = 1.02, "bottom", fontsize_titled + else: + text = f"{letter})" + y, va, fs = 1.0, "top", fontsize_solo + ax.text( + x, y, text, + transform=ax.transAxes, + ha="left", va=va, + fontsize=fs, fontweight="bold", + color="#1a1a1a", + ) + + # Shared drought-window shading + def _shade(ax): + year = int(rain_daily_drought.index[0].year) + xmin = pd.Timestamp(year=year, month=1, day=1) + pd.Timedelta(days=DROUGHT_DOY_START - 1) + xmax = pd.Timestamp(year=year, month=1, day=1) + pd.Timedelta(days=DROUGHT_DOY_END - 1) + ax.axvspan(xmin, xmax, color="#d94922", alpha=0.06, lw=0) + return xmin, xmax + + # Section labels down the left-hand side for visual grouping. + def _side_label(ax, text, colour): + ax.text( + -0.075, 0.5, text, + transform=ax.transAxes, rotation=90, ha="center", va="center", + fontsize=9, color=colour, fontweight="bold", + ) + + # Panel A -- moisture input (rainfall) + ax = axes[0] + ax.plot( + rain_daily_baseline.index, rain_daily_baseline.values, + color="#5a9db5", lw=0.9, alpha=0.7, label="Baseline forcing", + ) + ax.plot( + rain_daily_drought.index, rain_daily_drought.values, + color="#d94922", lw=0.9, + label=f"Drought forcing (DOY {DROUGHT_DOY_START}-{DROUGHT_DOY_END} rainfall zeroed)", + ) + _shade(ax) + ax.set_ylabel("Daily rainfall\n[mm]") + ax.legend(loc="upper right", fontsize=9, frameon=False) + ax.grid(True, alpha=0.25) + + # Panel B -- moisture state (SMD + thresholds) + ax = axes[1] + ax.plot(smd_grass.index, smd_grass.values, color="#b06e21", lw=1.2, label="SMD (grass surface)") + _shade(ax) + ax.set_ylabel("SMD [mm]\n(drier upward)") + smdcap = 150.0 + smd_w_opt = smdcap * (1.0 - DEMO_MOISTURE_PARAMS["w_opt"]) + smd_w_wilt = smdcap * (1.0 - DEMO_MOISTURE_PARAMS["w_wilt"]) + ax.axhline(smd_w_opt, color="#6da33a", lw=0.8, ls="--", alpha=0.7) + ax.axhline(smd_w_wilt, color="#d94922", lw=0.8, ls="--", alpha=0.7) + ax.text( + rain_daily_drought.index[-40], smd_w_opt + 1, + f"w_opt = {DEMO_MOISTURE_PARAMS['w_opt']} (SMD ~ {smd_w_opt:.0f} mm)", + color="#6da33a", fontsize=8, va="bottom", ha="right", + ) + ax.text( + rain_daily_drought.index[-40], smd_w_wilt + 1, + f"w_wilt = {DEMO_MOISTURE_PARAMS['w_wilt']} (SMD ~ {smd_w_wilt:.0f} mm)", + color="#d94922", fontsize=8, va="bottom", ha="right", + ) + ax.legend(loc="upper left", fontsize=9, frameon=False) + ax.grid(True, alpha=0.25) + + # Panel C -- thermal state (air temperature + GDD / SDD on twin axis) + ax = axes[2] + ax.plot(tair_daily.index, tair_daily.values, + color="#5a9db5", lw=1.0, alpha=0.8, label="Daily mean air T (T2)") + _shade(ax) + ax.set_ylabel("Air temperature\n[\u00b0C]", color="#5a9db5") + ax.tick_params(axis="y", labelcolor="#5a9db5") + ax.axhline(0, color="#5a9db5", lw=0.4, alpha=0.35) + ax.grid(True, alpha=0.25) + + ax_right = ax.twinx() + ax_right.plot(gdd_grass.index, gdd_grass.values, + color="#6da33a", lw=1.3, label="GDD (grass)") + ax_right.plot(sdd_grass.index, sdd_grass.values, + color="#a34b6d", lw=1.3, label="SDD (grass)") + ax_right.set_ylabel("GDD / SDD\n[\u00b0C day]", color="#4a3a2a") + ax_right.axhline(0, color="#6c6456", lw=0.4, alpha=0.35) + lines1, labels1 = ax.get_legend_handles_labels() + lines2, labels2 = ax_right.get_legend_handles_labels() + ax.legend(lines1 + lines2, labels1 + labels2, loc="upper left", fontsize=9, frameon=False, ncol=3) + + # Panel D -- outcome (LAI) with phase annotations + ax = axes[3] + ax.plot(lai_l0.index, lai_l0.values, color="#3a342d", lw=1.5, label="LAIType=0 (thermal only)") + ax.plot(lai_l2.index, lai_l2.values, color="#c4841d", lw=1.5, label="LAIType=2 (moisture-aware)") + _shade(ax) + ax.set_ylabel("LAI\n[m$^2$ m$^{-2}$]") + ax.set_xlabel("Date") + ax.legend(loc="lower center", fontsize=9, frameon=False, ncol=2) + ax.grid(True, alpha=0.25) + ax.xaxis.set_major_locator(mdates.MonthLocator(bymonth=(1, 3, 5, 7, 9, 11))) + ax.xaxis.set_major_formatter(mdates.DateFormatter("%b")) + + # --- Phase band on Panel D --- + # Keep the top band (visual anchor); the per-phase rationale lives in the + # caption (see dashboard.html) so the figure itself stays clean. + year = int(rain_daily_drought.index[0].year) + def _ts(doy): + return pd.Timestamp(year=year, month=1, day=1) + pd.Timedelta(days=doy - 1) + + # Four phases, boundaries picked from the simulation itself: + # P1 winter (DOY 1 - 95): T below BaseT_GDD, no leaf-on from either scheme + # P2 spring-drought (95 - DROUGHT_DOY_END): GDD crosses BaseT_GDD but moisture latch closed + # P3 recovery (DROUGHT_DOY_END - 220): rain returns, latch reopens, LAI2 catches up + # P4 summer -> autumn (220 - end): both at LAImax, SDD drives senescence + phase_bounds = [1, 95, DROUGHT_DOY_END, 220, 365] + phase_short = [ + "P1 winter", + "P2 spring drought", + "P3 recovery", + "P4 summer -> autumn", + ] + phase_colours = ["#5a9db5", "#d94922", "#6da33a", "#a07050"] + y_lo, y_hi = ax.get_ylim() + # Make room for the label band above the data. + ax.set_ylim(y_lo, y_hi + 0.95) + label_y = y_hi + 0.32 + band_y = y_hi + 0.78 + + # Vertical dotted boundaries behind the data. + for doy in phase_bounds[1:-1]: + ax.axvline(_ts(doy), color="#6c6456", lw=0.6, ls=":", alpha=0.7, zorder=0) + + # Coloured phase band at the top with short PX labels. + for (lo, hi), short, col in zip( + zip(phase_bounds[:-1], phase_bounds[1:]), phase_short, phase_colours + ): + ax.plot([_ts(lo), _ts(hi)], [band_y, band_y], color=col, lw=5, solid_capstyle="butt") + ax.text( + _ts((lo + hi) / 2), label_y, short, + ha="center", va="bottom", fontsize=9, + color=col, fontweight="bold", + ) + + # Section labels down the left-hand side for visual grouping. + _side_label(axes[0], "MOISTURE INPUT", "#5a9db5") + _side_label(axes[1], "MOISTURE STATE", "#b06e21") + _side_label(axes[2], "THERMAL STATE", "#6da33a") + _side_label(axes[3], "LAI OUTCOME", "#c4841d") + + # Panel labels A)/B)/C)/D) at the upper-left corner of each panel, outside + # the axes (aligned roughly with the y-label's left edge). Matches the + # convention used by moisture_phenology_sensitivity_summary.py. + for letter, axN in zip("ABCD", axes): + panel_label(axN, letter) + + # Figure-level caption -- describes the whole four-panel composition. + fig.suptitle( + "GH-1292 moisture-aware LAI phenology: two control branches on the London sample", + fontsize=11, color="#2a2a2a", x=0.06, y=0.995, ha="left", + ) + + # Tight margins so panel letters hug their panels rather than floating in + # the inter-panel gap. + fig.subplots_adjust(left=0.13, right=0.95, top=0.96, bottom=0.06, hspace=0.13) + png_path = OUT_DIR / "moisture_control.png" + fig.savefig(png_path, dpi=150, facecolor="white") + plt.close(fig) + + # --- Numeric summary ------------------------------------------------------ + drought_mask = (lai_l0.index.dayofyear >= DROUGHT_DOY_START) & ( + lai_l0.index.dayofyear <= DROUGHT_DOY_END + ) + summary = { + "drought_window_doy": [DROUGHT_DOY_START, DROUGHT_DOY_END], + "smd_mean_outside_drought_mm": float(smd_grass[~drought_mask].mean()), + "smd_mean_inside_drought_mm": float(smd_grass[drought_mask].mean()), + "lai_mean_laitype0_inside_drought": float(lai_l0[drought_mask].mean()), + "lai_mean_laitype2_inside_drought": float(lai_l2[drought_mask].mean()), + "lai_mean_diff_inside_drought": float(lai_l0[drought_mask].mean() - lai_l2[drought_mask].mean()), + "lai_max_laitype0": float(lai_l0.max()), + "lai_max_laitype2": float(lai_l2.max()), + } + (OUT_DIR / "moisture_control.json").write_text( + json.dumps(summary, indent=2), encoding="utf-8" + ) + + print(f"[drought-demo] Wrote {png_path}") + print(f"[drought-demo] Summary: {json.dumps(summary, indent=2)}") + return 0 + + +if __name__ == "__main__": + raise SystemExit(main()) diff --git a/scripts/verify/moisture_phenology_sensitivity_summary.py b/scripts/verify/moisture_phenology_sensitivity_summary.py new file mode 100644 index 000000000..fd017b273 --- /dev/null +++ b/scripts/verify/moisture_phenology_sensitivity_summary.py @@ -0,0 +1,222 @@ +"""GH-1292 PR3 consolidated sensitivity figure. + +Replaces the six-panel sweep grid with a single, self-explanatory +figure that answers one question: *which of the six moisture-aware +LAI parameters actually move the needle on the dry-start London +sample, and in what direction?* + +Design (figure-craft principles) +-------------------------------- +* One clear message: rank the parameters by sensitivity. +* Two stacked panels with a shared narrative: + - Top: horizontal bar chart of |mean LAI deviation from baseline| + across the parameter's sweep range -- a ranking at a glance. + - Bottom: dose-response lines for the three parameters that matter. + X-axis is the normalised sweep position (0 = lowest value in the + sweep, 1 = highest), labelled with the actual numeric values + underneath. +* Legends pulled outside the plot area; no overlap with data. +* Consistent typography and muted palette tied to the dashboard. + +Expects JSON artefacts under ``.context/gh1292//sweep_.json`` +produced by ``moisture_phenology_sweep.py --all --dry-start``. +""" + +from __future__ import annotations + +import json +from pathlib import Path +from typing import Dict, List, Tuple + +import matplotlib +import matplotlib.pyplot as plt +import numpy as np + +matplotlib.use("Agg") + +SITE = "London" +OUT_DIR = Path(".context/gh1292") / SITE +PARAM_ORDER = ["w_wilt", "w_opt", "w_off", "w_on", "f_shape", "tau_w"] +PARAM_LABELS = { + "w_wilt": "w_wilt", + "w_opt": "w_opt", + "w_off": "w_off", + "w_on": "w_on", + "f_shape": "f_shape", + "tau_w": "tau_w", +} + + +def _load_sweeps() -> Tuple[Dict[str, List[dict]], float]: + sweeps: Dict[str, List[dict]] = {} + baseline_mean = None + for p in PARAM_ORDER: + path = OUT_DIR / f"sweep_{p}.json" + if not path.exists(): + raise FileNotFoundError( + f"Missing sweep artefact {path}. " + "Run `uv run python scripts/verify/moisture_phenology_sweep.py --all --dry-start` first." + ) + payload = json.loads(path.read_text(encoding="utf-8")) + sweeps[p] = payload["points"] + baseline_mean = float(payload["baseline_mean_lai"]) + assert baseline_mean is not None + return sweeps, baseline_mean + + +def _sensitivity(points: List[dict], baseline_mean: float) -> float: + """Return max |mean LAI - baseline| over the sweep, i.e. worst-case deviation.""" + devs = [abs(pt["mean_lai"] - baseline_mean) for pt in points] + return max(devs) if devs else 0.0 + + +def main() -> int: + sweeps, baseline_mean = _load_sweeps() + sensitivities = {p: _sensitivity(sweeps[p], baseline_mean) for p in PARAM_ORDER} + ranked = sorted(sensitivities.items(), key=lambda kv: -kv[1]) + + # Split into sensitive vs silent using a fixed threshold (0.02 m2 m-2 ≈ RMSE noise floor). + THRESHOLD = 0.02 + sensitive = [p for p, s in ranked if s >= THRESHOLD] + silent = [p for p, s in ranked if s < THRESHOLD] + + fig = plt.figure(figsize=(10.5, 7.2), facecolor="white") + gs = fig.add_gridspec( + 2, 1, height_ratios=[1.1, 1.6], hspace=0.32, left=0.13, right=0.83, top=0.92, bottom=0.10, + ) + + def panel_label(ax, letter, description=None, x=-0.06, fontsize_solo=13, fontsize_titled=11): + """Panel label with left edge aligned to the y-axis label's left edge. + + Two scenarios, one alignment principle (``x`` always fixed): + + * **Solo** (``description=None``) -- letter alone at the upper-left + corner, top edge on the axes top spine (``va="top"``, ``y=1.0``). + Use when the panel has no descriptive title. + * **Labelled title** (``description`` given) -- letter concatenated + with the description into one bold line placed above the axes + (``va="bottom"``, ``y=1.02``). Use when the panel's title + doubles as the panel label. + + Shared signature with moisture_phenology_drought_demo.py so the + two figures look identical under both scenarios. + """ + if description: + text = f"{letter}) {description}" + y, va, fs = 1.02, "bottom", fontsize_titled + else: + text = f"{letter})" + y, va, fs = 1.0, "top", fontsize_solo + ax.text( + x, y, text, + transform=ax.transAxes, + ha="left", va=va, + fontsize=fs, fontweight="bold", + color="#1a1a1a", + ) + + # ---------------- Panel A: sensitivity ranking bar chart ---------------- + ax_bar = fig.add_subplot(gs[0]) + names = [p for p, _ in ranked] + vals = [s for _, s in ranked] + colours = ["#c4841d" if s >= THRESHOLD else "#b7ae9e" for _, s in ranked] + bars = ax_bar.barh( + range(len(names)), vals, color=colours, height=0.6, edgecolor="none", + ) + ax_bar.set_yticks(range(len(names))) + ax_bar.set_yticklabels( + [PARAM_LABELS[p] for p in names], fontsize=10, fontfamily="monospace", + ) + ax_bar.invert_yaxis() + ax_bar.set_xlabel( + "max |mean LAI deviation from baseline| [m$^2$ m$^{-2}$]", + fontsize=9, + ) + panel_label( + ax_bar, "A", + description="Parameter sensitivity ranking (London dry-start drought sweep)", + ) + ax_bar.axvline(THRESHOLD, color="#6c6456", lw=0.8, ls=":", alpha=0.8) + ax_bar.text( + THRESHOLD, len(names) - 0.25, + f" noise floor ({THRESHOLD:.2f})", + fontsize=8, color="#6c6456", va="center", ha="left", + ) + for bar, val in zip(bars, vals): + ax_bar.text( + val + 0.005, bar.get_y() + bar.get_height() / 2, + f"{val:.3f}", + va="center", ha="left", fontsize=8.5, color="#2a2a2a", + ) + ax_bar.spines["top"].set_visible(False) + ax_bar.spines["right"].set_visible(False) + ax_bar.grid(True, axis="x", alpha=0.25) + + # Panel-level legend (outside plot, right-hand gutter). + legend_ax = fig.add_axes([0.845, 0.63, 0.14, 0.22]) + legend_ax.axis("off") + legend_ax.add_patch(plt.Rectangle((0.05, 0.60), 0.18, 0.08, color="#c4841d")) + legend_ax.text(0.28, 0.64, "sensitive", fontsize=9, va="center") + legend_ax.add_patch(plt.Rectangle((0.05, 0.42), 0.18, 0.08, color="#b7ae9e")) + legend_ax.text(0.28, 0.46, "silent", fontsize=9, va="center") + legend_ax.text( + 0.05, 0.28, + f"n = {len(sensitive)} above\nnoise, {len(silent)} below.", + fontsize=8, color="#6c6456", va="top", + ) + legend_ax.text( + 0.05, 0.08, + f"baseline mean LAI\n{baseline_mean:.3f} m$^2$ m$^{{-2}}$", + fontsize=7.5, color="#6c6456", va="top", style="italic", + ) + + # ---------------- Panel B: dose-response lines for the sensitive trio ---- + ax = fig.add_subplot(gs[1]) + palette = {"w_wilt": "#d94922", "w_opt": "#c4841d", "w_off": "#8a4a8a", "w_on": "#5a9db5", "f_shape": "#6da33a", "tau_w": "#a07050"} + # Only plot the sensitive ones in full colour; draw the silent set as faint grey references. + for p in PARAM_ORDER: + pts = sweeps[p] + xs = np.linspace(0, 1, len(pts)) + ys = [pt["mean_lai"] for pt in pts] + is_sensitive = sensitivities[p] >= THRESHOLD + col = palette[p] if is_sensitive else "#d0c9b8" + lw = 1.8 if is_sensitive else 1.0 + alpha = 1.0 if is_sensitive else 0.55 + ax.plot(xs, ys, "o-", color=col, lw=lw, markersize=5, label=PARAM_LABELS[p], alpha=alpha) + + ax.axhline( + baseline_mean, color="#3a342d", lw=0.9, ls="--", alpha=0.75, + label="LAIType=0 baseline", + ) + ax.set_xticks([0, 0.25, 0.5, 0.75, 1.0]) + ax.set_xticklabels(["min", "", "mid", "", "max"], fontsize=9) + ax.set_xlabel( + "sweep position (per-parameter value range normalised, 0 = low end, 1 = high end)", + fontsize=9, + ) + ax.set_ylabel("mean annual LAI [m$^2$ m$^{-2}$]", fontsize=9) + panel_label( + ax, "B", + description="Dose-response across the sweep: tightening the sensitive trio pulls mean LAI below baseline", + ) + ax.grid(True, alpha=0.25) + ax.spines["top"].set_visible(False) + ax.spines["right"].set_visible(False) + + # Legend for dose-response (outside right). + ax.legend( + loc="center left", bbox_to_anchor=(1.02, 0.5), + fontsize=9, frameon=False, borderaxespad=0.0, + ) + + out_path = OUT_DIR / "sensitivity_summary.png" + fig.savefig(out_path, dpi=160, bbox_inches="tight", facecolor="white") + plt.close(fig) + print(f"[sensitivity-summary] Wrote {out_path}") + print(f"[sensitivity-summary] Sensitive parameters: {sensitive}") + print(f"[sensitivity-summary] Silent parameters : {silent}") + return 0 + + +if __name__ == "__main__": + raise SystemExit(main()) diff --git a/scripts/verify/moisture_phenology_site.py b/scripts/verify/moisture_phenology_site.py new file mode 100644 index 000000000..90bea5659 --- /dev/null +++ b/scripts/verify/moisture_phenology_site.py @@ -0,0 +1,282 @@ +"""GH-1292 forcing-driven diagnostic harness for moisture-aware LAI. + +Usage +----- + uv run python scripts/verify/moisture_phenology_site.py --site AU-ASM + +Purpose +------- +Runs SUEWS under the bundled sample configuration while swapping in a +FLUXNET2015 forcing series for all three vegetation ``LAItype`` values +(0 = thermal original, 1 = thermal high-latitude, 2 = moisture-aware) +and writes a diagnostic under +``.context/gh1292//``: + +* ``lai_timeseries.png`` -- simulated LAI per variant, plus the observed + MODIS LAI carried in the forcing file's ``lai`` column. +* ``metrics.json`` -- RMSE (vs MODIS), seasonal amplitude, green-up and + brown-down day-of-year estimates per variant and per vegetation + surface. +* ``summary.txt`` -- one-line verdict for the forcing-driven diagnostic. + The forcing comes from the requested FLUXNET site, but the static + site/state configuration remains the bundled sample unless the script + grows explicit per-site configs in a future follow-up. + +The forcing text file is read from Ting's local FLUXNET2015 archive +(``/Users/tingsun/Dropbox (Personal)/6.Repos/SUEWS-FLUXNET2015/``, +read-only per the ``ref_dropbox-legacy-repos`` memory); the script +copies the file into ``.context/gh1292//`` on first run and uses +the cached copy thereafter. + +This script is intentionally kept out of pytest: it depends on external +data that does not ship with the repo and is meant to be run by hand. +""" + +from __future__ import annotations + +import argparse +import json +import shutil +import sys +from pathlib import Path +from typing import Dict + +import numpy as np +import pandas as pd + +FLUXNET_ARCHIVE = Path( + "/Users/tingsun/Dropbox (Personal)/6.Repos/SUEWS-FLUXNET2015/data" +) +DEFAULT_SITE = "AU-ASM" +LAI_TYPE_VARIANTS = [0, 1, 2] +LAI_TYPE_COLS = [ + ("laitype", "(0,)"), + ("laitype", "(1,)"), + ("laitype", "(2,)"), +] +VEG_SURFACE_NAMES = ["evetr", "dectr", "grass"] + + +def parse_args() -> argparse.Namespace: + parser = argparse.ArgumentParser(description=__doc__.splitlines()[0]) + parser.add_argument( + "--site", + default=DEFAULT_SITE, + help="FLUXNET2015 forcing-site identifier (e.g. AU-ASM, AU-DaS, US-MMS)", + ) + parser.add_argument( + "--year", + type=int, + default=None, + help="Optional year filter; default uses the first full year in the forcing record", + ) + parser.add_argument( + "--archive", + type=Path, + default=FLUXNET_ARCHIVE, + help="FLUXNET2015 archive root (defaults to Ting's local Dropbox path)", + ) + return parser.parse_args() + + +def ensure_cached_forcing(site: str, archive: Path) -> Path: + """Copy the FLUXNET forcing text file into .context/gh1292// if not cached.""" + + source = archive / "suews_forcing" / f"FLUXNET2015_LST_{site}.txt" + if not source.exists(): + raise FileNotFoundError( + f"Expected FLUXNET forcing file missing: {source}" + ) + cache_dir = Path(".context/gh1292") / site + cache_dir.mkdir(parents=True, exist_ok=True) + dest = cache_dir / source.name + if not dest.exists(): + shutil.copy2(source, dest) + return dest + + +def _datetime_index(df: pd.DataFrame) -> pd.DatetimeIndex: + """Return the datetime axis regardless of whether forcing uses a flat or multi-index.""" + + if isinstance(df.index, pd.MultiIndex): + return df.index.get_level_values("datetime") + return pd.DatetimeIndex(df.index) + + +def _infer_first_full_year(dt_index: pd.DatetimeIndex) -> int: + """Pick the first full calendar year, falling back to the first year present.""" + + years = sorted(dt_index.year.unique()) + for year in years: + year_index = dt_index[dt_index.year == year] + if year_index.empty: + continue + if ( + year_index.min() <= pd.Timestamp(year=year, month=1, day=1) + and year_index.max() >= pd.Timestamp(year=year, month=12, day=31) + ): + return year + return int(years[0]) + + +def load_forcing(path: Path, year: int | None) -> pd.DataFrame: + """Read the FLUXNET forcing file using SuPy's datetime-aware parser.""" + + from supy.util import read_forcing + + df = read_forcing(str(path), tstep_mod=300).copy() + dt_index = _datetime_index(df) + target_year = year if year is not None else _infer_first_full_year(dt_index) + df = df[dt_index.year == target_year] + if df.empty: + raise ValueError( + "forcing file has no rows for the requested year; check --year value" + ) + return df + + +def run_scenario( + laitype: int, + forcing: pd.DataFrame, + start_date: pd.Timestamp, + end_date: pd.Timestamp, +) -> pd.DataFrame: + """Run SUEWS with the bundled sample state and the requested external forcing.""" + + # Imported lazily so this script can surface a clean error if supy is missing. + from supy import SUEWSSimulation + + sim = SUEWSSimulation.from_sample_data() + state = sim.state_init.copy() + for col in LAI_TYPE_COLS: + state.loc[:, col] = laitype + scenario = SUEWSSimulation.from_state(state) + scenario.update_forcing(forcing) + output = scenario.run(start_date=start_date, end_date=end_date) + df = output.SUEWS[["LAI"]].copy() + df.columns = ["LAI"] + return df + + +def compute_metrics(lai_sim: pd.Series, lai_obs: pd.Series) -> Dict[str, float]: + """Compute RMSE vs observed and seasonal amplitude / onset / offset DOY.""" + + aligned = pd.concat([lai_sim, lai_obs], axis=1).dropna() + rmse = float("nan") + if not aligned.empty: + sim_aligned = aligned.iloc[:, 0] + obs = aligned.iloc[:, 1] + rmse = float(np.sqrt(((sim_aligned - obs) ** 2).mean())) + + sim = lai_sim.dropna() + daily = sim.resample("1D").mean() + if daily.empty: + return { + "rmse_vs_obs": rmse, + "amplitude": float("nan"), + "green_up_doy": float("nan"), + "brown_down_doy": float("nan"), + } + peak = float(daily.max()) + trough = float(daily.min()) + amplitude = peak - trough + threshold = trough + 0.5 * amplitude + above = daily >= threshold + if above.any(): + green_doy = int(above.idxmax().dayofyear) + brown_doy = int(above[::-1].idxmax().dayofyear) + else: + green_doy = float("nan") + brown_doy = float("nan") + return { + "rmse_vs_obs": rmse, + "amplitude": amplitude, + "green_up_doy": float(green_doy), + "brown_down_doy": float(brown_doy), + } + + +def main() -> int: + args = parse_args() + archive = args.archive + forcing_path = ensure_cached_forcing(args.site, archive) + forcing = load_forcing(forcing_path, args.year) + start_date = forcing.index.min() + end_date = forcing.index.max() + obs_lai = forcing["lai"].replace(-999.0, np.nan) + obs_lai = obs_lai[~obs_lai.index.duplicated(keep="first")].dropna() + + results: Dict[str, pd.Series] = {} + metrics: Dict[str, Dict[str, float]] = {} + for laitype in LAI_TYPE_VARIANTS: + lai_df = run_scenario(laitype, forcing, start_date, end_date) + lai_series = lai_df["LAI"] + # Drop pandas multi-index grid level if present so we can align with observations. + if isinstance(lai_series.index, pd.MultiIndex): + lai_series = lai_series.droplevel(0) + lai_series = lai_series[~lai_series.index.duplicated(keep="first")] + results[f"laitype_{laitype}"] = lai_series + metrics[f"laitype_{laitype}"] = compute_metrics(lai_series, obs_lai) + + out_dir = Path(".context/gh1292") / args.site + out_dir.mkdir(parents=True, exist_ok=True) + + metrics["observed"] = { + "rmse_vs_obs": 0.0, + "amplitude": float(obs_lai.max() - obs_lai.min()), + "green_up_doy": float("nan"), + "brown_down_doy": float("nan"), + } + (out_dir / "metrics.json").write_text( + json.dumps(metrics, indent=2), encoding="utf-8" + ) + + try: + import matplotlib + + matplotlib.use("Agg") + import matplotlib.pyplot as plt + + fig, ax = plt.subplots(figsize=(10, 5)) + obs_daily = obs_lai.resample("1D").mean() + ax.plot(obs_daily.index, obs_daily.values, label="Observed (MODIS LAI)", color="k", lw=1.2) + colours = {0: "tab:blue", 1: "tab:orange", 2: "tab:green"} + for laitype in LAI_TYPE_VARIANTS: + series = results[f"laitype_{laitype}"].resample("1D").mean() + ax.plot( + series.index, + series.values, + label=f"Simulated (LAIType={laitype})", + color=colours[laitype], + lw=0.9, + alpha=0.8, + ) + ax.set_xlabel("Date") + ax.set_ylabel("LAI") + ax.set_title( + f"GH-1292 moisture-aware forcing diagnostic -- {args.site} forcing" + ) + ax.legend(loc="upper right") + fig.tight_layout() + fig.savefig(out_dir / "lai_timeseries.png", dpi=160) + plt.close(fig) + except ImportError: + print("matplotlib not available; skipping plot", file=sys.stderr) + + lai0 = results["laitype_0"] + lai2 = results["laitype_2"] + aligned = pd.concat([lai0, lai2], axis=1, join="inner").dropna() + max_abs_diff = float(np.max(np.abs(aligned.iloc[:, 0].values - aligned.iloc[:, 1].values))) + verdict = ( + "Forcing diagnostic OK: LAIType=2 reproduces LAIType=0 for this forcing." + if max_abs_diff == 0.0 + else f"Forcing diagnostic diverged: max |LAI0 - LAI2| = {max_abs_diff:.6e}." + ) + (out_dir / "summary.txt").write_text(verdict + "\n", encoding="utf-8") + print(f"[{args.site}] {verdict}") + print(f"Outputs under: {out_dir}") + return 0 + + +if __name__ == "__main__": + sys.exit(main()) diff --git a/scripts/verify/moisture_phenology_sweep.py b/scripts/verify/moisture_phenology_sweep.py new file mode 100644 index 000000000..c72c57c60 --- /dev/null +++ b/scripts/verify/moisture_phenology_sweep.py @@ -0,0 +1,310 @@ +"""GH-1292 PR3 parameter-sweep tool for moisture-aware LAI calibration. + +Usage +----- + uv run python scripts/verify/moisture_phenology_sweep.py \\ + --param w_opt --values 0.3 0.5 0.7 0.9 --dry-start + + # All-parameter sensitivity scan at the bundled sample with depleted soil: + uv run python scripts/verify/moisture_phenology_sweep.py --all --dry-start + +Purpose +------- +Scans one of the six moisture-aware ``LAI_PRM`` parameters (``w_wilt``, +``w_opt``, ``f_shape``, ``w_on``, ``w_off``, ``tau_w``) over a +user-provided value list, runs SUEWS for each, and emits a dose-response +summary alongside the existing ``LAIType = 0`` baseline. Outputs go to +``.context/gh1292//sweep_.{json,png}`` and are the +foundation for eventual FLUXNET-driven calibration. + +The tool is a calibration **scaffold**, not a full parameter fit: + +* Single-site, single-parameter-per-invocation (serialised; nested + scans are available via ``--all`` or repeated calls). +* Uses the bundled London sample by default; ``--dry-start`` depletes + the vegetation soil store on day 1 so the moisture gate is actually + engaged during the year (London is otherwise well-watered). +* Reports RMSE and mean-difference against the ``LAIType = 0`` + baseline, plus seasonal amplitude and estimated green-up DOY. + +A proper per-IGBP calibration against FLUXNET2015 requires dedicated +per-site SUEWS configurations (latitude, land-cover fractions, +soil-store capacity, thermal fits from Omidvar et al. 2022 Table S1). +Those configs are not shipped with this tool; see the "Roadmap for +full calibration" appendix in +``dev-ref/design-notes/gh1292-moisture-phenology.md``. +""" + +from __future__ import annotations + +import argparse +import json +import sys +from dataclasses import dataclass +from pathlib import Path +from typing import Dict, Iterable, List, Optional, Sequence + +import numpy as np +import pandas as pd + +DEFAULT_SITE = "London" # bundled sample is central London (lat 51.51, lng -0.12; Ward et al. 2016 KCL benchmark) +PARAM_NAMES = ("w_wilt", "w_opt", "f_shape", "w_on", "w_off", "tau_w") +DEFAULT_SWEEP_VALUES: Dict[str, Sequence[float]] = { + "w_wilt": (0.05, 0.15, 0.30, 0.50, 0.70), + "w_opt": (0.25, 0.40, 0.55, 0.70, 0.90), + "f_shape": (0.5, 1.0, 1.5, 2.0, 3.0), + "w_on": (0.15, 0.30, 0.45, 0.60, 0.80), + "w_off": (0.05, 0.10, 0.20, 0.30, 0.45), + "tau_w": (3.0, 7.0, 15.0, 30.0, 60.0), +} +LAI_TYPE_COLS = [("laitype", f"({i},)") for i in range(3)] + + +@dataclass(frozen=True) +class SweepResult: + param: str + value: float + mean_lai: float + rmse_vs_baseline: float + amplitude: float + green_up_doy: Optional[float] + + +def parse_args() -> argparse.Namespace: + parser = argparse.ArgumentParser(description=__doc__.splitlines()[0]) + parser.add_argument( + "--param", + choices=PARAM_NAMES, + help="Single parameter to sweep (overrides --all)", + ) + parser.add_argument( + "--values", + type=float, + nargs="+", + help="Override default sweep values for the selected parameter", + ) + parser.add_argument( + "--all", + action="store_true", + help="Sweep every parameter using default value ranges", + ) + parser.add_argument( + "--site", + default=DEFAULT_SITE, + help="Site tag used in output paths (default London, i.e. the bundled sample based on the KCL King's Cross benchmark)", + ) + parser.add_argument( + "--dry-start", + action="store_true", + help=( + "Deplete soilstore_surf on vegetation indices (2, 3, 4) to 10 mm on " + "day 1 so the moisture gate engages at well-watered sites" + ), + ) + args = parser.parse_args() + if not args.all and not args.param: + parser.error("either --param or --all is required") + return args + + +def _apply_dry_start(state: pd.DataFrame) -> pd.DataFrame: + for veg_idx in (2, 3, 4): + state.loc[:, ("soilstore_surf", f"({veg_idx},)")] = 10.0 + return state + + +def _set_laitype(state: pd.DataFrame, value: int) -> pd.DataFrame: + for col in LAI_TYPE_COLS: + state.loc[:, col] = value + return state + + +def _set_param(state: pd.DataFrame, param: str, value: float) -> pd.DataFrame: + for i in range(3): + state.loc[:, (param, f"({i},)")] = value + return state + + +def _co_adjust(param_overrides: Dict[str, float]) -> Dict[str, float]: + """Preserve the LAIParams validators (w_opt > w_wilt, w_on > w_off) during a sweep. + + When sweeping ``w_wilt`` we may exceed the default ``w_opt`` (0.40); similarly + when sweeping ``w_off`` we may exceed the default ``w_on`` (0.35). The validator + fires under laitype=2 (as designed), so we lift the paired threshold by a small + margin whenever the swept value would violate the constraint. + """ + + adjusted = dict(param_overrides) + margin = 0.05 + # w_wilt sweep: lift w_opt if necessary so w_opt > w_wilt. + if "w_wilt" in adjusted and adjusted["w_wilt"] + margin >= 0.40: + adjusted.setdefault("w_opt", min(0.98, adjusted["w_wilt"] + margin)) + # w_off sweep: lift w_on if necessary so w_on > w_off. + if "w_off" in adjusted and adjusted["w_off"] + margin >= 0.35: + adjusted.setdefault("w_on", min(0.98, adjusted["w_off"] + margin)) + # w_on sweep (low values): drop w_off if necessary so w_off < w_on. + if "w_on" in adjusted and adjusted["w_on"] - margin <= 0.20: + adjusted.setdefault("w_off", max(0.02, adjusted["w_on"] - margin)) + # w_opt sweep (low values): drop w_wilt if necessary so w_wilt < w_opt. + if "w_opt" in adjusted and adjusted["w_opt"] - margin <= 0.15: + adjusted.setdefault("w_wilt", max(0.02, adjusted["w_opt"] - margin)) + return adjusted + + +def _run_sample_scenario( + laitype: int, + param_overrides: Optional[Dict[str, float]] = None, + dry_start: bool = False, +) -> pd.Series: + from supy import SUEWSSimulation + + sim = SUEWSSimulation.from_sample_data() + state = sim.state_init.copy() + _set_laitype(state, laitype) + if dry_start: + _apply_dry_start(state) + if param_overrides: + for p, v in _co_adjust(param_overrides).items(): + _set_param(state, p, v) + + scenario = SUEWSSimulation.from_state(state) + scenario.update_forcing(sim.forcing) + output = scenario.run() + return output.SUEWS["LAI"].copy() + + +def _compute_metrics(lai: pd.Series, baseline: pd.Series) -> Dict[str, float]: + diff = lai.values - baseline.values + rmse = float(np.sqrt(np.mean(diff**2))) + mean = float(lai.mean()) + amplitude = float(lai.max() - lai.min()) + + green_up = float("nan") + # Daily aggregation; detect first DOY where LAI exceeds trough + 0.5 * amplitude. + try: + daily = lai.groupby(lai.index.get_level_values("datetime").dayofyear).mean() + if len(daily) > 0 and amplitude > 1.0e-6: + threshold = float(daily.min()) + 0.5 * amplitude + above = daily >= threshold + if above.any(): + green_up = float(above.idxmax()) + except Exception: # noqa: BLE001 - robust to missing index layers + pass + + return { + "rmse_vs_baseline": rmse, + "mean_lai": mean, + "amplitude": amplitude, + "green_up_doy": green_up, + } + + +def sweep_parameter( + param: str, + values: Sequence[float], + baseline: pd.Series, + dry_start: bool, +) -> List[SweepResult]: + results: List[SweepResult] = [] + for value in values: + lai = _run_sample_scenario( + laitype=2, param_overrides={param: float(value)}, dry_start=dry_start + ) + metrics = _compute_metrics(lai, baseline) + results.append( + SweepResult( + param=param, + value=float(value), + mean_lai=metrics["mean_lai"], + rmse_vs_baseline=metrics["rmse_vs_baseline"], + amplitude=metrics["amplitude"], + green_up_doy=metrics["green_up_doy"], + ) + ) + return results + + +def _write_sweep_outputs( + out_dir: Path, param: str, results: Sequence[SweepResult], baseline_mean: float +) -> None: + out_dir.mkdir(parents=True, exist_ok=True) + payload = { + "parameter": param, + "baseline_mean_lai": baseline_mean, + "points": [ + { + "value": r.value, + "mean_lai": r.mean_lai, + "rmse_vs_baseline": r.rmse_vs_baseline, + "amplitude": r.amplitude, + "green_up_doy": r.green_up_doy, + } + for r in results + ], + } + (out_dir / f"sweep_{param}.json").write_text( + json.dumps(payload, indent=2), encoding="utf-8" + ) + + try: + import matplotlib + + matplotlib.use("Agg") + import matplotlib.pyplot as plt + + fig, ax_left = plt.subplots(figsize=(7, 4)) + x = [r.value for r in results] + ax_left.plot(x, [r.mean_lai for r in results], "o-", color="tab:blue", label="mean LAI") + ax_left.axhline(baseline_mean, color="k", lw=0.7, ls="--", label="LAIType=0 baseline") + ax_left.set_xlabel(param) + ax_left.set_ylabel("mean LAI [m2/m2]", color="tab:blue") + ax_left.tick_params(axis="y", labelcolor="tab:blue") + + ax_right = ax_left.twinx() + ax_right.plot(x, [r.rmse_vs_baseline for r in results], "s-", color="tab:orange", label="RMSE vs baseline") + ax_right.set_ylabel("RMSE vs LAIType=0 [m2/m2]", color="tab:orange") + ax_right.tick_params(axis="y", labelcolor="tab:orange") + + ax_left.set_title(f"GH-1292 moisture-aware sensitivity: {param}") + lines1, labels1 = ax_left.get_legend_handles_labels() + lines2, labels2 = ax_right.get_legend_handles_labels() + ax_left.legend(lines1 + lines2, labels1 + labels2, loc="best", fontsize=8) + fig.tight_layout() + fig.savefig(out_dir / f"sweep_{param}.png", dpi=150) + plt.close(fig) + except ImportError: + print("matplotlib not available; skipping plot", file=sys.stderr) + + +def main() -> int: + args = parse_args() + out_dir = Path(".context/gh1292") / args.site + + baseline_series = _run_sample_scenario(laitype=0, dry_start=args.dry_start) + baseline_mean = float(baseline_series.mean()) + print(f"[{args.site}] LAIType=0 baseline mean LAI = {baseline_mean:.4f}") + + params_to_sweep: Iterable[str] + if args.all: + params_to_sweep = PARAM_NAMES + else: + params_to_sweep = [args.param] + + for param in params_to_sweep: + values = args.values if args.values and not args.all else DEFAULT_SWEEP_VALUES[param] + print(f" Sweeping {param} over {list(values)}") + results = sweep_parameter(param, values, baseline_series, dry_start=args.dry_start) + _write_sweep_outputs(out_dir, param, results, baseline_mean) + for r in results: + green_str = f"{r.green_up_doy:.0f}" if r.green_up_doy == r.green_up_doy else "nan" + print( + f" {param}={r.value:<6g} meanLAI={r.mean_lai:.4f} " + f"RMSE={r.rmse_vs_baseline:.4f} ampl={r.amplitude:.3f} greenDOY={green_str}" + ) + + print(f"Outputs under: {out_dir}") + return 0 + + +if __name__ == "__main__": + sys.exit(main()) diff --git a/src/suews/src/suews_phys_dailystate.f95 b/src/suews/src/suews_phys_dailystate.f95 index cf46300ca..76c8b0911 100644 --- a/src/suews/src/suews_phys_dailystate.f95 +++ b/src/suews/src/suews_phys_dailystate.f95 @@ -168,6 +168,10 @@ SUBROUTINE SUEWS_cal_DailyState( & LAI_id => phenState%LAI_id, & GDD_id => phenState%GDD_id, & SDD_id => phenState%SDD_id, & + ! --- GH-1292: moisture-aware phenology state (laitype=2) --- + wbar_id => phenState%wbar_id, & + w_id_prev => phenState%w_id_prev, & + leaf_on_permitted => phenState%leaf_on_permitted, & albDecTr_id => phenState%albDecTr_id, & albEveTr_id => phenState%albEveTr_id, & albGrass_id => phenState%albGrass_id, & @@ -336,8 +340,18 @@ SUBROUTINE SUEWS_cal_DailyState( & BaseT, BaseTe, & GDDFull, SDDFull, & LAIMin, LAIMax, LAIPower, LAIType, & + ! --- GH-1292: SMD inputs (mm); caller computes smd = SoilStoreCap - soilstore_surf for veg (surfaces 3-5) --- + SoilStoreCap(3:5) - soilstore_surf(3:5), SoilStoreCap(3:5), & + ! --- GH-1292: moisture-aware parameters (per-veg-surface arrays) --- + [evetrPrm%lai%w_wilt, dectrPrm%lai%w_wilt, grassPrm%lai%w_wilt], & + [evetrPrm%lai%w_opt, dectrPrm%lai%w_opt, grassPrm%lai%w_opt], & + [evetrPrm%lai%f_shape, dectrPrm%lai%f_shape, grassPrm%lai%f_shape], & + [evetrPrm%lai%w_on, dectrPrm%lai%w_on, grassPrm%lai%w_on], & + [evetrPrm%lai%w_off, dectrPrm%lai%w_off, grassPrm%lai%w_off], & + [evetrPrm%lai%tau_w, dectrPrm%lai%tau_w, grassPrm%lai%tau_w], & LAI_id_prev, & GDD_id, SDD_id, & !inout + wbar_id, w_id_prev, leaf_on_permitted, & !inout LAI_id) !output IF (supy_error_flag) RETURN @@ -532,8 +546,12 @@ SUBROUTINE update_GDDLAI( & BaseT_GDD, BaseT_SDD, & GDDFull, SDDFull, & LAIMin, LAIMax, LAIPower, LAIType, & + ! --- GH-1292 PR1: moisture-aware phenology scaffolding (no-op for laitype=2) --- + smd_veg_prev, smdcap_veg, & !input [mm] + w_wilt, w_opt, f_shape, w_on, w_off, tau_w, & !input (LAI_PRM moisture params; tau_w in days, rest dimensionless) LAI_id_prev, & GDD_id, SDD_id, & !inout + wbar_id, w_id_prev, leaf_on_permitted, & !inout (PHENOLOGY_STATE moisture-aware slots) LAI_id_next) !output IMPLICIT NONE @@ -561,10 +579,28 @@ SUBROUTINE update_GDDLAI( & REAL(KIND(1D0)), DIMENSION(nvegsurf), INTENT(IN) :: LAIMax !Max LAI [m2 m-2] REAL(KIND(1D0)), DIMENSION(4, nvegsurf), INTENT(IN) :: LAIPower !Coeffs for LAI equation: 1,2 - leaf growth; 3,4 - leaf off !! N.B. currently DecTr only, although input provided for all veg types - INTEGER, DIMENSION(nvegsurf), INTENT(IN) :: LAIType !LAI equation to use: original (0) or new (1) + INTEGER, DIMENSION(nvegsurf), INTENT(IN) :: LAIType !LAI equation to use: 0 original, 1 high-latitude, 2 moisture-aware (PR1 no-op) + + ! --- GH-1292 moisture-aware inputs (SMD framing, mm) ------------------------------------------- + REAL(KIND(1D0)), DIMENSION(nvegsurf), INTENT(IN) :: smd_veg_prev !Previous-day vegetation soil moisture deficit [mm] + REAL(KIND(1D0)), DIMENSION(nvegsurf), INTENT(IN) :: smdcap_veg !SMD capacity for vegetation surfaces [mm] + + ! --- GH-1292 moisture-aware parameters (thresholds on dimensionless w = 1 - smd/smdcap) -------- + REAL(KIND(1D0)), DIMENSION(nvegsurf), INTENT(IN) :: w_wilt !Jarvis lower bound on relative soil water [-] + REAL(KIND(1D0)), DIMENSION(nvegsurf), INTENT(IN) :: w_opt !Jarvis upper bound on relative soil water [-] + REAL(KIND(1D0)), DIMENSION(nvegsurf), INTENT(IN) :: f_shape !Jarvis shape exponent [-] + REAL(KIND(1D0)), DIMENSION(nvegsurf), INTENT(IN) :: w_on !Leaf-on persistence threshold [-] + REAL(KIND(1D0)), DIMENSION(nvegsurf), INTENT(IN) :: w_off !Leaf-off persistence threshold [-] + REAL(KIND(1D0)), DIMENSION(nvegsurf), INTENT(IN) :: tau_w !Persistence window for running-mean relative soil water [days] REAL(KIND(1D0)), DIMENSION(3), INTENT(INOUT) :: GDD_id !Growing Degree Days (see SUEWS_DailyState.f95) REAL(KIND(1D0)), DIMENSION(3), INTENT(INOUT) :: SDD_id !Senescence Degree Days (see SUEWS_DailyState.f95) + + ! --- GH-1292 moisture-aware state (relative soil water, dimensionless; logical latch) --------- + REAL(KIND(1D0)), DIMENSION(nvegsurf), INTENT(INOUT) :: wbar_id !Running-mean relative soil water [-] + REAL(KIND(1D0)), DIMENSION(nvegsurf), INTENT(INOUT) :: w_id_prev !Previous-day relative soil water snapshot [-] + LOGICAL, DIMENSION(nvegsurf), INTENT(INOUT) :: leaf_on_permitted !CLM5-style leaf-on latch + REAL(KIND(1D0)), DIMENSION(nvegsurf), INTENT(OUT) :: LAI_id_next !LAI for each veg surface [m2 m-2] REAL(KIND(1D0)), DIMENSION(nvegsurf), INTENT(IN) :: LAI_id_prev ! LAI of previous day @@ -574,6 +610,12 @@ SUBROUTINE update_GDDLAI( & REAL(KIND(1D0)), DIMENSION(3) :: GDD_id_prev ! GDD of previous day REAL(KIND(1D0)), DIMENSION(3) :: SDD_id_prev ! SDD of previous day + ! --- GH-1292 local moisture-aware working variables ------------------------------------------ + REAL(KIND(1D0)) :: w !Current-step relative soil water [-] + REAL(KIND(1D0)) :: f_w !Jarvis water-stress multiplier on delta_GDD [-] + LOGICAL :: persistence_ok !Whether wbar_id crosses the leaf-on threshold + LOGICAL :: thermal_ok !Whether delta_GDD > 0 (Tbar >= BaseT_GDD); companion to persistence_ok + INTEGER :: critDays INTEGER :: iv @@ -610,6 +652,61 @@ SUBROUTINE update_GDDLAI( & IF (delta_SDD > 0) delta_SDD = 0 !SDD cannot be positive + ! --- GH-1292 PR2: moisture-aware phenology (laitype=2, Design C) -------------------------- + ! Inputs are SMD (mm); internal thresholds operate on relative soil water w = 1 - smd/smdcap. + ! Hybrid scheme: Jarvis-style continuous stress factor on delta_GDD plus a CLM5-style + ! persistence latch (Lawrence et al. 2019 §14) that freezes thermal accumulation when the + ! tau_w-day running mean of relative soil water sits below the leaf-off threshold. + IF (LAIType(iv) == 2) THEN + ! Convert previous-day SMD (mm) to dimensionless relative soil water. + w = 1.0D0 - smd_veg_prev(iv)/MAX(smdcap_veg(iv), 1.0D-6) + IF (w < 0.0D0) w = 0.0D0 + IF (w > 1.0D0) w = 1.0D0 + + ! Seed the running mean on first activation so well-watered sites do not spend + ! tau_w days ramping up from an unset moisture history. Negative values are reserved + ! for this sentinel; w itself is clamped to the physical interval [0, 1]. + IF (w_id_prev(iv) < 0.0D0 .OR. wbar_id(iv) < 0.0D0) THEN + wbar_id(iv) = w + ELSE + wbar_id(iv) = wbar_id(iv) + (w - wbar_id(iv))/MAX(tau_w(iv), 1.0D0) + END IF + + ! Jarvis piecewise-power water-stress factor on delta_GDD. + IF (w <= w_wilt(iv)) THEN + f_w = 0.0D0 + ELSE IF (w >= w_opt(iv)) THEN + f_w = 1.0D0 + ELSE + f_w = ((w - w_wilt(iv))/MAX(w_opt(iv) - w_wilt(iv), 1.0D-6))**f_shape(iv) + END IF + + ! CLM5-style persistence latch: flips ON when the tau_w-day running mean crosses the + ! leaf-on threshold AND the thermal gate (delta_GDD > 0, i.e. Tbar >= BaseT_GDD) is + ! satisfied; flips OFF when the running mean drops below the leaf-off threshold. + thermal_ok = (delta_GDD > 0.0D0) + persistence_ok = (wbar_id(iv) >= w_on(iv)) + IF (leaf_on_permitted(iv)) THEN + IF (wbar_id(iv) < w_off(iv)) THEN + leaf_on_permitted(iv) = .FALSE. + END IF + ELSE + IF (persistence_ok .AND. thermal_ok) THEN + leaf_on_permitted(iv) = .TRUE. + END IF + END IF + + ! Apply the moisture gate. When the latch is open, delta_GDD is modulated continuously + ! by f_w; when the latch is closed (drought regime) thermal accumulation halts entirely. + IF (leaf_on_permitted(iv)) THEN + delta_GDD = delta_GDD*f_w + ELSE + delta_GDD = 0.0D0 + END IF + + w_id_prev(iv) = w + END IF + ! Calculate cumulative growing and senescence degree days GDD_id(iv) = GDD_id_prev(iv) + delta_GDD SDD_id(iv) = SDD_id_prev(iv) + delta_SDD @@ -642,7 +739,9 @@ SUBROUTINE update_GDDLAI( & ! Set GDD zero in winter time IF (SDD_id(iv) < -critDays .AND. id > 170) GDD_id(iv) = 0 - IF (LAItype(iv) < 0.5) THEN !Original LAI type + ! GH-1292: LAItype=2 shares the original power-law (0). In PR1 this makes laitype=2 a no-op + ! equivalent to laitype=0; PR2 layers the SMD gate on top of the same curve. + IF (LAItype(iv) == 0 .OR. LAItype(iv) == 2 .OR. LAItype(iv) < 0) THEN !Original LAI type (and GH-1292 scaffold) IF (GDD_id(iv) > 0 .AND. GDD_id(iv) < GDDFull(iv)) THEN !Leaves can still grow LAI_id_next(iv) = (LAI_id_prev(iv)**LAIPower(1, iv)*GDD_id(iv)*LAIPower(2, iv)) + LAI_id_prev(iv) ELSEIF (SDD_id(iv) < 0 .AND. SDD_id(iv) > SDDFull(iv)) THEN !Start senescence @@ -650,7 +749,10 @@ SUBROUTINE update_GDDLAI( & ELSE LAI_id_next(iv) = LAI_id_prev(iv) END IF - ELSEIF (LAItype(iv) >= 0.5) THEN + ELSE + ! Preserve the pre-GH-1292 behaviour for unsupported positive integers: + ! any non-zero value other than LAItype=2 falls back to the legacy + ! high-latitude branch instead of leaving LAI_id_next unset. IF (GDD_id(iv) > 0 .AND. GDD_id(iv) < GDDFull(iv)) THEN !Leaves can still grow LAI_id_next(iv) = (LAI_id_prev(iv)**LAIPower(1, iv)*GDD_id(iv)*LAIPower(2, iv)) + LAI_id_prev(iv) !! Use day length to start senescence at high latitudes (N hemisphere) @@ -669,7 +771,8 @@ SUBROUTINE update_GDDLAI( & ! Set GDD zero in winter time IF (SDD_id(iv) < -critDays .AND. id < 250) GDD_id(iv) = 0 - IF (LAItype(iv) < 0.5) THEN !Original LAI type + ! GH-1292: LAItype=2 shares the original southern-hemisphere power-law (0) in PR1. + IF (LAItype(iv) == 0 .OR. LAItype(iv) == 2 .OR. LAItype(iv) < 0) THEN !Original LAI type (and GH-1292 scaffold) IF (GDD_id(iv) > 0 .AND. GDD_id(iv) < GDDFull(iv)) THEN LAI_id_next(iv) = (LAI_id_prev(iv)**LAIPower(1, iv)*GDD_id(iv)*LAIPower(2, iv)) + LAI_id_prev(iv) ELSEIF (SDD_id(iv) < 0 .AND. SDD_id(iv) > SDDFull(iv)) THEN @@ -678,6 +781,9 @@ SUBROUTINE update_GDDLAI( & LAI_id_next(iv) = LAI_id_prev(iv) END IF ELSE + ! Preserve the pre-GH-1292 behaviour for unsupported positive integers: + ! any non-zero value other than LAItype=2 falls back to the legacy + ! high-latitude branch instead of leaving LAI_id_next unset. IF (GDD_id(iv) > 0 .AND. GDD_id(iv) < GDDFull(iv)) THEN LAI_id_next(iv) = (LAI_id_prev(iv)**LAIPower(1, iv)*GDD_id(iv)*LAIPower(2, iv)) + LAI_id_prev(iv) !! Day length not used to start senescence in S hemisphere (not much land) diff --git a/src/suews/src/suews_type_vegetation.f95 b/src/suews/src/suews_type_vegetation.f95 index b6b8bfb89..c422ea97d 100644 --- a/src/suews/src/suews_type_vegetation.f95 +++ b/src/suews/src/suews_type_vegetation.f95 @@ -24,7 +24,15 @@ module module_type_vegetation REAL(KIND(1D0)) :: laimin = 0.0D0 ! Minimum LAI [m2 m-2] REAL(KIND(1D0)) :: laimax = 0.0D0 ! Maximum LAI [m2 m-2] REAL(KIND(1D0)), DIMENSION(4) :: laipower = 0.0D0 ! Coefficients for LAI equation: 1,2 - leaf growth; 3,4 - leaf off [-] - INTEGER :: laitype = 0 ! LAI equation to use: original (0) or new (1) [-] + INTEGER :: laitype = 0 ! LAI equation to use: 0 original, 1 high-latitude, 2 moisture-aware (GH-1292; PR1 no-op) [-] + ! --- Moisture-aware phenology parameters (laitype=2) ------------------------------------------- + ! Thresholds operate on dimensionless relative soil water w = 1 - smd/smdcap (not on SMD in mm). + REAL(KIND(1D0)) :: w_wilt = 0.15D0 ! Lower bound of Jarvis water-stress factor; relative soil water [-] + REAL(KIND(1D0)) :: w_opt = 0.40D0 ! Upper bound of Jarvis water-stress factor; relative soil water [-] + REAL(KIND(1D0)) :: f_shape = 1.0D0 ! Shape exponent on the Jarvis water-stress factor [-] + REAL(KIND(1D0)) :: w_on = 0.35D0 ! Leaf-on persistence threshold (CLM5-style); relative soil water [-] + REAL(KIND(1D0)) :: w_off = 0.20D0 ! Leaf-off persistence threshold; relative soil water [-] + REAL(KIND(1D0)) :: tau_w = 15.0D0 ! Persistence window for running-mean relative soil water [days] END TYPE LAI_PRM TYPE, PUBLIC :: PHENOLOGY_STATE @@ -50,6 +58,11 @@ module module_type_vegetation REAL(KIND(1D0)) :: g_ta = 0.0D0 ! Surface conductance function for air temperature [-] REAL(KIND(1D0)) :: g_smd = 0.0D0 ! Surface conductance function for soil moisture deficit [-] REAL(KIND(1D0)) :: g_lai = 0.0D0 ! Surface conductance function for LAI [-] + ! --- Moisture-aware phenology state (laitype=2) ------------------------------------------- + ! Carried as dimensionless relative soil water w = 1 - smd/smdcap; SMD->w conversion lives in update_GDDLAI. + REAL(KIND(1D0)), DIMENSION(nvegsurf) :: wbar_id = -1.0D0 ! Running-mean relative soil water for phenology [-]; negative marks an unset moisture history + REAL(KIND(1D0)), DIMENSION(nvegsurf) :: w_id_prev = -1.0D0 ! Previous-day relative soil water snapshot [-]; negative marks an unset moisture history + LOGICAL, DIMENSION(nvegsurf) :: leaf_on_permitted = .TRUE. ! CLM5-style leaf-on latch; default permitted so well-watered sites do not wait for the trigger ! flag for iteration safety - NO ! GDD_id, SDD_id are extensive quantities and thus cannot be used for iteration safety LOGICAL :: iter_safe = .FALSE. diff --git a/src/suews_bridge/c_api/driver.f95 b/src/suews_bridge/c_api/driver.f95 index 48a372509..8409f71ef 100644 --- a/src/suews_bridge/c_api/driver.f95 +++ b/src/suews_bridge/c_api/driver.f95 @@ -101,7 +101,7 @@ module module_c_api_driver integer(c_int), parameter :: SUEWS_CAPI_OHM_STATE_LEN = 53_c_int integer(c_int), parameter :: SUEWS_CAPI_SOLAR_STATE_LEN = 3_c_int integer(c_int), parameter :: SUEWS_CAPI_ATM_STATE_LEN = 39_c_int -integer(c_int), parameter :: SUEWS_CAPI_PHENOLOGY_STATE_LEN = 76_c_int +integer(c_int), parameter :: SUEWS_CAPI_PHENOLOGY_STATE_LEN = 85_c_int integer(c_int), parameter :: SUEWS_CAPI_SNOW_STATE_LEN = 79_c_int integer(c_int), parameter :: SUEWS_CAPI_HYDRO_STATE_BASE_LEN = 10_c_int * int(nsurf, c_int) + 34_c_int integer(c_int), parameter :: SUEWS_CAPI_HEAT_STATE_BASE_LEN = 7_c_int * int(nsurf, c_int) + 79_c_int @@ -1594,6 +1594,23 @@ subroutine pack_phenology_state(s, flat, n_flat, err) flat(idx) = s%g_ta; idx = idx + 1_c_int flat(idx) = s%g_smd; idx = idx + 1_c_int flat(idx) = s%g_lai; idx = idx + 1_c_int + + ! GH-1292 moisture-aware state: wbar_id, w_id_prev, leaf_on_permitted (per veg surface). + do i = 1_c_int, int(nvegsurf, c_int) + flat(idx) = s%wbar_id(i) + idx = idx + 1_c_int + end do + + do i = 1_c_int, int(nvegsurf, c_int) + flat(idx) = s%w_id_prev(i) + idx = idx + 1_c_int + end do + + do i = 1_c_int, int(nvegsurf, c_int) + flat(idx) = merge(1.0_c_double, 0.0_c_double, s%leaf_on_permitted(i)) + idx = idx + 1_c_int + end do + flat(idx) = merge(1.0_c_double, 0.0_c_double, s%iter_safe) err = SUEWS_CAPI_OK diff --git a/src/suews_bridge/c_api/lai.f95 b/src/suews_bridge/c_api/lai.f95 index 04bd3e061..f4ff58b67 100644 --- a/src/suews_bridge/c_api/lai.f95 +++ b/src/suews_bridge/c_api/lai.f95 @@ -15,8 +15,8 @@ MODULE module_c_api_lai PUBLIC :: SUEWS_CAPI_BAD_BUFFER PUBLIC :: SUEWS_CAPI_BAD_STATE - INTEGER(c_int), PARAMETER, PUBLIC :: SUEWS_CAPI_LAI_PRM_LEN = 11_c_int - INTEGER(c_int), PARAMETER, PUBLIC :: SUEWS_CAPI_LAI_PRM_SCHEMA_VERSION = 1_c_int + INTEGER(c_int), PARAMETER, PUBLIC :: SUEWS_CAPI_LAI_PRM_LEN = 17_c_int + INTEGER(c_int), PARAMETER, PUBLIC :: SUEWS_CAPI_LAI_PRM_SCHEMA_VERSION = 2_c_int TYPE :: lai_prm_shadow REAL(c_double) :: baset = 0.0_c_double @@ -27,6 +27,13 @@ MODULE module_c_api_lai REAL(c_double) :: laimax = 0.0_c_double REAL(c_double), DIMENSION(4) :: laipower = 0.0_c_double INTEGER(c_int) :: laitype = 0_c_int + ! GH-1292 moisture-aware phenology parameters (laitype=2); relative soil water thresholds + days. + REAL(c_double) :: w_wilt = 0.15_c_double + REAL(c_double) :: w_opt = 0.40_c_double + REAL(c_double) :: f_shape = 1.0_c_double + REAL(c_double) :: w_on = 0.35_c_double + REAL(c_double) :: w_off = 0.20_c_double + REAL(c_double) :: tau_w = 15.0_c_double END TYPE lai_prm_shadow PUBLIC :: suews_lai_prm_len @@ -98,6 +105,12 @@ SUBROUTINE lai_prm_pack(state, flat, n_flat, err) flat(9) = state%laipower(3) flat(10) = state%laipower(4) flat(11) = REAL(state%laitype, c_double) + flat(12) = state%w_wilt + flat(13) = state%w_opt + flat(14) = state%f_shape + flat(15) = state%w_on + flat(16) = state%w_off + flat(17) = state%tau_w err = SUEWS_CAPI_OK diff --git a/src/suews_bridge/c_api/lc_dectr_prm.f95 b/src/suews_bridge/c_api/lc_dectr_prm.f95 index 7ac011aa8..8229e0bdc 100644 --- a/src/suews_bridge/c_api/lc_dectr_prm.f95 +++ b/src/suews_bridge/c_api/lc_dectr_prm.f95 @@ -16,8 +16,8 @@ module module_c_api_lc_dectr_prm public :: SUEWS_CAPI_BAD_BUFFER public :: SUEWS_CAPI_BAD_STATE -integer(c_int), parameter, public :: SUEWS_CAPI_LC_DECTR_PRM_LEN = 61_c_int -integer(c_int), parameter, public :: SUEWS_CAPI_LC_DECTR_PRM_SCHEMA_VERSION = 1_c_int +integer(c_int), parameter, public :: SUEWS_CAPI_LC_DECTR_PRM_LEN = 67_c_int +integer(c_int), parameter, public :: SUEWS_CAPI_LC_DECTR_PRM_SCHEMA_VERSION = 2_c_int type :: ohm_coef_lc_shadow real(c_double) :: summer_dry = 0.0_c_double @@ -61,6 +61,13 @@ module module_c_api_lc_dectr_prm real(c_double) :: laimax = 0.0_c_double real(c_double), dimension(4) :: laipower = 0.0_c_double integer(c_int) :: laitype = 0_c_int + ! GH-1292 moisture-aware phenology parameters (laitype=2). + real(c_double) :: w_wilt = 0.15_c_double + real(c_double) :: w_opt = 0.40_c_double + real(c_double) :: f_shape = 1.0_c_double + real(c_double) :: w_on = 0.35_c_double + real(c_double) :: w_off = 0.20_c_double + real(c_double) :: tau_w = 15.0_c_double end type lai_prm_shadow type :: water_dist_prm_shadow @@ -208,6 +215,13 @@ subroutine lc_dectr_prm_pack(state, flat, n_flat, err) flat(idx) = state%lai%laipower(i); idx = idx + 1 end do flat(idx) = real(state%lai%laitype, c_double); idx = idx + 1 + ! GH-1292 moisture-aware phenology parameters (laitype=2). + flat(idx) = state%lai%w_wilt; idx = idx + 1 + flat(idx) = state%lai%w_opt; idx = idx + 1 + flat(idx) = state%lai%f_shape; idx = idx + 1 + flat(idx) = state%lai%w_on; idx = idx + 1 + flat(idx) = state%lai%w_off; idx = idx + 1 + flat(idx) = state%lai%tau_w; idx = idx + 1 flat(idx) = state%waterdist%to_paved; idx = idx + 1 flat(idx) = state%waterdist%to_bldg; idx = idx + 1 @@ -292,6 +306,13 @@ subroutine lc_dectr_prm_unpack(flat, n_flat, state, err) idx = idx + 1_c_int end do state%lai%laitype = int(nint(flat(idx))); idx = idx + 1_c_int + ! GH-1292 moisture-aware phenology parameters (laitype=2). + state%lai%w_wilt = flat(idx); idx = idx + 1_c_int + state%lai%w_opt = flat(idx); idx = idx + 1_c_int + state%lai%f_shape = flat(idx); idx = idx + 1_c_int + state%lai%w_on = flat(idx); idx = idx + 1_c_int + state%lai%w_off = flat(idx); idx = idx + 1_c_int + state%lai%tau_w = flat(idx); idx = idx + 1_c_int state%waterdist%to_paved = flat(idx); idx = idx + 1_c_int state%waterdist%to_bldg = flat(idx); idx = idx + 1_c_int diff --git a/src/suews_bridge/c_api/lc_evetr_prm.f95 b/src/suews_bridge/c_api/lc_evetr_prm.f95 index e408df1f8..c1b9c6f1e 100644 --- a/src/suews_bridge/c_api/lc_evetr_prm.f95 +++ b/src/suews_bridge/c_api/lc_evetr_prm.f95 @@ -16,8 +16,8 @@ module module_c_api_lc_evetr_prm public :: SUEWS_CAPI_BAD_BUFFER public :: SUEWS_CAPI_BAD_STATE -integer(c_int), parameter, public :: SUEWS_CAPI_LC_EVETR_PRM_LEN = 57_c_int -integer(c_int), parameter, public :: SUEWS_CAPI_LC_EVETR_PRM_SCHEMA_VERSION = 1_c_int +integer(c_int), parameter, public :: SUEWS_CAPI_LC_EVETR_PRM_LEN = 63_c_int +integer(c_int), parameter, public :: SUEWS_CAPI_LC_EVETR_PRM_SCHEMA_VERSION = 2_c_int type :: ohm_coef_lc_shadow real(c_double) :: summer_dry = 0.0_c_double @@ -61,6 +61,13 @@ module module_c_api_lc_evetr_prm real(c_double) :: laimax = 0.0_c_double real(c_double), dimension(4) :: laipower = 0.0_c_double integer(c_int) :: laitype = 0_c_int + ! GH-1292 moisture-aware phenology parameters (laitype=2). + real(c_double) :: w_wilt = 0.15_c_double + real(c_double) :: w_opt = 0.40_c_double + real(c_double) :: f_shape = 1.0_c_double + real(c_double) :: w_on = 0.35_c_double + real(c_double) :: w_off = 0.20_c_double + real(c_double) :: tau_w = 15.0_c_double end type lai_prm_shadow type :: water_dist_prm_shadow @@ -200,6 +207,13 @@ subroutine lc_evetr_prm_pack(state, flat, n_flat, err) flat(idx) = state%lai%laipower(i); idx = idx + 1 end do flat(idx) = real(state%lai%laitype, c_double); idx = idx + 1 + ! GH-1292 moisture-aware phenology parameters (laitype=2). + flat(idx) = state%lai%w_wilt; idx = idx + 1 + flat(idx) = state%lai%w_opt; idx = idx + 1 + flat(idx) = state%lai%f_shape; idx = idx + 1 + flat(idx) = state%lai%w_on; idx = idx + 1 + flat(idx) = state%lai%w_off; idx = idx + 1 + flat(idx) = state%lai%tau_w; idx = idx + 1 flat(idx) = state%waterdist%to_paved; idx = idx + 1 flat(idx) = state%waterdist%to_bldg; idx = idx + 1 @@ -280,6 +294,13 @@ subroutine lc_evetr_prm_unpack(flat, n_flat, state, err) idx = idx + 1_c_int end do state%lai%laitype = int(nint(flat(idx))); idx = idx + 1_c_int + ! GH-1292 moisture-aware phenology parameters (laitype=2). + state%lai%w_wilt = flat(idx); idx = idx + 1_c_int + state%lai%w_opt = flat(idx); idx = idx + 1_c_int + state%lai%f_shape = flat(idx); idx = idx + 1_c_int + state%lai%w_on = flat(idx); idx = idx + 1_c_int + state%lai%w_off = flat(idx); idx = idx + 1_c_int + state%lai%tau_w = flat(idx); idx = idx + 1_c_int state%waterdist%to_paved = flat(idx); idx = idx + 1_c_int state%waterdist%to_bldg = flat(idx); idx = idx + 1_c_int diff --git a/src/suews_bridge/c_api/lc_grass_prm.f95 b/src/suews_bridge/c_api/lc_grass_prm.f95 index a61d3e558..ffb7ae942 100644 --- a/src/suews_bridge/c_api/lc_grass_prm.f95 +++ b/src/suews_bridge/c_api/lc_grass_prm.f95 @@ -16,8 +16,8 @@ module module_c_api_lc_grass_prm public :: SUEWS_CAPI_BAD_BUFFER public :: SUEWS_CAPI_BAD_STATE -integer(c_int), parameter, public :: SUEWS_CAPI_LC_GRASS_PRM_LEN = 55_c_int -integer(c_int), parameter, public :: SUEWS_CAPI_LC_GRASS_PRM_SCHEMA_VERSION = 1_c_int +integer(c_int), parameter, public :: SUEWS_CAPI_LC_GRASS_PRM_LEN = 61_c_int +integer(c_int), parameter, public :: SUEWS_CAPI_LC_GRASS_PRM_SCHEMA_VERSION = 2_c_int type :: ohm_coef_lc_shadow real(c_double) :: summer_dry = 0.0_c_double @@ -61,6 +61,13 @@ module module_c_api_lc_grass_prm real(c_double) :: laimax = 0.0_c_double real(c_double), dimension(4) :: laipower = 0.0_c_double integer(c_int) :: laitype = 0_c_int + ! GH-1292 moisture-aware phenology parameters (laitype=2). + real(c_double) :: w_wilt = 0.15_c_double + real(c_double) :: w_opt = 0.40_c_double + real(c_double) :: f_shape = 1.0_c_double + real(c_double) :: w_on = 0.35_c_double + real(c_double) :: w_off = 0.20_c_double + real(c_double) :: tau_w = 15.0_c_double end type lai_prm_shadow type :: water_dist_prm_shadow @@ -196,6 +203,13 @@ subroutine lc_grass_prm_pack(state, flat, n_flat, err) flat(idx) = state%lai%laipower(i); idx = idx + 1 end do flat(idx) = real(state%lai%laitype, c_double); idx = idx + 1 + ! GH-1292 moisture-aware phenology parameters (laitype=2). + flat(idx) = state%lai%w_wilt; idx = idx + 1 + flat(idx) = state%lai%w_opt; idx = idx + 1 + flat(idx) = state%lai%f_shape; idx = idx + 1 + flat(idx) = state%lai%w_on; idx = idx + 1 + flat(idx) = state%lai%w_off; idx = idx + 1 + flat(idx) = state%lai%tau_w; idx = idx + 1 flat(idx) = state%waterdist%to_paved; idx = idx + 1 flat(idx) = state%waterdist%to_bldg; idx = idx + 1 @@ -274,6 +288,13 @@ subroutine lc_grass_prm_unpack(flat, n_flat, state, err) idx = idx + 1_c_int end do state%lai%laitype = int(nint(flat(idx))); idx = idx + 1_c_int + ! GH-1292 moisture-aware phenology parameters (laitype=2). + state%lai%w_wilt = flat(idx); idx = idx + 1_c_int + state%lai%w_opt = flat(idx); idx = idx + 1_c_int + state%lai%f_shape = flat(idx); idx = idx + 1_c_int + state%lai%w_on = flat(idx); idx = idx + 1_c_int + state%lai%w_off = flat(idx); idx = idx + 1_c_int + state%lai%tau_w = flat(idx); idx = idx + 1_c_int state%waterdist%to_paved = flat(idx); idx = idx + 1_c_int state%waterdist%to_bldg = flat(idx); idx = idx + 1_c_int diff --git a/src/suews_bridge/c_api/phenology.f95 b/src/suews_bridge/c_api/phenology.f95 index 7802202d3..476008a6d 100644 --- a/src/suews_bridge/c_api/phenology.f95 +++ b/src/suews_bridge/c_api/phenology.f95 @@ -17,8 +17,8 @@ module module_c_api_phenology public :: SUEWS_CAPI_BAD_BUFFER public :: SUEWS_CAPI_BAD_STATE -integer(c_int), parameter, public :: SUEWS_CAPI_PHENOLOGY_STATE_LEN = 76_c_int -integer(c_int), parameter, public :: SUEWS_CAPI_PHENOLOGY_STATE_SCHEMA_VERSION = 1_c_int +integer(c_int), parameter, public :: SUEWS_CAPI_PHENOLOGY_STATE_LEN = 85_c_int +integer(c_int), parameter, public :: SUEWS_CAPI_PHENOLOGY_STATE_SCHEMA_VERSION = 2_c_int type :: phenology_state_shadow real(c_double), dimension(nsurf) :: alb = 0.0_c_double @@ -43,6 +43,10 @@ module module_c_api_phenology real(c_double) :: g_ta = 0.0_c_double real(c_double) :: g_smd = 0.0_c_double real(c_double) :: g_lai = 0.0_c_double + ! GH-1292 moisture-aware phenology state (laitype=2). leaf_on_permitted encoded as 0.0/1.0 on the wire. + real(c_double), dimension(nvegsurf) :: wbar_id = -1.0_c_double + real(c_double), dimension(nvegsurf) :: w_id_prev = -1.0_c_double + logical, dimension(nvegsurf) :: leaf_on_permitted = .true. logical :: iter_safe = .false. end type phenology_state_shadow @@ -153,6 +157,22 @@ subroutine phenology_state_pack(state, flat, n_flat, err) flat(idx) = state%g_ta; idx = idx + 1 flat(idx) = state%g_smd; idx = idx + 1 flat(idx) = state%g_lai; idx = idx + 1 + + do i = 1, nvegsurf + flat(idx) = state%wbar_id(i) + idx = idx + 1 + end do + + do i = 1, nvegsurf + flat(idx) = state%w_id_prev(i) + idx = idx + 1 + end do + + do i = 1, nvegsurf + flat(idx) = merge(1.0_c_double, 0.0_c_double, state%leaf_on_permitted(i)) + idx = idx + 1 + end do + flat(idx) = merge(1.0_c_double, 0.0_c_double, state%iter_safe) err = SUEWS_CAPI_OK @@ -222,6 +242,22 @@ subroutine phenology_state_unpack(flat, n_flat, state, err) state%g_ta = flat(idx); idx = idx + 1_c_int state%g_smd = flat(idx); idx = idx + 1_c_int state%g_lai = flat(idx); idx = idx + 1_c_int + + do i = 1_c_int, int(nvegsurf, c_int) + state%wbar_id(i) = flat(idx) + idx = idx + 1_c_int + end do + + do i = 1_c_int, int(nvegsurf, c_int) + state%w_id_prev(i) = flat(idx) + idx = idx + 1_c_int + end do + + do i = 1_c_int, int(nvegsurf, c_int) + state%leaf_on_permitted(i) = flat(idx)>=0.5_c_double + idx = idx + 1_c_int + end do + state%iter_safe = flat(idx)>=0.5_c_double err = SUEWS_CAPI_OK diff --git a/src/suews_bridge/src/lai.rs b/src/suews_bridge/src/lai.rs index 58f7333e3..60abc95c2 100644 --- a/src/suews_bridge/src/lai.rs +++ b/src/suews_bridge/src/lai.rs @@ -2,8 +2,8 @@ use crate::codec::{validate_flat_len, StateCodec, TypeSchema, ValuesPayload}; use crate::error::BridgeError; use crate::ffi; -pub const LAI_PRM_FLAT_LEN: usize = 11; -pub const LAI_PRM_SCHEMA_VERSION: u32 = 1; +pub const LAI_PRM_FLAT_LEN: usize = 17; +pub const LAI_PRM_SCHEMA_VERSION: u32 = 2; pub type LaiPrmSchema = crate::codec::SimpleSchema; @@ -19,6 +19,14 @@ pub struct LaiPrm { pub laimax: f64, pub laipower: [f64; 4], pub laitype: i32, + // GH-1292 moisture-aware phenology parameters (laitype=2); thresholds operate on + // dimensionless relative soil water w = 1 - smd/smdcap, except tau_w which is days. + pub w_wilt: f64, + pub w_opt: f64, + pub f_shape: f64, + pub w_on: f64, + pub w_off: f64, + pub tau_w: f64, } impl Default for LaiPrm { @@ -32,6 +40,12 @@ impl Default for LaiPrm { laimax: 0.0, laipower: [0.0; 4], laitype: 0, + w_wilt: 0.15, + w_opt: 0.40, + f_shape: 1.0, + w_on: 0.35, + w_off: 0.20, + tau_w: 15.0, } } } @@ -65,6 +79,12 @@ impl LaiPrm { laimax: flat[5], laipower: [flat[6], flat[7], flat[8], flat[9]], laitype: decode_int(flat[10])?, + w_wilt: flat[11], + w_opt: flat[12], + f_shape: flat[13], + w_on: flat[14], + w_off: flat[15], + tau_w: flat[16], }) } @@ -81,6 +101,12 @@ impl LaiPrm { self.laipower[2], self.laipower[3], self.laitype as f64, + self.w_wilt, + self.w_opt, + self.f_shape, + self.w_on, + self.w_off, + self.tau_w, ] } } @@ -117,6 +143,12 @@ pub fn lai_prm_field_names() -> Vec { "laipower_3".to_string(), "laipower_4".to_string(), "laitype".to_string(), + "w_wilt".to_string(), + "w_opt".to_string(), + "f_shape".to_string(), + "w_on".to_string(), + "w_off".to_string(), + "tau_w".to_string(), ] } @@ -187,4 +219,29 @@ mod tests { let err = LaiPrm::from_flat(&flat).expect_err("fractional integer should fail"); assert_eq!(err, BridgeError::BadState); } + + #[test] + fn moisture_fields_roundtrip() { + // GH-1292 PR1: verify the six new LAI_PRM moisture slots serialise through to_flat/from_flat. + let original = LaiPrm { + w_wilt: 0.12, + w_opt: 0.48, + f_shape: 1.3, + w_on: 0.38, + w_off: 0.18, + tau_w: 20.0, + ..LaiPrm::default() + }; + let flat = original.to_flat(); + assert_eq!(flat.len(), LAI_PRM_FLAT_LEN); + let decoded = LaiPrm::from_flat(&flat).expect("round-trip should succeed"); + assert_eq!(original, decoded); + } + + #[test] + fn schema_version_is_two() { + // GH-1292 PR1: schema bumped because LAI_PRM layout changed. + assert_eq!(LAI_PRM_SCHEMA_VERSION, 2); + assert_eq!(LAI_PRM_FLAT_LEN, 17); + } } diff --git a/src/suews_bridge/src/lc_dectr_prm.rs b/src/suews_bridge/src/lc_dectr_prm.rs index 16467d954..ac6572e5e 100644 --- a/src/suews_bridge/src/lc_dectr_prm.rs +++ b/src/suews_bridge/src/lc_dectr_prm.rs @@ -7,8 +7,10 @@ use crate::ohm_prm::{ohm_prm_field_names, OhmPrm}; use crate::soil::{soil_prm_field_names, SoilPrm}; use crate::water_dist::{water_dist_prm_field_names, WaterDistPrm}; -pub const LC_DECTR_PRM_FLAT_LEN: usize = 61; -pub const LC_DECTR_PRM_SCHEMA_VERSION: u32 = 1; +// GH-1292 PR1: LAI_PRM grew by six moisture fields (11 -> 17), so the embedded LAI slice widens +// LC_DECTR_PRM by 6. Schema version bumped accordingly. +pub const LC_DECTR_PRM_FLAT_LEN: usize = 67; +pub const LC_DECTR_PRM_SCHEMA_VERSION: u32 = 2; pub type LcDectrPrmSchema = crate::codec::SimpleSchema; @@ -85,8 +87,8 @@ impl LcDectrPrm { wetthresh: flat[32], bioco2: BioCo2Prm::from_flat(&flat[33..41])?, maxconductance: flat[41], - lai: LaiPrm::from_flat(&flat[42..53])?, - waterdist: WaterDistPrm::from_flat(&flat[53..61])?, + lai: LaiPrm::from_flat(&flat[42..59])?, + waterdist: WaterDistPrm::from_flat(&flat[59..67])?, }) } diff --git a/src/suews_bridge/src/lc_evetr_prm.rs b/src/suews_bridge/src/lc_evetr_prm.rs index 0d727afe9..989fee621 100644 --- a/src/suews_bridge/src/lc_evetr_prm.rs +++ b/src/suews_bridge/src/lc_evetr_prm.rs @@ -7,8 +7,9 @@ use crate::ohm_prm::{ohm_prm_field_names, OhmPrm}; use crate::soil::{soil_prm_field_names, SoilPrm}; use crate::water_dist::{water_dist_prm_field_names, WaterDistPrm}; -pub const LC_EVETR_PRM_FLAT_LEN: usize = 57; -pub const LC_EVETR_PRM_SCHEMA_VERSION: u32 = 1; +// GH-1292 PR1: LAI_PRM grew by six moisture fields (11 -> 17). +pub const LC_EVETR_PRM_FLAT_LEN: usize = 63; +pub const LC_EVETR_PRM_SCHEMA_VERSION: u32 = 2; pub type LcEvetrPrmSchema = crate::codec::SimpleSchema; @@ -73,8 +74,8 @@ impl LcEvetrPrm { wetthresh: flat[28], bioco2: BioCo2Prm::from_flat(&flat[29..37])?, maxconductance: flat[37], - lai: LaiPrm::from_flat(&flat[38..49])?, - waterdist: WaterDistPrm::from_flat(&flat[49..57])?, + lai: LaiPrm::from_flat(&flat[38..55])?, + waterdist: WaterDistPrm::from_flat(&flat[55..63])?, }) } diff --git a/src/suews_bridge/src/lc_grass_prm.rs b/src/suews_bridge/src/lc_grass_prm.rs index 603f214c2..2799ebdb3 100644 --- a/src/suews_bridge/src/lc_grass_prm.rs +++ b/src/suews_bridge/src/lc_grass_prm.rs @@ -7,8 +7,9 @@ use crate::ohm_prm::{ohm_prm_field_names, OhmPrm}; use crate::soil::{soil_prm_field_names, SoilPrm}; use crate::water_dist::{water_dist_prm_field_names, WaterDistPrm}; -pub const LC_GRASS_PRM_FLAT_LEN: usize = 55; -pub const LC_GRASS_PRM_SCHEMA_VERSION: u32 = 1; +// GH-1292 PR1: LAI_PRM grew by six moisture fields (11 -> 17). +pub const LC_GRASS_PRM_FLAT_LEN: usize = 61; +pub const LC_GRASS_PRM_SCHEMA_VERSION: u32 = 2; pub type LcGrassPrmSchema = crate::codec::SimpleSchema; @@ -67,8 +68,8 @@ impl LcGrassPrm { wetthresh: flat[26], bioco2: BioCo2Prm::from_flat(&flat[27..35])?, maxconductance: flat[35], - lai: LaiPrm::from_flat(&flat[36..47])?, - waterdist: WaterDistPrm::from_flat(&flat[47..55])?, + lai: LaiPrm::from_flat(&flat[36..53])?, + waterdist: WaterDistPrm::from_flat(&flat[53..61])?, }) } diff --git a/src/suews_bridge/src/phenology.rs b/src/suews_bridge/src/phenology.rs index 65e50126f..4c13d941b 100644 --- a/src/suews_bridge/src/phenology.rs +++ b/src/suews_bridge/src/phenology.rs @@ -5,8 +5,9 @@ use crate::NSURF; pub const PHENOLOGY_STATE_NVEGSURF: usize = 3; pub const PHENOLOGY_STATE_STORE_DRAIN_ROWS: usize = 6; -pub const PHENOLOGY_STATE_FLAT_LEN: usize = 76; -pub const PHENOLOGY_STATE_SCHEMA_VERSION: u32 = 1; +// GH-1292 PR1: +3 per-veg-surface slots (wbar_id, w_id_prev, leaf_on_permitted) -> 76 + 9 = 85. +pub const PHENOLOGY_STATE_FLAT_LEN: usize = 85; +pub const PHENOLOGY_STATE_SCHEMA_VERSION: u32 = 2; pub type PhenologyStateSchema = crate::codec::SimpleSchema; @@ -36,6 +37,11 @@ pub struct PhenologyState { pub g_ta: f64, pub g_smd: f64, pub g_lai: f64, + // GH-1292 moisture-aware phenology state (laitype=2); relative soil water (dimensionless) + // and a per-surface logical latch. Wire format: each bool is encoded as f64 0.0/1.0. + pub wbar_id: [f64; PHENOLOGY_STATE_NVEGSURF], + pub w_id_prev: [f64; PHENOLOGY_STATE_NVEGSURF], + pub leaf_on_permitted: [bool; PHENOLOGY_STATE_NVEGSURF], pub iter_safe: bool, } @@ -64,6 +70,11 @@ impl Default for PhenologyState { g_ta: 0.0, g_smd: 0.0, g_lai: 0.0, + wbar_id: [-1.0; PHENOLOGY_STATE_NVEGSURF], + w_id_prev: [-1.0; PHENOLOGY_STATE_NVEGSURF], + // GH-1292 PR2: default latch to true so well-watered sites do not + // wait for the tau-day running-mean trigger to ramp up. + leaf_on_permitted: [true; PHENOLOGY_STATE_NVEGSURF], iter_safe: false, } } @@ -122,6 +133,19 @@ impl PhenologyState { state.g_ta = next(); state.g_smd = next(); state.g_lai = next(); + + for i in 0..PHENOLOGY_STATE_NVEGSURF { + state.wbar_id[i] = next(); + } + + for i in 0..PHENOLOGY_STATE_NVEGSURF { + state.w_id_prev[i] = next(); + } + + for i in 0..PHENOLOGY_STATE_NVEGSURF { + state.leaf_on_permitted[i] = next() >= 0.5; + } + state.iter_safe = next() >= 0.5; Ok(state) @@ -157,6 +181,13 @@ impl PhenologyState { flat.push(self.g_ta); flat.push(self.g_smd); flat.push(self.g_lai); + + flat.extend_from_slice(&self.wbar_id); + flat.extend_from_slice(&self.w_id_prev); + for permitted in &self.leaf_on_permitted { + flat.push(if *permitted { 1.0 } else { 0.0 }); + } + flat.push(if self.iter_safe { 1.0 } else { 0.0 }); flat @@ -233,6 +264,19 @@ pub fn phenology_state_field_names() -> Vec { names.push("g_ta".to_string()); names.push("g_smd".to_string()); names.push("g_lai".to_string()); + + for veg_surface in veg_surface_names() { + names.push(format!("wbar_id.{veg_surface}")); + } + + for veg_surface in veg_surface_names() { + names.push(format!("w_id_prev.{veg_surface}")); + } + + for veg_surface in veg_surface_names() { + names.push(format!("leaf_on_permitted.{veg_surface}")); + } + names.push("iter_safe".to_string()); names @@ -303,4 +347,26 @@ mod tests { .expect_err("payload with schema mismatch should fail"); assert_eq!(err, BridgeError::BadState); } + + #[test] + fn moisture_aware_state_roundtrip() { + // GH-1292 PR1: verify wbar_id, w_id_prev and leaf_on_permitted round-trip through to_flat/from_flat. + let original = PhenologyState { + wbar_id: [0.30, 0.45, 0.12], + w_id_prev: [0.33, 0.48, 0.09], + leaf_on_permitted: [true, false, true], + ..PhenologyState::default() + }; + let flat = original.to_flat(); + assert_eq!(flat.len(), PHENOLOGY_STATE_FLAT_LEN); + let decoded = PhenologyState::from_flat(&flat).expect("round-trip should succeed"); + assert_eq!(original, decoded); + } + + #[test] + fn schema_version_is_two() { + // GH-1292 PR1: schema bumped because PHENOLOGY_STATE gained three per-veg-surface slots. + assert_eq!(PHENOLOGY_STATE_SCHEMA_VERSION, 2); + assert_eq!(PHENOLOGY_STATE_FLAT_LEN, 85); + } } diff --git a/src/suews_bridge/src/yaml_config.rs b/src/suews_bridge/src/yaml_config.rs index c029e6921..7192c7427 100644 --- a/src/suews_bridge/src/yaml_config.rs +++ b/src/suews_bridge/src/yaml_config.rs @@ -492,6 +492,25 @@ fn apply_lai_overrides(lai: &mut crate::lai::LaiPrm, lc_root: &Value) { if let Some(v) = read_i32(lc_root, &["lai", "laitype"]) { lai.laitype = v; } + // GH-1292: moisture-aware phenology parameters (laitype=2). + if let Some(v) = read_numeric(lc_root, &["lai", "w_wilt"]) { + lai.w_wilt = v; + } + if let Some(v) = read_numeric(lc_root, &["lai", "w_opt"]) { + lai.w_opt = v; + } + if let Some(v) = read_numeric(lc_root, &["lai", "f_shape"]) { + lai.f_shape = v; + } + if let Some(v) = read_numeric(lc_root, &["lai", "w_on"]) { + lai.w_on = v; + } + if let Some(v) = read_numeric(lc_root, &["lai", "w_off"]) { + lai.w_off = v; + } + if let Some(v) = read_numeric(lc_root, &["lai", "tau_w"]) { + lai.tau_w = v; + } } fn apply_bioco2_overrides(bioco2: &mut crate::bioco2::BioCo2Prm, lc_root: &Value) { diff --git a/src/supy/_check.py b/src/supy/_check.py index 77aed582a..ee6498536 100644 --- a/src/supy/_check.py +++ b/src/supy/_check.py @@ -514,6 +514,9 @@ def check_state(df_state: pd.DataFrame, fix=True) -> List: list_col_rule = set(dict_rules_indiv.keys()).difference([ x.lower() for x in list_col_forcing ]) + list_col_rule_required = { + key for key in list_col_rule if not dict_rules_indiv[key].get("optional", False) + } # check the following: # 0. mandatory variables for the Fortran interface @@ -526,7 +529,7 @@ def check_state(df_state: pd.DataFrame, fix=True) -> List: # 1. correct columns col_df_state = df_state.columns.get_level_values("var") # 1.1 if all columns are present - set_diff = set(list_col_rule).difference(col_df_state) + set_diff = set(list_col_rule_required).difference(col_df_state) if len(set_diff) > 0: str_issue = f"Mandatory columns missing from df_state: {set_diff}" list_issues.append(str_issue) diff --git a/src/supy/checker_rules_indiv.json b/src/supy/checker_rules_indiv.json index 0391c49db..96b49b8e0 100644 --- a/src/supy/checker_rules_indiv.json +++ b/src/supy/checker_rules_indiv.json @@ -954,6 +954,48 @@ "param": {}, "unit": "1" }, + "w_wilt": { + "cat": "grid", + "logic": "NA", + "optional": true, + "param": {}, + "unit": "1" + }, + "w_opt": { + "cat": "grid", + "logic": "NA", + "optional": true, + "param": {}, + "unit": "1" + }, + "f_shape": { + "cat": "grid", + "logic": "NA", + "optional": true, + "param": {}, + "unit": "1" + }, + "w_on": { + "cat": "grid", + "logic": "NA", + "optional": true, + "param": {}, + "unit": "1" + }, + "w_off": { + "cat": "grid", + "logic": "NA", + "optional": true, + "param": {}, + "unit": "1" + }, + "tau_w": { + "cat": "grid", + "logic": "NA", + "optional": true, + "param": {}, + "unit": "1" + }, "lat": { "cat": "grid", "logic": "range", @@ -2231,4 +2273,4 @@ }, "unit": "1" } -} \ No newline at end of file +} diff --git a/src/supy/data_model/core/site.py b/src/supy/data_model/core/site.py index ae44cb689..4cf960f3e 100644 --- a/src/supy/data_model/core/site.py +++ b/src/supy/data_model/core/site.py @@ -400,15 +400,134 @@ class LAIParams(BaseModel): ) laitype: FlexibleRefValue(int) = Field( default=0, - description="LAI calculation choice (0: original, 1: new high latitude)", + description=( + "LAI calculation choice: 0 original (thermal-only), 1 new high latitude, " + "2 moisture-aware (Jarvis stress factor + CLM5-style persistence trigger; " + "PR1 scaffolding, no-op until PR2). See GH-1292." + ), json_schema_extra={ "unit": "dimensionless", "display_name": "LAI Calculation Method", }, ) + # --- GH-1292 moisture-aware phenology parameters (laitype=2) ------------------------------- + # Thresholds operate on dimensionless relative soil water w = 1 - smd/smdcap, not on SMD (mm). + w_wilt: Optional[FlexibleRefValue(float)] = Field( + default=None, + description=( + "Lower bound of Jarvis water-stress factor on relative soil water (moisture-aware LAI). " + "Effective only when laitype == 2; PR1 no-op." + ), + json_schema_extra={ + "unit": "dimensionless", + "display_name": "Jarvis lower relative soil water (w_wilt)", + }, + ) + w_opt: Optional[FlexibleRefValue(float)] = Field( + default=None, + description=( + "Upper bound of Jarvis water-stress factor on relative soil water. " + "Effective only when laitype == 2; PR1 no-op." + ), + json_schema_extra={ + "unit": "dimensionless", + "display_name": "Jarvis upper relative soil water (w_opt)", + }, + ) + f_shape: Optional[FlexibleRefValue(float)] = Field( + default=None, + description=( + "Shape exponent on the Jarvis water-stress factor. " + "Effective only when laitype == 2; PR1 no-op." + ), + json_schema_extra={ + "unit": "dimensionless", + "display_name": "Jarvis shape exponent (f_shape)", + }, + ) + w_on: Optional[FlexibleRefValue(float)] = Field( + default=None, + description=( + "Leaf-on persistence threshold on relative soil water (CLM5-style stress-deciduous). " + "Effective only when laitype == 2; PR1 no-op." + ), + json_schema_extra={ + "unit": "dimensionless", + "display_name": "Leaf-on persistence threshold (w_on)", + }, + ) + w_off: Optional[FlexibleRefValue(float)] = Field( + default=None, + description=( + "Leaf-off persistence threshold on relative soil water. " + "Effective only when laitype == 2; PR1 no-op." + ), + json_schema_extra={ + "unit": "dimensionless", + "display_name": "Leaf-off persistence threshold (w_off)", + }, + ) + tau_w: Optional[FlexibleRefValue(float)] = Field( + default=None, + description=( + "Persistence window for running-mean relative soil water. " + "Effective only when laitype == 2; PR1 no-op." + ), + json_schema_extra={ + "unit": "days", + "display_name": "Soil-water persistence window (tau_w)", + }, + ) ref: Optional[Reference] = None + # Defaults for the six moisture-aware fields; mirrored in Fortran LAI_PRM and Rust LaiPrm. + MOISTURE_DEFAULTS: ClassVar[Dict[str, float]] = { + "w_wilt": 0.15, + "w_opt": 0.40, + "f_shape": 1.0, + "w_on": 0.35, + "w_off": 0.20, + "tau_w": 15.0, + } + + @model_validator(mode="after") + def _check_moisture_thresholds(self) -> "LAIParams": + """GH-1292: moisture-aware threshold consistency; only enforced when laitype == 2.""" + + laitype_value = self.laitype.value if isinstance(self.laitype, RefValue) else self.laitype + if int(laitype_value) != 2: + return self + + def _resolve(field_name: str) -> float: + raw = getattr(self, field_name) + if raw is None: + return self.MOISTURE_DEFAULTS[field_name] + return raw.value if isinstance(raw, RefValue) else raw + + w_wilt_val = _resolve("w_wilt") + w_opt_val = _resolve("w_opt") + w_on_val = _resolve("w_on") + w_off_val = _resolve("w_off") + + if w_opt_val <= w_wilt_val: + raise ValueError( + f"LAIParams (laitype=2): w_opt ({w_opt_val}) must be greater than w_wilt ({w_wilt_val})" + ) + if w_off_val >= w_on_val: + raise ValueError( + f"LAIParams (laitype=2): w_off ({w_off_val}) must be less than w_on ({w_on_val})" + ) + + missing = [name for name in self.MOISTURE_DEFAULTS if getattr(self, name) is None] + if missing: + logger_supy.info( + "LAIParams (laitype=2): using defaults for unset moisture fields: %s", + ", ".join(missing), + ) + + return self + def to_df_state(self, grid_id: int, surf_idx: int) -> pd.DataFrame: """Convert LAI parameters to DataFrame state format. @@ -422,7 +541,7 @@ def to_df_state(self, grid_id: int, surf_idx: int) -> pd.DataFrame: # Adjust index for vegetation surfaces (surface index - 2) veg_idx = surf_idx - 2 - # Set basic LAI parameters + # Set basic LAI parameters (plus GH-1292 moisture-aware fields when laitype == 2). lai_params = { "baset": self.baset, "gddfull": self.gddfull, @@ -431,6 +550,12 @@ def to_df_state(self, grid_id: int, surf_idx: int) -> pd.DataFrame: "laimin": self.laimin, "laimax": self.laimax, "laitype": self.laitype, + "w_wilt": self.w_wilt, + "w_opt": self.w_opt, + "f_shape": self.f_shape, + "w_on": self.w_on, + "w_off": self.w_off, + "tau_w": self.tau_w, } cols = {("gridiv", "0"): grid_id} @@ -446,6 +571,7 @@ def to_df_state(self, grid_id: int, surf_idx: int) -> pd.DataFrame: "sddfull": 100.0, "laimax": self.LAIMAX_DF_DEFAULT, "laitype": 0, + **self.MOISTURE_DEFAULTS, } val = defaults.get(param, 0.0) cols[(param, f"({veg_idx},)")] = val @@ -498,6 +624,12 @@ def get_df_value(col_name: str, indices: Union[Tuple, int]) -> float: # Convert scalar parameters to RefValue lai_params = {key: RefValue(value) for key, value in lai_params.items()} + # GH-1292: pick up moisture-aware fields if present in df_state; legacy columns may be absent. + for field_name in ("w_wilt", "w_opt", "f_shape", "w_on", "w_off", "tau_w"): + col = (field_name, f"({veg_idx},)") + if col in df.columns: + lai_params[field_name] = RefValue(df.loc[grid_id, col]) + # Extract LAI power coefficients laipower = LAIPowerCoefficients.from_df_state(df, grid_id, veg_idx) diff --git a/test/data_model/test_laiparams_moisture_fields.py b/test/data_model/test_laiparams_moisture_fields.py new file mode 100644 index 000000000..e3ad192aa --- /dev/null +++ b/test/data_model/test_laiparams_moisture_fields.py @@ -0,0 +1,118 @@ +"""GH-1292 PR1 data-model tests for the moisture-aware LAI parameters. + +These tests lock the PR1 contract: (a) the six new optional fields +(`w_wilt`, `w_opt`, `f_shape`, `w_on`, `w_off`, `tau_w`) are serialised +through `to_df_state`/`from_df_state` with tolerant defaults; (b) the +threshold consistency validator fires only for `laitype == 2`; and (c) +`checker_rules_indiv.json` covers every new df_state column so that +`check_state()` does not reject a run. +""" + +from __future__ import annotations + +import json +from pathlib import Path + +import pytest +from pydantic import ValidationError + +from supy import SUEWSSimulation +from supy._check import check_state +from supy.data_model.core.site import LAIParams +from supy.data_model.core.type import RefValue + +pytestmark = pytest.mark.api + +VEG_IDX_EVETR = 0 # LAIParams.to_df_state uses surf_idx - 2, so veg_idx=0 maps to EVETR. + + +def _roundtrip(lai: LAIParams) -> LAIParams: + df = lai.to_df_state(grid_id=1, surf_idx=2) + return LAIParams.from_df_state(df, grid_id=1, surf_idx=2) + + +def test_defaults_loaded_for_legacy_config() -> None: + """LAIParams(laitype=0) without moisture fields round-trips to defaults.""" + + lai = LAIParams(laitype=0) + df = lai.to_df_state(grid_id=1, surf_idx=2) + for field in ("w_wilt", "w_opt", "f_shape", "w_on", "w_off", "tau_w"): + col = (field, f"({VEG_IDX_EVETR},)") + assert col in df.columns, f"missing df_state column {col}" + assert df.loc[1, col] == LAIParams.MOISTURE_DEFAULTS[field] + + recovered = LAIParams.from_df_state(df, grid_id=1, surf_idx=2) + assert recovered.laitype.value == 0 + for field in LAIParams.MOISTURE_DEFAULTS: + value = getattr(recovered, field) + assert value is not None + actual = value.value if isinstance(value, RefValue) else value + assert actual == LAIParams.MOISTURE_DEFAULTS[field] + + +def test_round_trip_with_explicit_fields() -> None: + """Explicit laitype=2 with non-default moisture fields survives round-trip.""" + + lai = LAIParams( + laitype=2, + w_wilt=0.2, + w_opt=0.45, + f_shape=1.5, + w_on=0.4, + w_off=0.15, + tau_w=20.0, + ) + recovered = _roundtrip(lai) + assert recovered.laitype.value == 2 + assert pytest.approx(0.2) == recovered.w_wilt.value + assert pytest.approx(0.45) == recovered.w_opt.value + assert pytest.approx(1.5) == recovered.f_shape.value + assert pytest.approx(0.4) == recovered.w_on.value + assert pytest.approx(0.15) == recovered.w_off.value + assert pytest.approx(20.0) == recovered.tau_w.value + + +def test_validator_fires_only_for_laitype_2() -> None: + """w_opt <= w_wilt raises only under laitype=2; legacy laitype=0 ignores it.""" + + with pytest.raises(ValidationError, match="w_opt"): + LAIParams(laitype=2, w_wilt=0.4, w_opt=0.3) + + with pytest.raises(ValidationError, match="w_off"): + LAIParams(laitype=2, w_on=0.3, w_off=0.5) + + # Under laitype=0 the validator must stay silent. + LAIParams(laitype=0, w_wilt=0.4, w_opt=0.3, w_on=0.3, w_off=0.5) + + +def test_checker_rules_cover_new_columns() -> None: + """src/supy/checker_rules_indiv.json must list every new moisture column.""" + + rules_path = ( + Path(__file__).resolve().parents[2] + / "src" + / "supy" + / "checker_rules_indiv.json" + ) + rules = json.loads(rules_path.read_text(encoding="utf-8")) + for field in ("w_wilt", "w_opt", "f_shape", "w_on", "w_off", "tau_w"): + assert field in rules, f"checker_rules_indiv.json missing entry for {field}" + + +def test_check_state_accepts_legacy_df_without_moisture_columns() -> None: + """Legacy df_state inputs may omit the new moisture columns without failing validation.""" + + sim = SUEWSSimulation.from_sample_data() + df_state = sim._config.to_df_state() + moisture_cols = [ + (field, f"({veg_idx},)") + for field in ("w_wilt", "w_opt", "f_shape", "w_on", "w_off", "tau_w") + for veg_idx in range(3) + ] + df_legacy = df_state.drop(columns=moisture_cols) + + issues = check_state(df_legacy, fix=False) + issues_text = "\n".join(issues) + assert "Mandatory columns missing from df_state" not in issues_text + for field in ("w_wilt", "w_opt", "f_shape", "w_on", "w_off", "tau_w"): + assert field not in issues_text, issues_text diff --git a/test/test_moisture_phenology_site.py b/test/test_moisture_phenology_site.py new file mode 100644 index 000000000..439869388 --- /dev/null +++ b/test/test_moisture_phenology_site.py @@ -0,0 +1,157 @@ +"""Regression tests for the moisture-aware FLUXNET diagnostic harness.""" + +from __future__ import annotations + +import importlib.util +import sys +import types +from pathlib import Path + +import numpy as np +import pandas as pd +import pytest + +pytestmark = pytest.mark.api + +REPO_ROOT = Path(__file__).resolve().parents[1] +SCRIPT_PATH = REPO_ROOT / "scripts" / "verify" / "moisture_phenology_site.py" +SCRIPT_SPEC = importlib.util.spec_from_file_location( + "moisture_phenology_site", SCRIPT_PATH +) +moisture_phenology_site = importlib.util.module_from_spec(SCRIPT_SPEC) +assert SCRIPT_SPEC.loader is not None +SCRIPT_SPEC.loader.exec_module(moisture_phenology_site) + + +def test_load_forcing_preserves_subhourly_index_and_filters_year(monkeypatch) -> None: + """load_forcing() should rely on SuPy parsing and keep minute resolution intact.""" + + idx = pd.DatetimeIndex( + ["2010-01-01 00:00:00", "2010-01-01 00:30:00", "2011-01-01 00:00:00"], + name="datetime", + ) + fake_df = pd.DataFrame({"lai": [1.0, 1.2, 0.8]}, index=idx) + + fake_util = types.ModuleType("supy.util") + fake_util.read_forcing = lambda path, tstep_mod=300: fake_df + fake_supy = types.ModuleType("supy") + fake_supy.util = fake_util + + monkeypatch.setitem(sys.modules, "supy", fake_supy) + monkeypatch.setitem(sys.modules, "supy.util", fake_util) + + loaded = moisture_phenology_site.load_forcing(Path("dummy.txt"), year=2010) + + assert loaded.index.tolist() == idx[:2].tolist() + assert loaded.index[1].minute == 30 + + +def test_load_forcing_defaults_to_first_full_year(monkeypatch) -> None: + """Without --year, load_forcing() should select the first full calendar year.""" + + idx = pd.DatetimeIndex( + [ + "2009-12-31 23:55:00", + "2010-01-01 00:00:00", + "2010-06-01 00:00:00", + "2010-12-31 23:55:00", + "2011-01-01 00:00:00", + ], + name="datetime", + ) + fake_df = pd.DataFrame({"lai": [0.5, 1.0, 1.2, 1.4, 0.8]}, index=idx) + + fake_util = types.ModuleType("supy.util") + fake_util.read_forcing = lambda path, tstep_mod=300: fake_df + fake_supy = types.ModuleType("supy") + fake_supy.util = fake_util + + monkeypatch.setitem(sys.modules, "supy", fake_supy) + monkeypatch.setitem(sys.modules, "supy.util", fake_util) + + loaded = moisture_phenology_site.load_forcing(Path("dummy.txt"), year=None) + + assert loaded.index.tolist() == idx[1:4].tolist() + assert loaded.index.min().year == 2010 + assert loaded.index.max().year == 2010 + + +def test_run_scenario_uses_requested_forcing_and_dates(monkeypatch) -> None: + """run_scenario() must pass the caller's forcing and run window through to SuPy.""" + + calls: dict[str, object] = {} + forcing = pd.DataFrame( + {"lai": [1.0, 1.1]}, + index=pd.DatetimeIndex( + ["2012-07-01 00:00:00", "2012-07-01 00:05:00"], name="datetime" + ), + ) + state_columns = pd.MultiIndex.from_tuples(moisture_phenology_site.LAI_TYPE_COLS) + + class FakeScenario: + def __init__(self, state_init: pd.DataFrame | None = None) -> None: + self.state_init = ( + state_init + if state_init is not None + else pd.DataFrame([[0, 0, 0]], columns=state_columns) + ) + + @classmethod + def from_sample_data(cls): + return cls() + + @classmethod + def from_state(cls, state: pd.DataFrame): + calls["state"] = state.copy() + return cls(state_init=state.copy()) + + def update_forcing(self, forcing_data: pd.DataFrame): + calls["forcing"] = forcing_data.copy() + return self + + def run(self, start_date=None, end_date=None): + calls["start_date"] = start_date + calls["end_date"] = end_date + output_index = pd.MultiIndex.from_product( + [[1], pd.DatetimeIndex(["2012-07-01 00:00:00"])], + names=["grid", "datetime"], + ) + return types.SimpleNamespace( + SUEWS=pd.DataFrame({"LAI": [1.5]}, index=output_index) + ) + + fake_supy = types.ModuleType("supy") + fake_supy.SUEWSSimulation = FakeScenario + monkeypatch.setitem(sys.modules, "supy", fake_supy) + + start_date = forcing.index.min() + end_date = forcing.index.max() + result = moisture_phenology_site.run_scenario( + 2, forcing, start_date=start_date, end_date=end_date + ) + + assert calls["forcing"].equals(forcing) + assert calls["start_date"] == start_date + assert calls["end_date"] == end_date + assert (calls["state"].loc[:, moisture_phenology_site.LAI_TYPE_COLS] == 2).all().all() + assert list(result.columns) == ["LAI"] + + +def test_compute_metrics_uses_full_simulation_series_for_seasonality() -> None: + """Amplitude and timing should come from the full simulated LAI series, not obs gaps.""" + + lai_sim = pd.Series( + [0.0, 5.0, 0.0], + index=pd.date_range("2012-01-01", periods=3, freq="1D"), + ) + lai_obs = pd.Series( + [0.0, 0.0], + index=pd.DatetimeIndex(["2012-01-01", "2012-01-03"]), + ) + + metrics = moisture_phenology_site.compute_metrics(lai_sim, lai_obs) + + assert metrics["rmse_vs_obs"] == 0.0 + assert metrics["amplitude"] == 5.0 + assert metrics["green_up_doy"] == 2.0 + assert np.isfinite(metrics["brown_down_doy"]) diff --git a/test/test_moisture_phenology_sweep.py b/test/test_moisture_phenology_sweep.py new file mode 100644 index 000000000..f46c22a49 --- /dev/null +++ b/test/test_moisture_phenology_sweep.py @@ -0,0 +1,57 @@ +"""GH-1292 PR3: smoke test for the moisture-aware LAI parameter-sweep tool. + +This test guards the sweep CLI's public contract without paying the full +cost of running SUEWS for many parameter values: it runs a single-value +sweep (equivalent to one extra SUEWS invocation past the baseline), and +asserts that the expected JSON artefact is produced with the documented +schema. The test is tagged slow because even one extra sample run takes +roughly a minute. +""" + +from __future__ import annotations + +import json +import subprocess +import sys +from pathlib import Path + +import pytest + +pytestmark = pytest.mark.api + +REPO_ROOT = Path(__file__).resolve().parents[1] +SWEEP_SCRIPT = REPO_ROOT / "scripts" / "verify" / "moisture_phenology_sweep.py" + + +@pytest.mark.slow +def test_sweep_cli_emits_expected_json() -> None: + out_path = REPO_ROOT / ".context" / "gh1292" / "London-test" / "sweep_w_opt.json" + if out_path.exists(): + out_path.unlink() + + result = subprocess.run( + [ + sys.executable, + str(SWEEP_SCRIPT), + "--param", + "w_opt", + "--values", + "0.40", + "--site", + "London-test", + "--dry-start", + ], + cwd=str(REPO_ROOT), + check=False, + capture_output=True, + text=True, + ) + assert result.returncode == 0, f"sweep CLI failed:\n{result.stderr}" + assert out_path.exists(), f"expected {out_path} to be written" + + payload = json.loads(out_path.read_text(encoding="utf-8")) + assert payload["parameter"] == "w_opt" + assert isinstance(payload["baseline_mean_lai"], (int, float)) + assert payload["points"] and payload["points"][0]["value"] == 0.40 + required_keys = {"value", "mean_lai", "rmse_vs_baseline", "amplitude", "green_up_doy"} + assert required_keys.issubset(payload["points"][0].keys()) diff --git a/test/test_phenology_lai_type_scaffolding.py b/test/test_phenology_lai_type_scaffolding.py new file mode 100644 index 000000000..0d8c762b1 --- /dev/null +++ b/test/test_phenology_lai_type_scaffolding.py @@ -0,0 +1,167 @@ +"""GH-1292 PR2 end-to-end contract for moisture-aware LAI phenology. + +The ``LAIType = 2`` branch implements Design C (see +``dev-ref/design-notes/gh1292-moisture-phenology.md``): a Jarvis +water-stress factor on ``delta_GDD`` combined with a CLM5-style +persistence latch on the tau_w-day running mean of relative soil +water. Two invariants define correctness: + +1. **Graceful degradation at well-watered sites.** When the soil + stays above ``w_opt`` and ``wbar_id`` stays above ``w_off`` at all + times, Jarvis returns 1.0 and the latch never flips, so the + trajectory must equal the ``LAIType = 0`` baseline. +2. **Engaged moisture gate under aggressive thresholds.** When the + Jarvis interval or persistence thresholds are tightened, the gate + must strictly reduce simulated LAI. + +Under PR1 the branch was a documented no-op; PR2 activates the +numerics and the bit-identity assertion now holds only because the +bundled London sample is well-watered. +""" + +from __future__ import annotations + +import numpy as np +import pandas as pd +import pytest + +from supy import SUEWSSimulation + +pytestmark = pytest.mark.physics + +LAI_TYPE_COLS = [ + ("laitype", "(0,)"), + ("laitype", "(1,)"), + ("laitype", "(2,)"), +] +MOISTURE_PARAM_COLS = { + "w_wilt": [("w_wilt", f"({i},)") for i in range(3)], + "w_opt": [("w_opt", f"({i},)") for i in range(3)], + "w_on": [("w_on", f"({i},)") for i in range(3)], + "w_off": [("w_off", f"({i},)") for i in range(3)], +} + + +def _run_with_state_override(laitype: int, moisture_overrides: dict | None = None) -> pd.DataFrame: + """Return the simulated LAI DataFrame with all veg surfaces set to ``laitype``. + + ``moisture_overrides`` is an optional dict mapping parameter name + (``w_wilt``, ``w_opt``, ``w_on``, ``w_off``) to a scalar value that + is broadcast across all three vegetation surfaces. + """ + + sim = SUEWSSimulation.from_sample_data() + state = sim.state_init.copy() + for col in LAI_TYPE_COLS: + state.loc[:, col] = laitype + if moisture_overrides: + for param, value in moisture_overrides.items(): + for col in MOISTURE_PARAM_COLS[param]: + state.loc[:, col] = value + + scenario = SUEWSSimulation.from_state(state) + scenario.update_forcing(sim.forcing) + output = scenario.run() + return output.SUEWS[["LAI"]].copy() + + +@pytest.fixture(scope="module") +def lai_baseline_0() -> pd.DataFrame: + return _run_with_state_override(0) + + +def test_laitype_2_matches_laitype_0_in_well_watered_regime( + lai_baseline_0: pd.DataFrame, +) -> None: + """Default moisture thresholds + well-watered London sample -> bit-identical to LAIType=0.""" + + lai_variant_2 = _run_with_state_override(2) + assert lai_baseline_0.shape == lai_variant_2.shape + np.testing.assert_array_equal(lai_baseline_0.values, lai_variant_2.values) + + +def test_laitype_2_reduces_lai_under_aggressive_moisture_thresholds( + lai_baseline_0: pd.DataFrame, +) -> None: + """Aggressive Jarvis thresholds force f_w < 1 for part of the year -> lower mean LAI.""" + + aggressive = { + "w_wilt": 0.80, + "w_opt": 0.95, + "w_on": 0.90, + "w_off": 0.70, + } + lai_aggressive = _run_with_state_override(2, aggressive) + + baseline_mean = float(lai_baseline_0.mean().iloc[0]) + aggressive_mean = float(lai_aggressive.mean().iloc[0]) + + assert aggressive_mean < baseline_mean, ( + f"aggressive LAIType=2 mean ({aggressive_mean:.4f}) should be strictly " + f"below LAIType=0 baseline ({baseline_mean:.4f})" + ) + # Must still run to completion without NaN and within physical LAI bounds. + assert not lai_aggressive.isna().any().any() + assert (lai_aggressive >= 0.0).all().all() + + +def test_sample_run_laitype_1_completes_without_nan() -> None: + """Legacy laitype=1 path must still run to completion without NaNs.""" + + sim = SUEWSSimulation.from_sample_data() + output = sim.run() + lai = output.SUEWS[["LAI"]] + assert len(lai) > 0 + assert not lai.isna().any().any() + + +def test_unsupported_positive_laitype_falls_back_to_legacy_nonzero_branch() -> None: + """Unsupported positive laitype values should retain the pre-GH-1292 fallback semantics.""" + + lai_type_1 = _run_with_state_override(1) + lai_type_3 = _run_with_state_override(3) + + assert lai_type_1.shape == lai_type_3.shape + np.testing.assert_array_equal(lai_type_1.values, lai_type_3.values) + + +def test_laitype_2_dry_start_delays_spring_green_up() -> None: + """Zeroing the vegetation soil store at init forces the latch off; spring green-up must lag. + + Design C expectation (CLM5 stress-deciduous): a dry initial root-zone keeps + ``leaf_on_permitted`` false until the tau_w-day running mean of relative soil + water climbs above ``w_on``. Cumulative LAI over the first months of the + sample year is therefore strictly lower than the LAIType=0 baseline. + """ + + sim = SUEWSSimulation.from_sample_data() + state_baseline = sim.state_init.copy() + for col in LAI_TYPE_COLS: + state_baseline.loc[:, col] = 0 + + state_dry = sim.state_init.copy() + for col in LAI_TYPE_COLS: + state_dry.loc[:, col] = 2 + # Deplete the soil store on the vegetation surfaces (indices 2, 3, 4) so initial SMD is + # close to capacity; the data model enforces soilstore >= 10 mm, so we cannot use exactly + # zero, but 10 mm gives w ~ 0.07, well below the default w_wilt (0.15) and w_off (0.20). + for veg_idx in (2, 3, 4): + state_dry.loc[:, ("soilstore_surf", f"({veg_idx},)")] = 10.0 + + def _run(state): + scenario = SUEWSSimulation.from_state(state) + scenario.update_forcing(sim.forcing) + return scenario.run().SUEWS[["LAI"]] + + lai_baseline = _run(state_baseline) + lai_dry = _run(state_dry) + + # Compare cumulative LAI over the first 60 days; the dry-start variant should integrate + # to a strictly smaller total because the moisture latch blocks spring GDD accumulation. + baseline_total = float(lai_baseline.iloc[: 60 * 288].sum().iloc[0]) + dry_total = float(lai_dry.iloc[: 60 * 288].sum().iloc[0]) + assert dry_total < baseline_total, ( + f"dry-start LAIType=2 cumulative LAI ({dry_total:.1f}) should lag LAIType=0 " + f"baseline ({baseline_total:.1f}) when the latch is forced off on day 1" + ) + assert not lai_dry.isna().any().any()