{ "cells": [ { "attachments": {}, "cell_type": "markdown", "id": "9f197b14", "metadata": {}, "source": [ "# Modeling snow\n", "*R.A. Collenteur, University of Graz / Eawag, November 2021*\n", "\n", "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." ] }, { "cell_type": "code", "execution_count": null, "id": "591e3e93", "metadata": {}, "outputs": [], "source": [ "import matplotlib.pyplot as plt\n", "import pandas as pd\n", "\n", "import pastas as ps\n", "\n", "ps.set_log_level(\"ERROR\")\n", "ps.show_versions()" ] }, { "cell_type": "markdown", "id": "f30f37d4", "metadata": {}, "source": [ "## 1. Load data\n", "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 temperature regularly drops below zero in winter. Given this observation, we may expect precipitation to (partially) fall as snow during these periods." ] }, { "cell_type": "code", "execution_count": null, "id": "21de2bf9", "metadata": {}, "outputs": [], "source": [ "head = pd.read_csv(\"data/heby_head.csv\", index_col=0, parse_dates=True).squeeze()\n", "evap = pd.read_csv(\"data/heby_evap.csv\", index_col=0, parse_dates=True).squeeze()\n", "prec = pd.read_csv(\"data/heby_prec.csv\", index_col=0, parse_dates=True).squeeze()\n", "temp = pd.read_csv(\"data/heby_temp.csv\", index_col=0, parse_dates=True).squeeze()\n", "\n", "ps.plots.series(head=head, stresses=[prec, evap, temp]);" ] }, { "cell_type": "markdown", "id": "8f1e3750", "metadata": {}, "source": [ "## 2. Make a simple model\n", "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." ] }, { "cell_type": "code", "execution_count": null, "id": "fb943327", "metadata": {}, "outputs": [], "source": [ "# Settings\n", "tmin = \"1985\" # Needs warmup\n", "tmax = \"2010\"" ] }, { "cell_type": "code", "execution_count": null, "id": "026ee16b", "metadata": {}, "outputs": [], "source": [ "ml1 = ps.Model(head)\n", "sm1 = ps.RechargeModel(\n", " prec, evap, recharge=ps.rch.FlexModel(), rfunc=ps.Gamma(), name=\"rch\"\n", ")\n", "ml1.add_stressmodel(sm1)\n", "\n", "# Solve the Pastas model in two steps\n", "ml1.solve(tmin=tmin, tmax=tmax, fit_constant=False, report=False)\n", "ml1.add_noisemodel(ps.ArNoiseModel())\n", "ml1.set_parameter(\"rch_srmax\", vary=False)\n", "ml1.solve(tmin=tmin, tmax=tmax, fit_constant=False, initial=False)\n", "ml1.plot(figsize=(10, 3));" ] }, { "cell_type": "markdown", "id": "0f3db43d", "metadata": {}, "source": [ "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 occurrence of snowfall and melt impacts the results.\n", "\n", "## 3. Account for snowfall and snow melt\n", "\n", "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." ] }, { "cell_type": "code", "execution_count": null, "id": "b55ecf45", "metadata": {}, "outputs": [], "source": [ "ml2 = ps.Model(head)\n", "sm2 = ps.RechargeModel(\n", " prec,\n", " evap,\n", " recharge=ps.rch.FlexModel(snow=True),\n", " rfunc=ps.Gamma(),\n", " name=\"rch\",\n", " temp=temp,\n", ")\n", "ml2.add_stressmodel(sm2)\n", "\n", "# Solve the Pastas model in two steps\n", "ml2.solve(tmin=tmin, tmax=tmax, fit_constant=False, report=False)\n", "ml2.add_noisemodel(ps.ArNoiseModel())\n", "ml2.set_parameter(\"rch_srmax\", vary=False)\n", "ml2.solve(tmin=tmin, tmax=tmax, fit_constant=False, initial=False)" ] }, { "cell_type": "markdown", "id": "93d48f8f", "metadata": {}, "source": [ "## Compare results\n", "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." ] }, { "cell_type": "code", "execution_count": null, "id": "b22218cc", "metadata": {}, "outputs": [], "source": [ "ax = ml2.plot(figsize=(10, 3))\n", "ml1.simulate().plot(ax=ax)\n", "plt.legend(\n", " [\n", " \"Observations\",\n", " \"Model w Snow NSE={:.2f}\".format(ml2.stats.nse()),\n", " \"Model w/o Snow NSE={:.2f}\".format(ml1.stats.nse()),\n", " ],\n", " ncol=3,\n", ")" ] }, { "cell_type": "markdown", "id": "a365d601", "metadata": {}, "source": [ "## Extract the water balance (States & Fluxes)" ] }, { "cell_type": "code", "execution_count": null, "id": "407489b6", "metadata": {}, "outputs": [], "source": [ "df = ml2.stressmodels[\"rch\"].get_water_balance(\n", " ml2.get_parameters(\"rch\"), tmin=tmin, tmax=tmax\n", ")\n", "df.plot(subplots=True, figsize=(10, 10));" ] }, { "cell_type": "markdown", "id": "997e386f", "metadata": {}, "source": [ "## References\n", "\n", "- 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)." ] } ], "metadata": { "kernelspec": { "display_name": "Python 3 (ipykernel)", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.10.12" }, "vscode": { "interpreter": { "hash": "29475f8be425919747d373d827cb41e481e140756dd3c75aa328bf3399a0138e" } } }, "nbformat": 4, "nbformat_minor": 5 }