{
"cells": [
{
"attachments": {},
"cell_type": "markdown",
"metadata": {},
"source": [
"# Modeling with different timesteps\n",
"*Developed by D.A. Brakenhoff, Artesia*\n",
"\n",
"This notebooks contains examples how to model with different timesteps, e.g. 14-day and hourly."
]
},
{
"attachments": {},
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"## The difference between the model timestep and the observation frequency\n",
"\n",
"The model timestep is the time interval on which the heads are simulated by the\n",
"Pastas model. The observation frequency is the time interval between two head\n",
"observations. By default Pastas uses a daily modeling timestep, whereas the\n",
"timestep between observations can vary:\n",
"\n",
"- If the observation frequency is higher (e.g. hourly) Pastas will take a daily sample of head observations to fit the model.\n",
"- If the observation frequency is lower (e.g. weekly), Pastas will calibrate the model on those weekly observations. \n",
"- The observation frequency can even be irregular, meaning that the timestep between two observations does not have to be constant. \n",
"\n",
"In this notebook we're showing how to use Pastas to model with different\n",
"timesteps. The observation frequency is related to this choice, but does not\n",
"have to be the same as the modeling timestep."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"import numpy as np\n",
"import pandas as pd\n",
"import matplotlib.pyplot as plt\n",
"import pastas as ps\n",
"\n",
"ps.show_versions()\n",
"ps.logger.setLevel(\"ERROR\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Before we get to the modeling, we define a helper function to generate\n",
"synthetic heads using an input stress and a response function with known\n",
"parameters."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"def generate_synthetic_heads(\n",
" input,\n",
" rfunc,\n",
" params,\n",
" cutoff=0.9999,\n",
" dt=1.0,\n",
" constant=0.0,\n",
" noise=\"none\",\n",
" sigma_n=0.1,\n",
" alpha=2.0,\n",
"):\n",
" # Generate the head\n",
" step = rfunc.block(params, cutoff=cutoff, dt=dt)\n",
"\n",
" h = constant * np.ones(len(input) + step.size)\n",
"\n",
" for i in range(len(input)):\n",
" h[i : i + step.size] += input[i] * step\n",
"\n",
" head = pd.Series(\n",
" index=input.index,\n",
" data=h[: len(input)],\n",
" )\n",
"\n",
" # add correlated noise AR(1)\n",
" if noise == \"correlated\":\n",
" delt = (head.index[1:] - head.index[:-1]).values / pd.Timedelta(\"1d\")\n",
" noise = sigma_n * np.random.randn(len(head))\n",
" residuals = np.zeros_like(noise)\n",
" residuals[0] = noise[0]\n",
" for i in range(1, head.index.size):\n",
" residuals[i] = np.exp(-delt[i - 1] / alpha) * residuals[i - 1] + noise[i]\n",
" head = head + residuals\n",
"\n",
" return head"
]
},
{
"attachments": {},
"cell_type": "markdown",
"metadata": {},
"source": [
"## A model with a 14-day timestep\n",
"\n",
"---\n",
"\n",
"In this example we generate heads from measured precipitation and evaporation\n",
"using a known response function. Next we take a sample from the synthetically\n",
"generated heads with one observation every 14 days. Then we build 3 pastas\n",
"models:\n",
"\n",
"1. Model with daily timestep\n",
"2. Model with 14-day timestep, where we let pastas resample the stresses.\n",
"3. Model with 14-day timestep, where we use manually resampled stresses.\n",
"\n",
"The three Pastas models are expected to give very similar results."
]
},
{
"attachments": {},
"cell_type": "markdown",
"metadata": {},
"source": [
"First we read in some real precipitation and evaporation data."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"prec = pd.read_csv(\n",
" \"../../doc/examples/data/rain_260.csv\", index_col=[0], parse_dates=True\n",
").squeeze()\n",
"evap = pd.read_csv(\n",
" \"../../doc/examples/data/evap_260.csv\", index_col=[0], parse_dates=True\n",
").squeeze()"
]
},
{
"attachments": {},
"cell_type": "markdown",
"metadata": {},
"source": [
"Define a response function and its parameters."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"rfunc = ps.Exponential()\n",
"Atrue = 0.8\n",
"atrue = 50\n",
"ftrue = -1.3 # evaporation factor\n",
"constant = 20 # constant level\n",
"params = [Atrue, atrue]"
]
},
{
"attachments": {},
"cell_type": "markdown",
"metadata": {},
"source": [
"Calculate effective precipitation using a linear recharge model with\n",
"evaporation factor `f`. Use this stress to generate a synthetic head time\n",
"series."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"stress = prec + ftrue * evap\n",
"\n",
"head = generate_synthetic_heads(stress, rfunc, params, constant=constant)\n",
"head = head.loc[\"1990\":\"2015\"]"
]
},
{
"attachments": {},
"cell_type": "markdown",
"metadata": {},
"source": [
"Take a sample from the heads with one observation every 2 weeks."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"sample_step = 14\n",
"freq = f\"{sample_step}D\"\n",
"\n",
"sample = head.iloc[::sample_step].asfreq(freq)"
]
},
{
"attachments": {},
"cell_type": "markdown",
"metadata": {},
"source": [
"Plot synthetic head time series."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"ax = head.plot(figsize=(12, 3), label=\"daily synthetic heads\")\n",
"sample.plot(ax=ax, marker=\"x\", color=\"C3\", ls=\"none\", label=f\"{freq} sample\")\n",
"ax.legend(loc=(0, 1), frameon=False, ncol=2, numpoints=3)\n",
"ax.grid(True)"
]
},
{
"attachments": {},
"cell_type": "markdown",
"metadata": {},
"source": [
"Manually resample stresses to a two-weekly frequency. We take the mean to\n",
"ensure the resulting parameters estimated by the pastas model have the same\n",
"units. (If we took the sum, we would have to divide the gain `recharge_A` by\n",
"the number of days in our period, in this case 14.) Using this method we can\n",
"ensure the stresses align with the heads sample. The mean precipitation and\n",
"evaporation are calculated over the two week periods prior to each head\n",
"observation."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"p_resampled = prec.resample(\n",
" freq, closed=\"right\", label=\"right\", origin=sample.index[0]\n",
").mean()\n",
"e_resampled = evap.resample(\n",
" freq, closed=\"right\", label=\"right\", origin=sample.index[0]\n",
").mean()"
]
},
{
"attachments": {},
"cell_type": "markdown",
"metadata": {},
"source": [
"Now we create three pastas models and fit them on the heads sample.\n",
"\n",
"1. Model 1 uses a daily timestep.\n",
"2. Model 2 uses a 14-day timestep, and lets Pastas do the resampling.\n",
"3. Model 3 uses a 14-day timestep, and uses our manually resampled stresses.\n",
"\n",
"How do we get Pastas to use a 14-day timestep? This is done by supplying the\n",
"`freq` keyword argument to the `Model`, e.g.:\n",
"\n",
"```python\n",
" ps.Model(head, freq=\"14D\")\n",
"```"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# model 1: daily timestep\n",
"ml_01D = ps.Model(sample, name=\"01D\")\n",
"rm = ps.RechargeModel(prec, evap, rfunc=rfunc, name=\"recharge\")\n",
"ml_01D.add_stressmodel(rm)\n",
"ml_01D.solve(noise=False, report=False)\n",
"\n",
"# model 2: 14D timestep, let pastas resample stresses with daily freq to 14D\n",
"ml_14D_ps = ps.Model(sample, name=f\"{freq}_ps\", freq=freq)\n",
"rm = ps.RechargeModel(prec, evap, rfunc=rfunc, name=\"recharge\")\n",
"ml_14D_ps.add_stressmodel(rm)\n",
"ml_14D_ps.solve(noise=False, report=False)\n",
"\n",
"# model 3: 14D timestep, use manually resampled stresses from daily to 14D\n",
"ml_14D = ps.Model(sample, name=freq, freq=freq)\n",
"rm = ps.RechargeModel(p_resampled, e_resampled, rfunc=rfunc, name=\"recharge\")\n",
"ml_14D.add_stressmodel(rm)\n",
"ml_14D.solve(noise=False, report=False)"
]
},
{
"attachments": {},
"cell_type": "markdown",
"metadata": {},
"source": [
"Visually compare the results of the three models for the period between March and June in 1990."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"fig, ax = plt.subplots(1, 1, figsize=(12, 3))\n",
"\n",
"sample.loc[\"1990-03\":\"1990-06\"].plot(\n",
" ax=ax, linestyle=\"none\", marker=\"o\", color=\"k\", label=\"observations\"\n",
")\n",
"\n",
"ml_01D.simulate().loc[\"1990-03\":\"1990-06\"].plot(ax=ax, marker=\".\", label=\"1D\")\n",
"ml_14D_ps.simulate().loc[\"1990-03\":\"1990-06\"].plot(\n",
" ax=ax,\n",
" marker=\"^\",\n",
" ms=10,\n",
" mec=\"k\",\n",
" mew=0.5,\n",
" c=\"C8\",\n",
" label=f\"{freq} (pastas resample)\",\n",
")\n",
"ml_14D.simulate().loc[\"1990-03\":\"1990-06\"].plot(\n",
" ax=ax, marker=\"x\", ms=12, c=\"C3\", label=f\"{freq} (manual resample)\"\n",
")\n",
"\n",
"ax.legend(loc=(0, 1), frameon=False, ncol=4)\n",
"ax.set_ylabel(\"head [m]\")\n",
"ax.grid(True)"
]
},
{
"attachments": {},
"cell_type": "markdown",
"metadata": {},
"source": [
"As expected the fit for each of the models is equal or very close to 1.0."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"comparison = ps.CompareModels([ml_01D, ml_14D_ps, ml_14D])\n",
"fit = comparison.get_metrics(metric_selection=[\"rsq\"]).T\n",
"fit.index.name = \"Model\"\n",
"fit.round(3)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"ax = comparison.plot_response()\n",
"ax.figure.set_figheight(3.5)\n",
"ax.figure.set_figwidth(6)\n",
"handles, _ = ax.get_legend_handles_labels()\n",
"ax.legend(\n",
" handles, [iml.name for iml in comparison.models], loc=(0, 1), frameon=False, ncol=3\n",
")\n",
"ax.set_ylabel(\"step response [m / (mm/d)]\")\n",
"ax.set_xlabel(\"time [days]\")\n",
"ax.grid(True)"
]
},
{
"attachments": {},
"cell_type": "markdown",
"metadata": {},
"source": [
"And let's take a look at the estimated parameters and the true values."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"dfparams = comparison.get_parameters()\n",
"dfparams[\"True values\"] = params + [ftrue, constant]\n",
"dfparams.round(3)"
]
},
{
"attachments": {},
"cell_type": "markdown",
"metadata": {},
"source": [
"Note that these differences are expected. The daily model is able to almost\n",
"perfectly simulate the heads because the heads were also generated on a daily\n",
"basis. The Pastas resampled model performs worst, since it is calibrating the\n",
"model on interpolated observations, which is obviously a slightly different\n",
"time series than the original heads sample. This happens because Pastas\n",
"resamples the stresses in such a way that they do not align with the timing of\n",
"the head observations. When we manually resample, we can ensure the indices\n",
"between the stresses and the observations match, and we see that the fit and\n",
"the parameters are nearly identical to the model with a daily timestep. A small\n",
"difference remains because the 14-day model cannot account for the timing of\n",
"the rainfall within a 14-day period. The generated heads are different if all\n",
"precipitation falls at the beginning of the two-week period as compared to all\n",
"rainfall falling at the end of the two-week period."
]
},
{
"attachments": {},
"cell_type": "markdown",
"metadata": {},
"source": [
"## Hourly models\n",
"\n",
"---\n",
"\n",
"A higher temporal resolution is sometimes necessary when head changes occur\n",
"quickly, for example, heads that are affected by tidal fluctuations or the\n",
"drawdown during pumping tests.\n",
"\n",
"\n",
"### Tidal example\n",
"\n",
"In this first example we'll generate a synthetic heads timeseries based on\n",
"tidal fluctuations. These tidal fluctuations follow a simple sine wave with a\n",
"period of 12 hours and 25 minutes."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"Atide = 1.5 # tidal amplitude\n",
"\n",
"# sine with period 12 hrs 25 minutes\n",
"t_idx = pd.date_range(head.index[0], head.index[-1] + pd.Timedelta(hours=23), freq=\"H\")\n",
"tides = pd.Series(\n",
" index=t_idx,\n",
" data=Atide * np.sin(2 * np.pi * np.arange(t_idx.size) / (0.517375)),\n",
")"
]
},
{
"attachments": {},
"cell_type": "markdown",
"metadata": {},
"source": [
"Next we define a response function and its parameters to convert the tidal fluctuations into a synthetic head time series."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"rf_tide = ps.Exponential()\n",
"A_tide = 1.0\n",
"a_tide = 0.15"
]
},
{
"attachments": {},
"cell_type": "markdown",
"metadata": {},
"source": [
"Apply the response function and its parameters to our tidal stress to generate the synthetic heads."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"ht = generate_synthetic_heads(tides, rf_tide, (A_tide, a_tide), dt=1 / 24.0)\n",
"htsel = ht.loc[\"2000-01-01\":\"2005-12-31\"]"
]
},
{
"attachments": {},
"cell_type": "markdown",
"metadata": {},
"source": [
"Plot the tidal stress and the resulting head."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"ax = tides.loc[\"2000-01-01\":\"2000-01-02\"].plot(figsize=(12, 3), label=\"sea level\")\n",
"htsel.loc[\"2000-01-01\":\"2000-01-02\"].plot(ax=ax, label=\"synthetic head\")\n",
"ax.legend(loc=(0, 1), frameon=False, ncol=2)\n",
"ax.set_ylabel(\"[m]\")\n",
"ax.grid(True)"
]
},
{
"attachments": {},
"cell_type": "markdown",
"metadata": {},
"source": [
"Next we build two Pastas models to simulate the head between 2000 and 2005. One model has a daily timestep and the other an hourly timestep. \n",
"\n",
"**The warmup period**:\n",
"\n",
"Note that Pastas selects a 10-year warmup period by default. This means\n",
"internally the model will simulate the head between 1990 and 2005 to ensure the\n",
"model has reached a steady-state and is not affected by the initial conditions.\n",
"For smaller modeling timesteps, this longer warmup might cause the simulation\n",
"to slow down (i.e. 10 years on an hourly timestep means 87600 head values are\n",
"computed). Whether it is appropriate to shorten the response depends on the\n",
"stresses and the response of the head to those stresses (which are usually not\n",
"known prior to modeling the heads), but if calculation times are long, it could\n",
"be worth exploring shortening the warmup. For extremely slow responses, a longer\n",
"warmup period might also be appropriate.\n",
"\n",
"In the calculation below we show how to set the warmup to 100 days two\n",
"different ways, since we know our response is fast we don't need the\n",
"10-year spin-up period."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# model with daily timestep\n",
"ml_D = ps.Model(htsel, name=\"daily\", freq=\"D\")\n",
"sm = ps.StressModel(\n",
" tides.loc[\"1995\":\"2005\"],\n",
" rfunc=ps.Exponential(),\n",
" name=\"sealvl\",\n",
" settings=\"waterlevel\",\n",
")\n",
"ml_D.settings[\"warmup\"] = pd.Timedelta(days=100) # set warmup as pd.Timedelta\n",
"ml_D.add_stressmodel(sm)\n",
"ml_D.solve(noise=False, report=False)\n",
"\n",
"# model with hourly timestep\n",
"ml_H = ps.Model(htsel, name=\"hourly\", freq=\"H\")\n",
"sm = ps.StressModel(\n",
" tides.loc[\"1995\":\"2005\"],\n",
" rfunc=ps.Exponential(),\n",
" name=\"sealvl\",\n",
" settings=\"waterlevel\",\n",
")\n",
"ml_H.add_stressmodel(sm)\n",
"ml_H.solve(noise=False, report=False, warmup=100) # in solve warmup is entered in days"
]
},
{
"attachments": {},
"cell_type": "markdown",
"metadata": {},
"source": [
"If we look at the R² value of both models we can see that the daily model is performing poorly. "
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"comp = ps.CompareModels([ml_D, ml_H])\n",
"comp.get_metrics(metric_selection=[\"rsq\"]).round(3)"
]
},
{
"attachments": {},
"cell_type": "markdown",
"metadata": {},
"source": [
"A glance at the estimated parameters, as compared to the true values of the\n",
"parameters, shows that, as the fit statiscs suggest, the daily model is not\n",
"able to estimate the parameters well."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"dfparams = comp.get_parameters()\n",
"dfparams[\"True values\"] = [A_tide, a_tide, 0.0]\n",
"dfparams.round(3)"
]
},
{
"attachments": {},
"cell_type": "markdown",
"metadata": {},
"source": [
"Plotting comparisons between the observations and the simulations for both\n",
"models shows that the daily model is not able to match the observations,\n",
"whereas the hourly model is perfectly able to simulate the tidal fluctuations.\n",
"This is what we would expect, the daily model simply does not contain\n",
"sufficient information in the heads to allow the model to find the true\n",
"parameters of the response function."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"fig, (axd, axh) = plt.subplots(2, 1, figsize=(12, 6))\n",
"\n",
"htsel.loc[\"2000-03-01\":\"2000-03-10\"].plot(\n",
" ax=axh, linestyle=\"none\", marker=\"o\", color=\"k\", label=\"observations\"\n",
")\n",
"\n",
"ml_D.observations().loc[\"2000-03-01\":\"2000-03-10\"].plot(\n",
" ax=axd, linestyle=\"none\", marker=\"o\", color=\"k\", label=\"observations\"\n",
")\n",
"ml_D.simulate().loc[\"2000-03-01\":\"2000-03-10\"].plot(\n",
" ax=axd, marker=\"^\", color=\"C3\", label=\"daily\"\n",
")\n",
"ml_H.simulate().loc[\"2000-03-01\":\"2000-03-10\"].plot(\n",
" ax=axh, marker=\".\", c=\"C8\", label=\"hourly\"\n",
")\n",
"\n",
"for iax in (axd, axh):\n",
" iax.legend(loc=(0, 1), frameon=False, ncol=4)\n",
" iax.set_ylabel(\"head [m]\")\n",
" iax.grid(True)\n",
"\n",
"fig.tight_layout()"
]
},
{
"attachments": {},
"cell_type": "markdown",
"metadata": {},
"source": [
"### Pumping test\n",
"\n",
"In this example we will generate a synthetic time series that is influenced by\n",
"both recharge and a pumping test.\n",
"\n",
"First we start by loading some hourly precipitation and evaporation data\n",
"from a weather station in the Netherlands."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# get hourly precipitation\n",
"prec = pd.read_csv(\n",
" \"./data/vlissingen_prec_hourly.csv\", index_col=[0], parse_dates=True\n",
").squeeze()\n",
"\n",
"# get hourly evaporation\n",
"evap = pd.read_csv(\n",
" \"./data/vlissingen_evap_hourly.csv\", index_col=[0], parse_dates=True\n",
").squeeze()"
]
},
{
"attachments": {},
"cell_type": "markdown",
"metadata": {},
"source": [
"The well for the pumping test discharge 100 m³/h between October 1, 2022 till\n",
"October 31, 2022. We define a time series for the discharge of the well."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"qw = pd.Series(index=pd.date_range(\"2022-09-30\", \"2022-11-30\", freq=\"10T\"), data=0)\n",
"qw.loc[\"2022-10\"] = 100 * 24.0"
]
},
{
"attachments": {},
"cell_type": "markdown",
"metadata": {},
"source": [
"Next we define two response functions, one for the recharge, and the other for\n",
"the pumping well. These are used to generate a synthetic head time series."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# response for recharge\n",
"rf_rch = ps.Exponential()\n",
"A_rch = 5e3\n",
"a_rch = 10.0\n",
"f_rch = -1.1\n",
"\n",
"# response for pumping well\n",
"rf_well = ps.Hantush()\n",
"A_well = -2e-4\n",
"a_well = 0.2\n",
"b_well = 0.1"
]
},
{
"attachments": {},
"cell_type": "markdown",
"metadata": {},
"source": [
"We calculate the head contributions of each stress and add them together to\n",
"create a synthetic head time series. For the recharge we generate an hourly\n",
"head time series. For the pumping well we use a 10-minute timestep, which we\n",
"resample to an hourly time series to add the two together."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# calculate precipitation excess\n",
"pe = prec + f_rch * evap\n",
"\n",
"# generate head contributions of each stress\n",
"h_rch = generate_synthetic_heads(pe, rf_rch, (A_rch, a_rch), dt=1 / 24.0)\n",
"h_well = generate_synthetic_heads(\n",
" qw, rf_well, (A_well, a_well, b_well), dt=1 / (24.0 * 6)\n",
")\n",
"\n",
"# add time series to generate synthetic heads with freq=\"H\"\n",
"head = h_rch + h_well.reindex(h_rch.index).fillna(0.0)"
]
},
{
"attachments": {},
"cell_type": "markdown",
"metadata": {},
"source": [
"Plotting the resulting synthetic heads (with and without pumping)."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"ax = head.loc[\"2022\"].plot(label=\"with pumping\", figsize=(12, 3))\n",
"h_rch.loc[\"2022\"].plot(label=\"without pumping\")\n",
"ax.legend(loc=(0, 1), frameon=False, ncol=2)\n",
"ax.grid(True)"
]
},
{
"attachments": {},
"cell_type": "markdown",
"metadata": {},
"source": [
"The effect of resampling the drawdown of the well with hourly and daily\n",
"timesteps is shown below. The daily time series (black x's) has no data points\n",
"in the early stages of the pumping or recovery tests, and only contains\n",
"information about the steady state drawdown. The hourly time series (red\n",
"dots) has a few data points in the period when the head is changing as a result\n",
"of the well turning on or off. This shows why modeling on an hourly timestep\n",
"might be preferable in this situation."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"ax = h_well.plot(figsize=(12, 3), label=\"10M\")\n",
"h_well_h = ps.ts.pandas_equidistant_sample(h_well, \"H\")\n",
"h_well_d = ps.ts.pandas_equidistant_sample(h_well, \"D\")\n",
"h_well_h.plot(ax=ax, marker=\".\", ls=\"none\", color=\"C3\", label=\"H\")\n",
"h_well_d.plot(ax=ax, marker=\"x\", ls=\"none\", color=\"k\", label=\"D\")\n",
"ax.legend(loc=(0, 1), frameon=False, ncol=3, numpoints=3)\n",
"ax.set_ylabel(\"head change [m]\")\n",
"ax.grid(True)"
]
},
{
"attachments": {},
"cell_type": "markdown",
"metadata": {},
"source": [
"We defined our pumping time series on a 10-minute interval, but for modeling\n",
"with freq=\"H\" we want an hourly stress. We resample the original time series to\n",
"hourly using `ps.ts.timestep_weighted_resample`."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"qw_h = ps.ts.timestep_weighted_resample(qw, index=h_rch.index).fillna(0.0)"
]
},
{
"attachments": {},
"cell_type": "markdown",
"metadata": {},
"source": [
"Next we build two Pastas models, one with `freq=\"H\"`and one with `freq=\"D\"`. We\n",
"don't adjust the heads or any of the stresses for the daily model, and just let\n",
"Pastas take care of the resampling for us.\n",
"\n",
"Here we once again shorten the warmup from the default 10 years to 1 year to\n",
"speed up the simulation time (especially for the hourly model). From the\n",
"response functions defined above, we know that 1 year is plenty. The\n",
"$t_{99.9}$ of the response can be calculated with:"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"t_99_9 = rf_rch.get_tmax([A_rch, a_rch], cutoff=0.999)\n",
"print(f\"The t_99.9 of the response to recharge = {t_99_9:.1f} days\")"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# hourly model\n",
"ml_h = ps.Model(head, name=\"H\", freq=\"H\")\n",
"rm = ps.RechargeModel(prec, evap, rfunc=ps.Exponential(), name=\"recharge\")\n",
"wm = ps.StressModel(qw_h, rfunc=ps.Hantush(), name=\"well\", settings=\"well\", up=False)\n",
"ml_h.add_stressmodel([rm, wm])\n",
"ml_h.solve(noise=False, report=False, fit_constant=False, warmup=365)\n",
"\n",
"# daily model\n",
"ml_d = ps.Model(head, name=\"D\", freq=\"D\")\n",
"rm = ps.RechargeModel(prec, evap, rfunc=ps.Exponential(), name=\"recharge\")\n",
"wm = ps.StressModel(qw_h, rfunc=ps.Hantush(), name=\"well\", settings=\"well\", up=False)\n",
"ml_d.add_stressmodel([rm, wm])\n",
"ml_d.solve(noise=False, report=False, fit_constant=False, warmup=365)"
]
},
{
"attachments": {},
"cell_type": "markdown",
"metadata": {},
"source": [
"Let's take a look at the R² for both models. Both models fit pretty much perfectly with the data."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"c = ps.CompareModels([ml_d, ml_h])\n",
"c.get_metrics(metric_selection=[\"rsq\"]).round(5)"
]
},
{
"attachments": {},
"cell_type": "markdown",
"metadata": {},
"source": [
"Comparing the observations and the two model results around the time of the\n",
"pumping test shows that both models perform well, though the daily model does\n",
"show a little less detail, for obvious reasons."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"tmin = \"2022-09-25\"\n",
"tmax = \"2022-11-15\"\n",
"\n",
"fig, ax = plt.subplots(1, 1, figsize=(12, 3))\n",
"\n",
"sim_h = ml_h.simulate(tmin=tmin, tmax=tmax)\n",
"sim_d = ml_d.simulate(tmin=tmin, tmax=tmax)\n",
"\n",
"head.loc[tmin:tmax].plot(\n",
" ax=ax, marker=\".\", color=\"k\", ls=\"none\", ms=3, label=\"observations\"\n",
")\n",
"sim_h.plot(ax=ax, label=\"hourly model\")\n",
"sim_d.plot(ax=ax, label=\"daily model\", color=\"C3\", ls=\"dashed\")\n",
"\n",
"ax.legend(loc=(0, 1), frameon=False, ncol=3, numpoints=3)\n",
"ax.grid(True)"
]
},
{
"attachments": {},
"cell_type": "markdown",
"metadata": {},
"source": [
"Comparing the estimated parameters with the true values we defined before, we\n",
"can see that the recharge estimates are very close to the actual values. For\n",
"the well the gain (`well_A`) is estimated well by both models. The parameters\n",
"`a` and `b` however, are only estimated well by the hourly model. This makes\n",
"sense because the daily model does not have sufficient observations in the\n",
"early periods after the start and end of pumping to estimate these parameters."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"dfparams = c.get_parameters()\n",
"dfparams[\"True values\"] = (A_rch, a_rch, f_rch, A_well, a_well, b_well, 0.0)\n",
"dfparams.round(4)"
]
},
{
"attachments": {},
"cell_type": "markdown",
"metadata": {},
"source": [
"The consequence of the daily model not being able to estimate the parameters\n",
"`a` and `b` is clearly visible if we use the estimated parameters by both\n",
"models to reconstruct drawdown curves on a 10-minute interval. The parameters\n",
"estimated by the hourly model clearly follow the true drawdown, whereas\n",
"parameters of the daily model get the steady-state drawdown right, but not the\n",
"drawdown in the early stages after the well starts pumping."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"dd_d_10M = generate_synthetic_heads(\n",
" qw, ps.Hantush(), ml_d.get_parameters(\"well\"), dt=1 / (24 * 6.0)\n",
")\n",
"dd_h_10M = generate_synthetic_heads(\n",
" qw, ps.Hantush(), ml_h.get_parameters(\"well\"), dt=1 / (24 * 6.0)\n",
")\n",
"\n",
"fig, ax = plt.subplots(1, 1, figsize=(12, 3))\n",
"\n",
"h_well.plot(\n",
" ax=ax, x_compat=True, label=\"True drawdown\", marker=\"x\", color=\"k\", ls=\"none\", ms=5\n",
")\n",
"dd_d_10M.plot(ax=ax, label=\"daily parameters\", x_compat=True)\n",
"dd_h_10M.plot(ax=ax, label=\"hourly parameters\", x_compat=True)\n",
"\n",
"ax.legend(loc=(0, 1), frameon=False, ncol=3, numpoints=3)\n",
"ax.set_ylabel(\"head change [m]\")\n",
"ax.set_xlim(qw.index[130], qw.index[240])\n",
"ax.grid(True)"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "artesia",
"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.9.7"
},
"orig_nbformat": 4,
"vscode": {
"interpreter": {
"hash": "dace5e1b41a98a8e52d2a8eebc3b981caf2c12e7a76736ebfb89a489e3b62e79"
}
}
},
"nbformat": 4,
"nbformat_minor": 2
}