PySnow

A modular Python snow model with switchable parameterizations for rapid testing and development.

PySnow combines snow physics from CLM (Community Land Model), VIC, and SNOWCLIM into a flexible framework where different parameterization schemes can be mixed and matched.

Tested accuracy: 2-3% error vs SNOTEL observations at multiple Western US sites.

Quick Start

Installation

git clone https://github.com/reedmaxwell/pysnow.git
cd pysnow
pip install -e .

# Optional: for SNOTEL site runner
pip install hf_hydrodata subsettools

Basic Usage

from pysnow import PySnowModel, SnowConfig
import numpy as np

# Configure with wet-bulb rain-snow partitioning (recommended)
config = SnowConfig(
    rain_snow_method='wetbulb',
    albedo_method='snicar_lite',
    phase_change_method='hierarchy',
    thin_snow_damping=0.5,  # Protects early-season thin snow
)

# Create model and run
model = PySnowModel(config)
outputs = model.run(forcing_data, dt=3600.0)

# Results
print(f"Peak SWE: {np.max(outputs['swe']):.1f} mm")

Run Against Any SNOTEL Site

# List available sites
python scripts/run_snotel_site.py --list-sites --state CO

# Run model comparison for a site
python scripts/run_snotel_site.py --site "Berthoud Summit" --state CO --wy 2023

# Run by site ID
python scripts/run_snotel_site.py --site_id "428:CA:SNTL" --wy 2024

This fetches CW3E forcing from HydroData, runs 5 rain-snow methods, and compares to SNOTEL observations.

Rain-Snow Partitioning Methods

Method	Description	Best For	Reference
wetbulb	Wet-bulb temperature threshold	Dry mountains (recommended)	Wang et al. 2019
wetbulb_sigmoid	Wet-bulb with smooth transition	Dry mountains	Wang et al. 2019
jennings	Humidity-dependent logistic	General use	Jennings et al. 2018
dai	Hyperbolic tangent	Global applications	Dai 2008
linear	Linear ramp (CLM default)	Humid regions	CLM
threshold	Simple temperature cutoff	Testing only	-

Why Wet-Bulb Methods?

In dry Western US mountains, falling hydrometeors cool via evaporation. Their surface temperature is closer to wet-bulb temperature (Tw) than air temperature (Ta). At Ta=2°C with 50% RH, Tw≈-1°C, so snow persists even above freezing.

# Recommended configuration for Western US
config = SnowConfig(
    rain_snow_method='wetbulb',
    rain_snow_params={'tw_threshold': 274.15},  # 1°C wet-bulb threshold
)

Thin Snow Damping

Early-season thin snowpacks can melt unrealistically fast during brief warm spells. The thin_snow_damping parameter provides thermal buffering:

config = SnowConfig(
    thin_snow_damping=0.5,      # 0-1, fraction of energy absorbed
    thin_snow_threshold=50.0,   # mm SWE, applies below this depth
)

• 0.0 = No damping (original behavior)
• 0.5 = Recommended (50% energy reduction for very thin snow)
• 0.7 = Aggressive protection for sites with known early-season issues

Configuration Reference

Full Configuration Options

from pysnow import SnowConfig

config = SnowConfig(
    # Rain-snow partitioning
    rain_snow_method='wetbulb',    # 'threshold', 'linear', 'jennings', 'dai', 'wetbulb', 'wetbulb_sigmoid'
    rain_snow_params={},           # Method-specific parameters

    # Albedo
    albedo_method='snicar_lite',   # 'constant', 'clm_decay', 'vic', 'essery', 'tarboton', 'snicar_lite'
    albedo_params={},

    # Phase change
    phase_change_method='hierarchy',  # 'hierarchy' (SNOWCLIM), 'clm'
    cold_content_tax=False,           # SNOWCLIM cold content limitation
    thin_snow_damping=0.5,            # Energy damping for thin snow [0-1]
    thin_snow_threshold=50.0,         # Threshold for damping [mm]

    # Thermal conductivity
    conductivity_method='jordan',     # 'jordan', 'sturm', 'calonne'

    # Compaction
    compaction_method='essery',       # 'essery', 'anderson', 'simple'

    # Ground heat flux
    ground_flux_method='constant',
    ground_flux_const=2.0,            # W/m²
    t_soil=278.15,                    # Soil temperature [K]

    # Turbulent fluxes
    stability_correction=True,
    z_wind=10.0,                      # Wind measurement height [m]
    z_temp=2.0,                       # Temperature measurement height [m]
    z_0=0.001,                        # Roughness length [m]

    # Liquid water
    max_water_frac=0.10,
    drain_rate=2.778e-5,              # m/s (10 cm/hr)
)

Forcing Data Format

PySnow expects forcing as a numpy array with shape (n_timesteps, 8):

Column	Variable	Units
0	Shortwave radiation (down)	W/m²
1	Longwave radiation (down)	W/m²
2	Precipitation rate	mm/hr
3	Air temperature	K
4	Wind U-component	m/s
5	Wind V-component	m/s
6	Air pressure	Pa
7	Specific humidity	kg/kg

Output Variables

Key	Description	Units
swe	Snow water equivalent	mm
depth	Snow depth	m
t_snow	Snow temperature	K
albedo	Surface albedo	-
density	Bulk snow density	kg/m³
melt	Melt rate	mm/timestep
runoff	Runoff	mm/timestep
snowfall	Snowfall	mm/timestep
rainfall	Rainfall	mm/timestep
sublimation	Sublimation (negative=loss)	mm/timestep
Rn, H, LE, G	Energy fluxes	W/m²

Example: Compare Methods

import numpy as np
from pysnow import PySnowModel, SnowConfig

# Load your forcing data
forcing = np.loadtxt('forcing.txt')

# Define configurations
configs = {
    'Wetbulb': SnowConfig(
        rain_snow_method='wetbulb',
        albedo_method='snicar_lite',
        thin_snow_damping=0.5,
    ),
    'CLM-linear': SnowConfig(
        rain_snow_method='linear',
        albedo_method='clm_decay',
    ),
    'Jennings': SnowConfig(
        rain_snow_method='jennings',
        albedo_method='snicar_lite',
    ),
}

# Run comparison
for name, config in configs.items():
    model = PySnowModel(config)
    result = model.run(forcing, dt=3600.0)
    print(f"{name}: Peak SWE = {np.max(result['swe']):.1f} mm")

Scripts

Script	Purpose
scripts/run_snotel_site.py	Run model against any SNOTEL site
scripts/run_comparison.py	Compare multiple configurations
scripts/diagnose_early_season.py	Debug early-season accumulation
scripts/diagnose_energy_balance.py	Detailed energy/mass balance

Validation Results

Tested against SNOTEL observations (Water Year 2023):

Site	Location	Observed	Wetbulb Model	Error
CSS Lab	Sierra Nevada, CA	1024 mm	1000 mm	-2.3%
Berthoud Summit	Front Range, CO	533 mm	522 mm	-2.1%
Canyon	Yellowstone, WY	389 mm	320 mm	-17.7%*

*Canyon error attributed to precipitation undercatch in CW3E forcing data.

Module Structure

pysnow/
├── __init__.py       # Package exports
├── constants.py      # Physical constants
├── rain_snow.py      # Rain-snow partitioning (6 methods)
├── albedo.py         # Albedo schemes (6 methods)
├── thermal.py        # Thermal properties & phase change
├── compaction.py     # Density evolution
├── turbulent.py      # Turbulent heat fluxes
├── radiation.py      # Radiation balance
└── model.py          # Main model driver

scripts/
├── run_snotel_site.py
├── run_comparison.py
├── diagnose_early_season.py
└── diagnose_energy_balance.py

docs/
└── SESSION_SUMMARY_2025-01-11.md

References

• Dai, A. (2008). Temperature and pressure dependence of the rain-snow phase transition. J. Climate.
• Jennings, K.S., et al. (2018). Spatial variation of the rain-snow temperature threshold. Nature Communications.
• Wang, Y., et al. (2019). Improving snow simulations with wet-bulb temperature. GRL.
• Stull, R. (2011). Wet-bulb temperature from relative humidity and air temperature. J. Appl. Meteor. Climatol.

License

MIT License

Contributing

Contributions welcome! Please open an issue or pull request.

Citation

PySnow: A modular Python snow model for parameterization testing.
https://github.com/reedmaxwell/pysnow