Modeling snow#
R.A. Collenteur, University of Graz / Eawag, November 2021
In this notebook it is shown how to account for snowfall and smowmelt on groundwater recharge and groundwater levels, using a degree-day snow model. This notebook is work in progress and will be extended in the future.
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy.signal import fftconvolve
import pastas as ps
ps.set_log_level("ERROR")
ps.show_versions()
/tmp/ipykernel_1234/2988175969.py:2: DeprecationWarning:
Pyarrow will become a required dependency of pandas in the next major release of pandas (pandas 3.0),
(to allow more performant data types, such as the Arrow string type, and better interoperability with other libraries)
but was not found to be installed on your system.
If this would cause problems for you,
please provide us feedback at https://github.com/pandas-dev/pandas/issues/54466
import pandas as pd
Python version: 3.11.6
NumPy version: 1.26.4
Pandas version: 2.2.0
SciPy version: 1.12.0
Matplotlib version: 3.8.3
Numba version: 0.59.0
LMfit version: 1.2.2
Latexify version: Not Installed
Pastas version: 1.4.0
1. Load data#
In this notebook we will look at some data for a well near Heby, Sweden. All the meteorological data is taken from the E-OBS database. As can be observed from the temperature time series, the temparature regularly drops below zero in winter. Given this observation, we may expect precipitation to (partially) fall as snow during these periods.
head = pd.read_csv("data/heby_head.csv", index_col=0, parse_dates=True).squeeze()
evap = pd.read_csv("data/heby_evap.csv", index_col=0, parse_dates=True).squeeze()
prec = pd.read_csv("data/heby_prec.csv", index_col=0, parse_dates=True).squeeze()
temp = pd.read_csv("data/heby_temp.csv", index_col=0, parse_dates=True).squeeze()
ps.plots.series(head=head, stresses=[prec, evap, temp]);
2. Make a simple model#
First we create a simple model with precipitation and potential evaporation as input, using the non-linear FlexModel to compute the recharge flux. We not not yet take snowfall into account, and thus assume all precipitation occurs as snowfall. The model is calibrated and the results are visualized for inspection.
# Settings
tmin = "1985" # Needs warmup
tmax = "2010"
ml1 = ps.Model(head)
sm1 = ps.RechargeModel(
prec, evap, recharge=ps.rch.FlexModel(), rfunc=ps.Gamma(), name="rch"
)
ml1.add_stressmodel(sm1)
# As the evaporation used is a very rough estimation, vary k_v
ml1.set_parameter("rch_kv", vary=True)
# Solve the Pastas model in two steps
ml1.solve(tmin=tmin, tmax=tmax, noise=False, fit_constant=False, report=False)
ml1.set_parameter("rch_srmax", vary=False)
ml1.solve(tmin=tmin, tmax=tmax, noise=True, fit_constant=False, initial=False)
ml1.plot(figsize=(10, 3));
/home/docs/checkouts/readthedocs.org/user_builds/pastas/envs/v1.4.0/lib/python3.11/site-packages/pastas/timeseries_utils.py:90: FutureWarning: Day.delta is deprecated and will be removed in a future version. Use pd.Timedelta(obj) instead
if hasattr(offset, "delta"):
/home/docs/checkouts/readthedocs.org/user_builds/pastas/envs/v1.4.0/lib/python3.11/site-packages/pastas/stressmodels.py:1302: FutureWarning: 'A' is deprecated and will be removed in a future version, please use 'YE' instead.
if self.prec.series.resample("A").sum().max() < 12:
Fit report Head Fit Statistics
===================================================
nfev 53 EVP 50.36
nobs 590 R2 0.50
noise True RMSE 0.12
tmin 1985-01-01 00:00:00 AIC -3281.41
tmax 2010-01-01 00:00:00 BIC -3250.75
freq D Obj 1.11
warmup 3650 days 00:00:00 ___
solver LeastSquares Interp. No
Parameters (7 optimized)
===================================================
optimal initial vary stderr
rch_A 2.663685 0.577011 True ±16.84%
rch_n 0.610386 2.444868 True ±3.46%
rch_a 668.765296 82.680615 True ±25.68%
rch_srmax 123.797696 123.797696 False nan
rch_lp 0.250000 0.250000 False nan
rch_ks 0.823636 207.261970 True ±6.91%
rch_gamma 0.526596 0.404507 True ±14.68%
rch_kv 0.692416 0.893974 True ±4.36%
rch_simax 2.000000 2.000000 False nan
constant_d 77.156538 0.000000 False nan
noise_alpha 94.865119 1.000000 True ±7.14%
The model fit with the data is not too bad, but we are clearly missing the highs and lows of the observed groundwater levels. This could have many causes, but in this case we may suspect that the occurence of snowfall and melt impacts the results.
3. Account for snowfall and snow melt#
A second model is now created that accounts for snowfall and melt through a degree-day snow model (see e.g., Kavetski & Kuczera (2007). To run the model we add the keyword snow=True to the FlexModel and provide a time series of mean daily temperature to the RechargeModel. The temperature time series is used to split the precipitation into snowfall and rainfall.
ml2 = ps.Model(head)
sm2 = ps.RechargeModel(
prec,
evap,
recharge=ps.rch.FlexModel(snow=True),
rfunc=ps.Gamma(),
name="rch",
temp=temp,
)
ml2.add_stressmodel(sm2)
# As the evaporation used is a very rough estimation, vary k_v
ml2.set_parameter("rch_kv", vary=True)
# Solve the Pastas model in two steps
ml2.solve(tmin=tmin, tmax=tmax, noise=False, fit_constant=False, report=False)
ml2.set_parameter("rch_srmax", vary=False)
ml2.solve(tmin=tmin, tmax=tmax, noise=True, fit_constant=False, initial=False)
/home/docs/checkouts/readthedocs.org/user_builds/pastas/envs/v1.4.0/lib/python3.11/site-packages/pastas/timeseries_utils.py:90: FutureWarning: Day.delta is deprecated and will be removed in a future version. Use pd.Timedelta(obj) instead
if hasattr(offset, "delta"):
/home/docs/checkouts/readthedocs.org/user_builds/pastas/envs/v1.4.0/lib/python3.11/site-packages/pastas/stressmodels.py:1302: FutureWarning: 'A' is deprecated and will be removed in a future version, please use 'YE' instead.
if self.prec.series.resample("A").sum().max() < 12:
Fit report Head Fit Statistics
===================================================
nfev 30 EVP 54.78
nobs 590 R2 0.55
noise True RMSE 0.11
tmin 1985-01-01 00:00:00 AIC -3339.02
tmax 2010-01-01 00:00:00 BIC -3299.60
freq D Obj 1.00
warmup 3650 days 00:00:00 ___
solver LeastSquares Interp. No
Parameters (9 optimized)
===================================================
optimal initial vary stderr
rch_A 0.918317 0.783740 True ±6.43%
rch_n 1.358060 1.553572 True ±1.49%
rch_a 152.123265 110.816491 True ±6.86%
rch_srmax 16.437531 16.437531 False nan
rch_lp 0.250000 0.250000 False nan
rch_ks 123.677918 123.314886 True ±33.11%
rch_gamma 8.487820 6.887808 True ±10.38%
rch_kv 1.520442 1.118175 True ±4.95%
rch_simax 2.000000 2.000000 False nan
rch_tt 1.514985 1.641662 True ±4.03%
rch_k 2.684309 2.990550 True ±3.12%
constant_d 78.013947 0.000000 False nan
noise_alpha 80.825172 1.000000 True ±5.45%
Compare results#
From the fit_report we can already observe that the model fit improved quite a bit. We can also visualize the results to see how the model improved.
ax = ml2.plot(figsize=(10, 3))
ml1.simulate().plot(ax=ax)
plt.legend(
[
"Observations",
"Model w Snow NSE={:.2f}".format(ml2.stats.nse()),
"Model w/o Snow NSE={:.2f}".format(ml1.stats.nse()),
],
ncol=3,
)
<matplotlib.legend.Legend at 0x7f5c7e9bbd10>
Extract the water balance (States & Fluxes)#
df = ml2.stressmodels["rch"].get_water_balance(
ml2.get_parameters("rch"), tmin=tmin, tmax=tmax
)
df.plot(subplots=True, figsize=(10, 10));
References#
Kavetski, D. and Kuczera, G. (2007). Model smoothing strategies to remove microscale discontinuities and spurious secondary optima in objective functions in hydrological calibration. Water Resources Research, 43(3).