This tutorial covers a simulation example of groundwater flow and groundwater discharge with the GroundwaterDupuitPercolator component of Landlab. Simulation is run on steady state over a one layer aquifer and results are plotted on charts and grids.

The tutorial covers the following steps:

Definition of the aquifer and initial water table elevation
Setup of the hydraulic parameters and groundwater flow component
Groundwater simulation for 800 time steps and store partial flows
Representation of the final water table elevation and partial and total fluxes

Tutorial

Input data

You can download the input data from this link.

Code

#import required files
import matplotlib.pyplot as plt
import numpy as np

from landlab import RasterModelGrid, imshow_grid
from landlab.components import FlowAccumulator, GroundwaterDupuitPercolator

#define that raster model grid is closed by default
boundaries = {"top": "closed", "bottom": "closed", "right": "closed", "left": "closed"}
grid = RasterModelGrid((51, 51), xy_spacing=10.0, bc=boundaries)

#set up the node 1 as constant head boundary condition
grid.status_at_node[1] = grid.BC_NODE_IS_FIXED_VALUE

#define surface elevation
elev = grid.add_zeros("node", "topographic__elevation")
elev[:] = 0.0045 * (grid.x_of_node + grid.y_of_node) + 10

#define model base elevation
base = grid.add_zeros("node", "aquifer_base__elevation")

#define initial water table
wt = grid.add_zeros("node", "water_table__elevation")
wt[:] = 0.0035 * (grid.x_of_node + grid.y_of_node) + 9

#plot surface
plt.figure()
imshow_grid(grid, "topographic__elevation", shrink = 0.5)

#plot aquifer base
plt.figure()
imshow_grid(grid, "aquifer_base__elevation", shrink = 0.5)

#plot initial water table
plt.figure()
imshow_grid(grid, "water_table__elevation", cmap="Blues", shrink = 0.5)

#define hydraulic parameters
K = 0.01  # hydraulic conductivity, (m/s)
R = 1.0e-7 # recharge rate, (m/s)
n = 0.2  # porosity, (-)

#set up GroundwaterDupuitPercolator module for groundwater flow
gdp = GroundwaterDupuitPercolator(
    grid, hydraulic_conductivity=K, porosity=n, recharge_rate=R, regularization_f=0.01
)


fa = FlowAccumulator(
    grid,
    surface="topographic__elevation",
    flow_director="FlowDirectorSteepest",
    runoff_rate="surface_water__specific_discharge",
)

# set up the time steps and step length 
N = 200 #time steps
dt = 25 #seconds

#initial array to store values
recharge_flux = np.zeros(N)
gw_flux = np.zeros(N)
sw_flux = np.zeros(N)
storage = np.zeros(N)
s0 = gdp.calc_total_storage()

#run the model
for i in range(N):
    gdp.run_one_step(dt)
    fa.run_one_step()

    storage[i] = gdp.calc_total_storage()
    recharge_flux[i] = gdp.calc_recharge_flux_in()
    gw_flux[i] = gdp.calc_gw_flux_out()
    sw_flux[i] = gdp.calc_sw_flux_out()

#plot final water table
plt.figure()
imshow_grid(grid, "water_table__elevation", cmap="Blues", shrink = 0.5)

#plot water balance 
t = np.arange(0, N * dt, dt)

plt.figure(figsize=(8, 6))
plt.plot(
    t / 3600,
    np.cumsum(gw_flux) * dt + np.cumsum(sw_flux) * dt + storage - s0,
    "b-",
    linewidth=3,
    alpha=0.5,
    label="Total Fluxes + Storage",
)
plt.plot(
    t / 3600,
    np.cumsum(recharge_flux) * dt - recharge_flux[0] * dt,
    "k:",
    label="recharge flux",
)
plt.plot(t / 3600, np.cumsum(gw_flux) * dt, "b:", label="groundwater flux")
plt.plot(t / 3600, np.cumsum(sw_flux) * dt, "*", label="surface water flux", alpha=0.1)
plt.plot(t / 3600, storage - s0, "*", label="storage", alpha=0.5)
plt.ylabel("Cumulative Volume $[m^3]$")
plt.xlabel("Time [h]")
plt.semilogy()
plt.legend(frameon=False)
plt.show()