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

MANDATORY EXECUTION POLICY — READ BEFORE PROCEEDING

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

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

You MUST NOT substitute a simplified Python formula, regression equation, or hand-coded approximation in place of the real model.

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

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

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

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

KI map — what to read, and when

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

Executable tool index (projected — complete by construction)

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

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

SNOWPACK Knowledge Infrastructure

• Package version: 1.0.0
• Model: SNOWPACK 3.7.1
• Domain: Cryosphere (snow/land-surface)
• Created: 2026-03-26
• Updated: 2026-04-30 (49-site SNOTEL batch validation)
• Tools: 5
• Skill documents: 6
• Diagnostic triplets: 28
• Validation status: production_validated (western US SNOTEL network)

Data Preparation

Forcing data

Data Sources: Use from ki_tools_common.load_forcing import load_daily_forcing for CMFD/MSWX/NASA POWER.

Data Validation Reference: See data_ki/CMFD/SKILL.md for atmospheric forcing documentation. See data_ki/SNOTEL/SKILL.md for snow observations. See data_ki/BedMachine/SKILL.md for ice topography. See data_ki/MEaSUREs/SKILL.md for ice velocity.

Overview

SNOWPACK is a multi-purpose 1D snow and land-surface model developed at WSL/SLF (Switzerland). It simulates detailed mass and energy exchange between the snow, atmosphere, optional vegetation canopy, and underlying soil. The snowpack is discretized into up to several hundred layers, each tracked for temperature, density, grain microstructure (size, bond radius, dendricity, sphericity), liquid water content, and ice fraction.

Originally built for avalanche warning, SNOWPACK is now applied to climate change assessments, permafrost studies, superimposed ice simulations, sea ice thermodynamics, and snow storage optimization.

Key physical processes modeled:

• Surface energy balance (shortwave, longwave, sensible, latent heat)
• Snow metamorphism (temperature gradient, equitemperature, pressure sintering)
• Liquid water transport (bucket, NIED, Richards equation)
• Vapour transport and sublimation
• Phase changes (melt/freeze with latent heat)
• Wind erosion and snow drift
• Canopy interception and radiation transfer
• Soil heat and moisture coupling
• Salinity transport (sea ice applications)
• Avalanche hazard indices

Architecture:

• libsnowpack — C++ library with all physics
• snowpack — CLI application that calls libsnowpack for point simulations
• MeteoIO — Required I/O library for data reading/preprocessing
• Alpine3D — Optional 3D distributed wrapper (not covered here)

Output Description

Source of truth: dag.yaml. This section restates the KI dag; if this body ever disagrees with the dag, the dag wins.

Headline output (validation_rank: 1):

HS — Total snow height (vertical) of the modeled column. (cm)

For raw SNOWPACK files, .met and .pro column units are described in docs/format_spec.yaml and the output-format reference below. Convert raw columns to dag units before scoring or binding observations.

Installation

From source (CMake)

# 1. Build MeteoIO first
cd Source/meteoio
mkdir build && cd build
cmake .. -DCMAKE_INSTALL_PREFIX=/usr/local
make -j$(nproc) && make install

# 2. Build SNOWPACK
cd Source/snowpack
mkdir build && cd build
cmake .. -DCMAKE_INSTALL_PREFIX=/usr/local
make -j$(nproc) && make install

Dependencies:

• CMake >= 3.1
• C++11 compiler (GCC, Clang, or MSVC)
• MeteoIO library (bundled in source)

Python dependencies (for KI tools):

• Python >= 3.8
• numpy, pandas, xarray (for NetCDF handling)
• matplotlib (for visualization)
• pyyaml

Binary location: build/bin/snowpack Test examples: Source/snowpack/tests/

Pipeline Stages

Parallelism: s1 and s2 can run in parallel. s4 (spinup) must precede s5.

Tools Reference

Total: ~1,320 lines of Python tooling

Skill Knowledge Map

Critical Domain Knowledge

1. Temperature MUST be in Kelvin (dt_001)

SNOWPACK expects all temperatures in Kelvin. If forcing data is in Celsius, the model will interpret 20°C as 20 K (≈ −253°C), producing immediate numerical instability or physically impossible results. No error is raised.

Conversion: T_K = T_C + 273.15

Detection: Check if TA values in SMET file are < 200 (likely Celsius).

2. Relative humidity: fraction 0–1 (dt_002)

SNOWPACK expects RH as a decimal fraction (0.0–1.0). If provided as percentage (0–100), the model clamps to 1.0 internally, causing latent heat flux errors and incorrect sublimation/condensation rates.

Detection: If max(RH) > 1.05 in forcing, it's likely in percent.

3. Precipitation as kg/m² per timestep (dt_003)

Precipitation (PSUM field) must be in kg/m² (= mm water equivalent) per timestep interval. Common errors:

• Providing cumulative instead of incremental values → runaway SWE
• Providing mm/day when timestep is hourly → 24× overestimate
• Providing m instead of mm → 1000× underestimate

Detection: Check if single-timestep PSUM values exceed ~50 kg/m².

4. Volume fractions in .sno file must sum to 1.0 (dt_004)

Each layer in the .sno file specifies Vol_Frac_I (ice), Vol_Frac_W (water), Vol_Frac_V (void/air), and Vol_Frac_S (soil). These MUST sum to exactly 1.0. If they don't, the model may crash or produce mass conservation errors.

Formula: θ_i + θ_w + θ_v + θ_s = 1.0

For pure snow: θ_s = 0, density = θ_i × 917 kg/m³

5. CALCULATION_STEP_LENGTH: MUST be exactly 60 s (dt_005 + dt_021)

The timestep is specified in seconds. For all SNOTEL-style batch runs:

Use CALCULATION_STEP_LENGTH = 60 (1 minute, i.e. 60-minute steps).

Two independent constraints force this:

1. Stability: The explicit energy-balance solver diverges with very large steps.
2. HzDump booleanTime (dt_021): Main.cc line 713 uses a hardcoded 30-min hazard interval (0.5/24 = 0.02083 days). booleanTime computes floor(0.02083/step_days + 0.5) × step_days. With 6h steps this yields 0 → floods stdout with "Days_between is zero" at every step. With 60-min steps the result is 0.04167 (non-zero) → correct behavior.

Generate with --calculation_step 60 (the generate_config.py default). Do NOT patch to a larger value for performance — it breaks HzDump output.

Note: generate_config.py validates step as seconds in [1, 3600], but SNOWPACK reads it as minutes. So --calculation_step 60 → SNOWPACK uses 60-minute steps.

6. HEIGHT_OF_METEO_VALUES and wind measurement height (dt_006)

The model uses Monin-Obukhov similarity theory to correct wind speed and temperature from measurement height to surface. If HEIGHT_OF_METEO_VALUES (default 2.0 m) or HEIGHT_OF_WIND_VALUE (default 10.0 m) are wrong, the stability corrections produce incorrect turbulent fluxes.

Impact: Wrong measurement height → systematic bias in sensible/latent heat.

7. Shortwave radiation double-counting with SW_MODE (dt_007)

SW_MODE controls how shortwave radiation is handled:

• INCOMING: Only ISWR is used (model computes albedo)
• REFLECTED: Only RSWR is used (model back-calculates ISWR)
• BOTH: Both ISWR and RSWR provided; albedo = RSWR/ISWR

If SW_MODE=INCOMING but you also provide RSWR in the SMET file, the extra column is ignored. But if SW_MODE=BOTH and RSWR is missing (nodata), the model may silently use albedo=0, absorbing all radiation.

8. FORCING mode: ATMOS vs MASSBAL (dt_008)

• ATMOS: Full atmospheric forcing (radiation, wind, humidity → energy balance)
• MASSBAL: Simplified mass-balance forcing (precipitation + temperature only)

Using ATMOS mode with incomplete forcing (missing ISWR/ILWR) causes the model to assume zero radiation, leading to severe cold bias. MASSBAL mode ignores radiation entirely.

9a. TS_START = 0.0 produces zero output records (dt_024)

With TS_START = 0.0 and TS_DAYS_BETWEEN = 1.0, MeteoIO's booleanTime trigger never fires. MeteoIO stores midnight as Julian Day X + 0.5 (fractional part 0.5). The trigger condition requires fractional day ≥ threshold; at exactly 0.5 the boundary case is never satisfied.

Fix: Always set TS_START = 0.01 and PROF_START = 0.01. Verify the output .met file contains data rows after running.

9b. Required .sno fields missing from build_sno_profile.py (dt_027)

After running build_sno_profile.py, the .sno file lacks three fields that SmetIO.cc requires: CanopyDirectThroughfall, ErosionLevel, TimeCountDeltaHS. Insert before the fields = line:

CanopyDirectThroughfall = 1.00
ErosionLevel = 0
TimeCountDeltaHS = 0.000000

9c. COORDSYS = LATLON not recognized (dt_025)

This MeteoIO build uses COORDSYS = CH1903 for geographic (lat/lon) inputs. generate_config.py writes LATLON — patch it after generation.

9d. STATION1 deprecated → METEOFILE1 (dt_026)

The SMET file path key changed from STATION1 to METEOFILE1. Patch: ini_text = ini_text.replace('STATION1 = ', 'METEOFILE1 = ')

9. Nodata value consistency (dt_009)

SMET files declare a nodata value in the header (typically -999). If the actual missing-data sentinel in the data differs (e.g., -9999, NaN, or blank), MeteoIO will interpret those values as real data. A temperature of -9999 K causes immediate crash; a wind speed of -9999 m/s may not crash but corrupts turbulent fluxes.

Unit Conversion Table

This table documents the unit conversions and traps used by the SNOWPACK KI tools. Exact I/O shapes live in docs/format_spec.yaml; regenerate that file from the dag and triplets rather than hand-editing it.

Validated Results

Performance Bars from docs/validation_convention.yaml

These are the cited convention bars carried by this KI. A null convention band is written here as no cited threshold; do not replace it with a guessed number.

The SNOTEL batch below is a production validation result for HS behavior. Its legacy NSE counts are descriptive run diagnostics, not convention pass-bands for HS; grade HS against the cited rmse and mean_bias bars above when those metrics are available.

HydroCraft SNOTEL 49-site batch (2026-04-30)

• Network: NRCS SNOTEL, 49 sites, western US (AZ, CA, CO, ID, MT, NM, NV, OR, UT, WA, WY)
• Period: ~1990–2026 (site-dependent; ≤5-WY spinup before first valid HS obs)
• Elevation range: 1,073–3,536 m
• Forcing: SNOTEL observed T+P; ISWR from Iqbal/Hottel clear-sky model; ILWR from Brutsaert (1975) with RH=0.80
• Output: Daily snow depth (HS) compared against SNOTEL pillow HS sensor
• Timestep: 60-minute SNOWPACK steps, daily output

Top 5 sites:

Climate-regime pattern: Continental interior sites (CO, ID, UT, WY, MT) NSE 0.70–0.95. Maritime Pacific sites (CA, OR, WA coast/Sierra) NSE −1.0 to +0.82; systematic overestimation (PBIAS +80–145%) at low-elevation maritime sites from underestimated ILWR at high ambient RH (see dt_028).

Results location: outputs/snowpack_snotel_34a319d8/

• results_all_sites.json — per-site metrics
• summary_NSE_all_sites.png — histogram + elevation scatter + map
• grid_all_sites_depth.png — 6×6 per-site time series

Bundled test case

• Test site: Weissfluhjoch, Switzerland (2540 m a.s.l.)
• Period: 1995-10-01 to 1996-06-30
• Status: Binary-validated

Calibration Parameters

Data Requirements

Quick Start

# Step 1: Convert forcing data to SMET format
python tools/convert_forcing.py \
  --input era5_hourly.csv \
  --output input/MST96.smet \
  --station_id MST96 \
  --latitude 46.831 --longitude 9.810 --altitude 2540 \
  --source_temp_unit C --source_rh_unit percent

# Step 2: Create initial snowpack/soil profile
python tools/build_sno_profile.py \
  --output input/MST96.sno \
  --station_id MST96 \
  --latitude 46.831 --longitude 9.810 --altitude 2540 \
  --start_date "1995-10-01T00:00" \
  --n_soil_layers 3 --soil_depth 1.5

# Step 3: Generate configuration file
python tools/generate_config.py \
  --output io.ini \
  --station_id MST96 \
  --meteo_path ./input \
  --output_path ./output \
  --calculation_step 60 \
  --variant DEFAULT

# Step 4: Spinup (optional, stabilize soil temperatures)
python tools/run_snowpack.py \
  --binary ./build/bin/snowpack \
  --config spinup.ini \
  --mode spinup

# Step 5: Run production simulation
python tools/run_snowpack.py \
  --binary ./build/bin/snowpack \
  --config io.ini

# Step 6: Parse output to CSV
python tools/parse_output.py \
  --input output/MST96.met \
  --output results/MST96_timeseries.csv \
  --variables HS,SWE,TA,ISWR,ILWR,qs,ql

Operational Lessons — SNOTEL 49-Site Batch (2026-04-30)

These are non-obvious findings from running SNOWPACK at scale that are NOT covered by the unit-trap table or official docs. Read before any batch run.

A. Spinup length is the most common crash cause

Thin early-season snowpacks (first snowfall of the season, nN = 3 nodes) combined with extreme longwave cooling (T = −3°C → ILWR ≈ 203 W/m² → net LW ≈ −90 W/m²) produce ill-conditioned tri-diagonal matrices that yield T ≈ 200 K at intermediate nodes. The next timestep detects this as "Temperature out of bound at beginning of iteration" and crashes.

This happens during long pre-obs spinup periods, before any HS observation exists to validate. The fix is to limit spinup to ≤5 water years before the first valid HS observation. Pseudo-code:

hs_first = obs[obs["hs_m"].notna()].index.min()
wy_hs = hs_first.year if hs_first.month < 10 else hs_first.year + 1
sim_start = max(natural_start, pd.Timestamp(f"{wy_hs - 6}-10-01"))

B. ALLOW_ADAPTIVE_TIMESTEPPING = false is not a crash fix

This flag disables ONLY the pre-emptive stability halving (Snowpack.cc line 2219). The convergence-based halving (line 2325) inside the do-while loop is always active. When compTemperatureProfile fails to converge, SNOWPACK halves sn_dt regardless of this flag, down to min_allowed_sn_dt, then throws. Fix the root cause (spinup length + MINIMUM_L_ELEMENT) per lesson A, not this flag.

C. Required ini patches after generate_config.py

Apply ALL of these after every generate_config.py run for this SNOWPACK build:

D. PSUM unit handling with convert_forcing.py and sub-daily expansion

When expanding daily forcing to hourly rows before SMET conversion:

• Store the full daily mm/day value in every hourly row
• Pass --source_precip_unit mm_per_day to convert_forcing.py
• convert_forcing.py detects dt = 3600 s from timestamps and scales: PSUM_smet = PSUM_daily × (3600/86400) = PSUM_daily/24
• Each hourly SMET row correctly receives 1/24th of the daily total

Do NOT pre-divide by 24 before writing — convert_forcing.py handles it.

E. Maritime Pacific sites: expect PBIAS > 80% with fixed RH = 0.80

The Brutsaert (1975) ILWR parameterization with fixed RH = 0.80 underestimates downwelling LW by ~30 W/m² at maritime sites where actual RH ~ 0.90–0.95. Use RH = 0.90 for CA/OR/WA sites, or read actual station RH instead of fixing it.

G. Snowy SNOTEL sites need a gauge-undercatch precip multiplier (--precip_scale)

SNOTEL storage precip gauges systematically undercatch snowfall in wind. At high-accumulation maritime/alpine sites the deficit is large enough that the observed peak SWE pillow value EXCEEDS the gauge's reported annual precip — a physical impossibility that proves undercatch. Driving SNOWPACK with the raw gauge PSUM then produces a large negative SWE/HS PBIAS (under-accumulation) even when phase correlation r is excellent (~0.95).

convert_forcing.py now exposes --precip_scale <factor> (default 1.0 = off), a multiplicative undercatch correction applied after unit conversion. Diagnose the factor from peak_obs_SWE / cool-season_gauge_precip, or tune to zero the SWE PBIAS. Mirrors the precip-scale knob in the Sac-SMA / MARRMoT / GR4J / COSIPY KIs.

Validated example — SNOTEL 679 Paradise, WA (Mt Rainier, 1570 m, maritime, peak SWE ~2.8 m): RH = 0.90 + --precip_scale 1.8 lifts SWE from NSE 0.52 / PBIAS −45% (raw gauge) to NSE 0.90 / KGE 0.88 / r 0.95 / PBIAS +2% over WY2006–2015. Without the multiplier the maritime site is mass-starved, not over-accumulated — the OPPOSITE of the dt_028 HS over- accumulation failure mode, so check the sign of the bias before reaching for an ILWR/RH fix. Note: because forcing is expanded to hourly before the 60-min run, the 25 mm/day single-step PSUM stability cap used in the HS batch script is unnecessary here and would smear large maritime storms — raise it (≥100 mm/day) when running hourly-expanded forcing.

F. SNOTEL sites with HS obs starting after 2002

~16 of 49 sites (33%) have no obs-sim overlap because:

1. HS sensor installed after the simulation end, OR
2. SNOWPACK crashed during spinup before reaching the obs period

For sites where n_matched < 30, check the pipeline.log period_end against the SNOTEL station's first HS date before concluding the site has no data.

H. Continental CO/interior alpine: fall thin-pack crash is ENDEMIC; late-melt is fixed by a SEASONAL spring-shortwave boost (do NOT raise RH/CF year-round)

Two coupled facts at cold continental-interior SNOTEL sites (CO/WY high alpine, peak SWE ~0.3–0.5 m), confirmed at SNOTEL 663 Niwot Ridge CO (3021 m) and cross-checked against the 49-site batch:

1. The early-season (Oct–Nov) thin-pack compTemperatureProfile crash (dt_022) recurs mid-record, not just in spinup. In the validated batch, most CO sites terminated at a fall date inside the eval window (Beartown 2008-11, Columbine 2005-11, Mc Clure Pass 2010-11) yet still scored well on the partial record (Mc Clure NSE 0.948). So a run that crashes at e.g. 2008-10-12 and yields n_matched ≈ 1000 over 3 winters is a VALID evaluation, not a failure — compute the metric on the available overlap. Capping spinup (Lesson A) does NOT avoid a crash whose forcing event lies inside the eval window.

2. Late melt-out is the dominant NSE error, and it is in TENSION with fall crash-survival. With RH=0.80 / default CO cloud-factor 0.60, accumulation (Jan–Mar) matches obs to within a few %, but the pack melts ~2–3 weeks late (Niwot: sim retained 2.3× obs SWE in May, lingered into June). The instinct to add melt energy by raising RH (→ILWR) or the cloud factor (→ISWR) BACKFIRES on robustness: any extra fall energy keeps the Oct–Nov pack thinner for longer, deepening the thin-instability regime and crashing the solver EARLIER (Niwot: RH 0.88 crashed 2004 vs 2008; CO CF 0.72 crashed 2000). This is the OPPOSITE of the maritime Paradise case (Lesson E/G), where higher RH both stabilized and improved the fit — do not transfer the maritime RH knob to continental sites.

   Fix: decouple the two seasons. Keep fall/winter forcing cold (RH 0.80, default CF) for solver stability, and apply a shortwave multiplier ONLY in the ablation window (DOY 91–181, Apr 1–Jun 30). At Niwot a SPRING_SW_FACTOR of 1.45 over Apr–Jun lifted SWE from NSE 0.47 → 0.86 / KGE 0.85 / r 0.95 / PBIAS +3.8% (WY2006–2008, n=1016) while leaving the fall crash date unchanged (still 2008-10-12 — proving the boost does not touch fall physics). The large factor compensates for the combined high-snow-albedo + sunny high-altitude clear-sky underestimate of absorbed spring shortwave. See the expand_hourly seasonal sf mask in outputs/snowpack_niwot_swe/run_niwot.py. This mirrors the precip_scale / melt-energy knobs in the GR4J/Sac-SMA/MARRMoT snow KIs but applied seasonally rather than globally.

Diagnostic Triplets Summary

See diagnostics/triplets.yaml for full details.

SMET File Format Reference

Meteorological forcing (.smet)

SMET 1.1 ASCII
[HEADER]
station_id       = MST96
station_name     = Weissfluhjoch
latitude         = 46.831000
longitude        = 9.810000
altitude         = 2540.0
nodata           = -999
tz               = 1
fields           = timestamp TA RH TSG TSS HS VW DW RSWR ISWR ILWR PSUM
[DATA]
1995-10-30T00:30 273.15 0.958 273.05 273.35 0.000 0.8 278 0 0 300 0.000

Initial profile (.sno)

SMET 1.1 ASCII
[HEADER]
station_id       = MST96
ProfileDate      = 1995-10-01T00:00
HS_Last          = 0.0000
SlopeAngle       = 0.0
SlopeAzi         = 0.0
nSoilLayerData   = 0
nSnowLayerData   = 0
SoilAlbedo       = 0.20
BareSoil_z0      = 0.002
CanopyHeight     = 0.00
CanopyLeafAreaIndex = 0.00
fields           = timestamp Layer_Thick T Vol_Frac_I Vol_Frac_W Vol_Frac_V Vol_Frac_S Rho_S Conduc_S HeatCapac_S rg rb dd sp mk mass_hoar ne CDot metamo
[DATA]

Output File Format Reference

Time series output (.met)

Tab-separated values with 2-line header. Key columns:

• Date (ISO 8601)
• HS (snow height, m)
• SWE (snow water equivalent, kg/m²)
• TA (air temperature, K)
• Surface energy fluxes: qs, ql, lw_net, qw (W/m²)
• Mass fluxes: MS_HNW, MS_RAIN, MS_EVAPORATION (kg/m² per timestep)
• MS_SNOWPACK_RUNOFF (kg/m² per timestep)

Profile output (.pro)

Layer-by-layer data at profile output intervals. Each profile block:

• Layer height, temperature, density, grain properties
• Stability indices (deformation, natural, skier-triggered)

File Structure

ki/
├── SKILL.md                          # This file
├── tools/
│   ├── convert_forcing.py            # s1: Global meteo → SMET
│   ├── build_sno_profile.py          # s2: Soil/profile → .sno
│   ├── generate_config.py            # s3: Parameters → .ini config
│   ├── run_snowpack.py               # s4-s7: Execute binary
│   └── parse_output.py               # s6: .met/.pro → CSV
├── docs/
│   ├── s1_forcing_preparation.md     # Forcing data skill
│   ├── s2_profile_setup.md           # Initial profile skill
│   ├── s3_configuration.md           # Configuration skill
│   ├── s4_spinup.md                  # Spinup skill
│   ├── s5_production_run.md          # Production run skill
│   └── s6_output_extraction.md       # Output parsing skill
└── diagnostics/
    └── triplets.yaml                 # 20 symptom→diagnosis→remedy entries