Modelling a river in a cross-section

This notebook shows how to model a river with tidal fluctuations, based on an example from an old program called WATEX, which was used to derive the head beneath levee structures for input into geotechnical stability calculations.

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd

import timflow
import timflow.transient as tft

timflow.show_versions()

timflow version    : 0.2.0.dev0

Python version     : 3.13.3
Numpy version      : 2.4.3
Numba version      : 0.64.0
Scipy version      : 1.17.1
Pandas version     : 3.0.1
Matplotlib version : 3.10.8
LmFit version      : 1.3.4

System description

The system we want to model is shown below. It consists of 4 sections:

1. River: the summer bed of the river cuts deeper into the subsurface than the other sections

2. Foreland: the winter bed of the river due to higher water levels

3. Hinterland: the polder area protected by a dike or levee

4. Far hinterland: same as hinterland, but optionally with different parameters.

Aquifer properties and parameters

We will define some aquifer properties and parameters for each of the sections based on subsurface data. Note that we’re modeling the hinterland and far-hinterland as one section in this example.

In a subsequent step the geohydrological parameters will be calibrated.

# cross-sections
L_river = 75.00  # m, width river, though it extends to -infinity
L_foreland = 28.00  # m, width foreland
L_hinterland = 4000.00  # m, width hinterland, extends to +infinity

x0 = -(L_river + L_foreland)  # set x=0 at river bank

# layer elevations
ztop_riverbed = 1.00  # surface level in riverbed
ztop_foreland = 8.0  # surface level in foreland
ztop_hinterland = 8.0  # surface level in hinterland
ztop_aquifer = 0.0  # top of aquifer
zbot = -30.0  # bottom of aquifer

# geohydrology
kh = 10.0  # m/d, hydraulic conductivity
Saq = 1e-4  # 1/m, specific storage

Sll_river = 1e-6  # 1/m, specific storage aquitard
Sll_hinterland = 1e-3  # 1/m, specific storage aquitard

c_river = 100.0  # d, resistance riverbed
c_foreland = 150.0  # d, resistance foreland
c_hinterland = 500.0  # d, resistance hinterland

# set z-arrays for model
z_river = [ztop_riverbed, ztop_aquifer, zbot]
z_foreland = [ztop_foreland, ztop_aquifer, zbot]
z_hinterland = [ztop_hinterland, ztop_aquifer, zbot]

Load river levels

Load data from csv file containing observed water levels relative to the mean observed water level for an 8-day period. The file contains the river stage and measurements from 2 piezometers next to the river. The water level at the polder is set to a constant value of 0.0 (i.e. no changes caused by river stage fluctuations).

data = pd.read_csv("data/watex_example.csv", index_col=[0])
data.head()

	rivier	kruin	binnenteen	polder
T [days]
0.000000	0.000	0.000	0.000	0
0.020833	0.150	0.173	0.091	0
0.041667	0.160	0.177	0.091	0
0.062500	0.187	0.182	0.089	0
0.083333	0.197	0.185	0.085	0

data.plot(figsize=(10, 2), grid=True)

<Axes: xlabel='T [days]'>

Note that the observed fluctuations do not fluctuate around 0 for each of these locations. In order to model this system with timflow, we need the variation to be around 0.0, therefore, we normalize each time series by subtracting the mean.

data_norm = data - data.mean()
data_norm.plot(figsize=(10, 2), grid=True);

Build the cross-section model

Now we can build the cross-section model.

First, we need to define the river levels for the top boundary condition in the river and foreland sections. We want a list of time waterlevel pairs for the river water level.

We’re dealing with a river under the influence of tides, which means we have to be careful when specifying the river level. timflow maintains a given water level from a specific time until the next specified water level.

The following plot shows how timflow inteprets the water level for different pre-processing steps on the input data. If we just pass in the river time series as is, we’re continually leading the observed water level (orange line). We want to make sure the timflow input meets the observed level about at the midpoint (black line) of each time interval.

data_norm.loc[:1, "rivier"].plot(figsize=(10, 3), grid=True)
plt.step(
    data_norm.loc[:1].index,
    data_norm.loc[:1, "rivier"],
    where="post",
    label="no shift",
)
plt.step(
    data_norm.loc[:1].index,
    data_norm.loc[:1, "rivier"].shift(-1),
    where="post",
    label="shift -1",
)
mid = (data_norm.loc[:, "rivier"] + data_norm.loc[:, "rivier"].shift(-1).values).divide(2)
plt.step(
    mid.loc[:1].index, mid.loc[:1], where="post", label="midpoint", color="k", lw=1.0
)
plt.legend(loc=(0, 1), frameon=False, ncol=4);

# tsandhstar = data_norm["rivier"].to_frame().to_records().tolist()  # original
tsandhstar = mid.dropna().to_frame().to_records().tolist()  # midpoint
# tsandhstar = (
#     data_norm["rivier"].shift(-1).dropna().to_frame().to_records().tolist()
# )  # shift-1

ml = tft.ModelXsection(naq=1, tmin=1e-4, tmax=1e2)

river = tft.XsectionMaq(
    ml,
    -np.inf,  # river extends to infinitiy
    x0 + L_river,
    z=z_river,
    kaq=kh,
    Saq=Saq,
    Sll=Sll_river,
    c=c_river,
    topboundary="semi",
    tsandhstar=tsandhstar,
    name="river",
)
foreland = tft.XsectionMaq(
    ml,
    x0 + L_river,
    x0 + L_river + L_foreland,
    kaq=kh,
    z=z_foreland,
    Saq=Saq,
    Sll=Sll_river,
    c=c_foreland,
    topboundary="semi",
    tsandhstar=tsandhstar,
    name="foreland",
)
hinterland = tft.XsectionMaq(
    ml,
    x0 + L_river + L_foreland,
    np.inf,  # hinterland extends to infinity
    kaq=kh,
    z=z_hinterland,
    Saq=Saq,
    Sll=Sll_hinterland,
    c=c_hinterland,
    topboundary="semi",
    name="hinterland",
)

ml.solve()

self.neq  4
solution complete

Check the sections

river, foreland, hinterland

(river: XsectionMaq [-inf, -28.0] with h*(t),
 foreland: XsectionMaq [-28.0, 0.0] with h*(t),
 hinterland: XsectionMaq [0.0, inf])

Note that all sections are also stored in a dictionary where the names are used that were specified with the name keyword (in this case they were the same as the variable in which each section was stored).

ml.aq.inhomdict

{'river': river: XsectionMaq [-inf, -28.0] with h*(t),
 'foreland': foreland: XsectionMaq [-28.0, 0.0] with h*(t),
 'hinterland': hinterland: XsectionMaq [0.0, inf]}

Plot a cross-section of the model. Note the water level is just schematic in this picture.

fig, ax = plt.subplots(1, 1, figsize=(10, 3))
ax.set_ylim(-30, 24)
ml.plots.xsection(
    xy=[(x0 - 0, 0), (60, 0)],
    labels=True,
    params=False,
    names=True,
    hstar=10,
    sep="\n",
    ax=ax,
    units={"k": "m/d", "Ss": "m$^{-1}$", "c": "d"},
    fontsize=8,
)
ax.set_xlabel("x-coordinate [m]")
ax.set_ylabel("elevation [m]")
fig.suptitle("Cross-section across a river of a semi-confined aquifer")
fig.savefig("river_cross_section.png", dpi=150, bbox_inches="tight")

Show the summary of aquifer parameters for each section:

ml.aquifer_summary()

		layer	layer_type	k_h	c	S_s
	#
river	0	0	leaky layer	NaN	100.0	0.000001
	1	0	aquifer	10.0	NaN	0.0001
foreland	0	0	leaky layer	NaN	150.0	0.000001
	1	0	aquifer	10.0	NaN	0.0001
hinterland	0	0	leaky layer	NaN	500.0	0.001
	1	0	aquifer	10.0	NaN	0.0001

Plot head over time at several x-locations.

t = np.linspace(0.01, 8, 150)
xlocs = [
    -75.0,  # in river section
    1.0,  # at observation well 1
    47.5,  # at observation well 2
    100.0,  # in hinterland section
]

f, ax = plt.subplots(1, 1, figsize=(10, 3))
for i in range(len(xlocs)):
    h = ml.head(xlocs[i], 0.0, t)
    ax.plot(t, h[0], label=f"x={xlocs[i]:.1f} m")

ax.plot(data_norm.index, data_norm["rivier"], color="k", label="river", lw=1.0)
ax.legend(loc=(0, 1), frameon=False, ncol=5, fontsize="small")
ax.set_xlabel("time [d]")
ax.set_ylabel("head [m]")
ax.set_title("Simulated heads", loc="right")
ax.grid()

Plot head along x for multiple times

t = np.array([0.1, 1.0, 2.0, 3.0, 4.0, 5.0])
x = np.linspace(x0, 1000.0, 101)
y = np.zeros_like(x)

h = ml.headalongline(x, y, t)

f, ax = plt.subplots(1, 1, figsize=(10, 3))
for i in range(len(t)):
    ax.plot(x, h[0, i], label="t = {:.1f}".format(t[i]))
ax.legend(loc=(0, 1), frameon=False, ncol=8, fontsize="small")
ax.set_xlabel("x [m]")
ax.set_ylabel("head [m]")
ax.grid()
ax.set_title("Simulated heads", loc="right")
ax.set_xlim(x0, x.max());

Compare the model results to the observed heads in the two piezometers at \(x = 1.0\) m and \(x = 47.5\) m.

# piezometer locations in cross-section
xpb_kr = 1.0
xpb_bite = 47.5
xlocs = [xpb_kr, xpb_bite]

t = data_norm.index.to_numpy()

f, axes = plt.subplots(len(xlocs), 1, figsize=(10, 6), sharex=True, sharey=True)

for i in range(len(xlocs)):
    h = ml.head(xlocs[i], 0.0, t)
    axes[i].plot(t, h[0], label="model")
    hobs = data_norm.iloc[:, i + 1]
    axes[i].plot(hobs.index, hobs, label="observed")
    axes[i].set_title(f"x = {xlocs[i]} m", loc="right")
    axes[i].legend(loc=(0, 1), frameon=False, ncol=2)
    axes[i].grid()

Calibrate model on observed heads

The model simulation is reasonable at the second observation well , but it clearly isn’t as good at the river bank. Let’s see if we can improve the fit of the model to the observations by calibrating the aquifer parameters.

Since the data at \(t=0\) shows some jumps, we want to give the model some spin-up time before we start fitting to the data. We set a tstart value to select only observations after this moment.

The following parameters will be calibrated:

• \(k_h\) and \(S_s\) of the aquifer, this parameter should be the same in each section.

• resistance \(c\) of the confining units, this parameter can vary within each section.

• specific storage (\(S_s\)) for the leaky layers of the river and foreland, and separately the specific storage of the leaky layer in the hinterland.

In total this means we have 7 parameters that can vary in this calibration step.

tstart = 2.0  # model spin-up time, calibrate on timeseries after 2 days

# create calibrate object
cal = tft.Calibrate(ml)

# set observation time series
cal.series(
    "ObsWell 1",
    xpb_kr,
    0.0,
    layer=0,
    t=data_norm.loc[tstart:, "kruin"].index.to_numpy(),
    h=data_norm.loc[tstart:, "kruin"].values,
)
cal.series(
    "ObsWell 2",
    xpb_bite,
    0.0,
    layer=0,
    t=data_norm.loc[tstart:, "binnenteen"].index.to_numpy(),
    h=data_norm.loc[tstart:, "binnenteen"].values,
)

# define parameters for calibration
cal.set_parameter(
    "kaq",
    layers=0,
    initial=river.kaq[0],
    pmin=1.0,
    pmax=30.0,
    inhoms=("river", "foreland", "hinterland"),
)
cal.set_parameter(
    "Saq",
    layers=0,
    initial=river.Saq[0],
    pmin=1e-6,
    pmax=1e-2,
    inhoms=("river", "foreland", "hinterland"),
)
cal.set_parameter(
    "c",
    layers=0,
    initial=river.c[0],
    pmin=0.0,
    pmax=100.0,
    inhoms="river",
)
cal.set_parameter(
    "c",
    layers=0,
    initial=foreland.c[0],
    pmin=10.0,
    pmax=500.0,
    inhoms="foreland",
)
cal.set_parameter(
    "c",
    layers=0,
    initial=hinterland.c[0],
    pmin=10.0,
    pmax=1000.0,
    inhoms="hinterland",
)
cal.set_parameter(
    "Sll",
    layers=0,
    initial=river.Sll[0],
    pmin=1e-8,
    pmax=1e-4,
    inhoms=("river", "foreland"),
)
cal.set_parameter(
    "Sll",
    layers=0,
    initial=hinterland.Sll[0],
    pmin=1e-8,
    pmax=1e-1,
    inhoms="hinterland",
)

# run the calibration
cal.fit(report=False)
# cal.fit_least_squares(report=False, method="trf")

# print the RMSE
print(f"RMSE: {cal.rmse():.5f}")

....
...
...
...
...
...
...
...
...
....
...
...
...
...
....
Fit succeeded.
RMSE: 0.01790

Show calibration parameters results

cal.parameters.iloc[:, :-2]

	layers	optimal	std	perc_std	pmin	pmax	initial
kaq_0_0_river_foreland_hinterland	0	2.487442	11.149055	4.482137e+02	1.000000e+00	30.0000	10.000000
Saq_0_0_river_foreland_hinterland	0	0.000008	0.000233	2.905725e+03	1.000000e-06	0.0100	0.000100
c_0_0_river	0	21.932360	3307.509898	1.508050e+04	0.000000e+00	100.0000	100.000000
c_0_0_foreland	0	40.300383	569.471646	1.413068e+03	1.000000e+01	500.0000	150.000000
c_0_0_hinterland	0	216.770260	1525.802286	7.038799e+02	1.000000e+01	1000.0000	500.000000
Sll_0_0_river_foreland	0	0.000006	12.063255	2.086629e+08	1.000000e-08	0.0001	0.000001
Sll_0_0_hinterland	0	0.022505	0.276014	1.226436e+03	1.000000e-08	0.1000	0.001000

Plot the locations of observation wells.

fig, ax = plt.subplots(1, 1, figsize=(10, 3))
river.plot(labels=False, params=True, names=True, ax=ax, sep="\n", fmt=".1f")
foreland.plot(ax=ax, params=True, names=True, sep="\n", fmt=".1f")
hinterland.plot(ax=ax, params=True, names=True, sep="\n", fmt=".1f")
ml.elementlist[0].plot(ax=ax, hstar=10.0)
ml.elementlist[1].plot(ax=ax, hstar=10.0)
ax.set_xlim(x0, 100.0)
ax.set_ylim(-30, 27)
cal.xsection(ax=ax)
ax.set_ylabel("elevation [m]")
ax.set_xlabel("x-coordinate [m]")

Text(0.5, 0, 'x-coordinate [m]')

Plot the modeled vs. observed time series for both observation wells after calibration.

# piezometer locations in cross-section
xpb_kr = 1.0
xpb_bite = 47.5
xlocs = [xpb_kr, xpb_bite]

t = data_norm.index.to_numpy()

f, axes = plt.subplots(len(xlocs), 1, figsize=(10, 5), sharex=True, sharey=True)

for i in range(len(xlocs)):
    h = ml.head(xlocs[i], 0.0, t)
    axes[i].plot(t, h[0], label="timflow")
    hobs = data_norm.loc[tstart:].iloc[:, i + 1]  # showing calibration selection
    axes[i].plot(hobs.index, hobs, label="observed")
    axes[i].set_title(f"Observation well {i + 1} (x = {xlocs[i]} m)", loc="right")
    axes[i].legend(loc=(0, 1), frameon=False, ncol=2)
    axes[i].grid()

for iax in axes:
    iax.set_ylabel("head change [m]")
axes[-1].set_xlabel("time [d]")
axes[-1].figure.suptitle(
    "Calibration result of timflow cross-section model compared to observed heads."
)
axes[-1].figure.savefig("calibrated_head_response.png", bbox_inches="tight", dpi=150)

Let’s take a look at the calibrated aquifer parameters.

(
    ml.aquifer_summary()
    .style.format(na_rep="-", precision=2)
    .format("{:.2e}", subset=["S_s"])
)

		layer	layer_type	k_h	c	S_s
	#
river	0	0	leaky layer	-	21.93	5.78e-06
	1	0	aquifer	2.49	-	8.02e-06
foreland	0	0	leaky layer	-	40.30	5.78e-06
	1	0	aquifer	2.49	-	8.02e-06
hinterland	0	0	leaky layer	-	216.77	2.25e-02
	1	0	aquifer	2.49	-	8.02e-06

Discussion

The calibration result shows that the specific storage of the leaky layer in the hinterland section is orders of magnitude higher than the other specific storage terms. Also the resistance of the confining unit is higher in the foreland than it is in the polder.

These results have to be analyzed further to determine whether these results make physical sense. The fit to the observed data is reasonable, but the problem is likely over-parametrized. It is not possible to estimate all parameters separately, since there are infinite combinations of parameters that fit the data equally well. This causes the problem to be highly sensitive to the starting values of parameters, as well as the selected optimization method.

Some other points that could influence the results:

• Only two piezometers in the cross-section.

• The downwards trend in the head observations in the second piezometer probably isn’t ideal as timflow will never be able to simulate that trend with the current input data.

• The time series are supposedly presented relative to a mean observed head/water level, but somehow do not fluctuate around 0.0, which mean they had to be normalized (again) for the timflow simulation.

• Neglecting the storage of the leaky layers yields a physically improbable solution where all resistance is put underneath the river section, and the foreland and hinterland have a resistance on the minimum parameter boundary. However, it is usually not expected that the storage in the leaky layers has such a large influence on the outcome. In many models this parameter is neglected altogether.