Modelamiento de flujo de agua subterránea usando la aproximación de Dupuit con Python y Landlab - Tutorial

Este tutorial cubre un ejemplo de simulación de flujo de agua subterránea y descarga de agua subterránea con el componente GroundwaterDupuitPercolator de Landlab. La simulación se ejecuta en régimen uniforme sobre un acuífero de una capa y los resultados se representan en gráficos y mapas.

El tutorial cubre los siguientes pasos:

Definición del acuífero y elevación inicial del nivel freático
Configuración de los parámetros hidráulicos y el componente de flujo de agua subterránea
Simulación de aguas subterráneas para múltiples pasos de tiempo y almacenar flujos parciales
Representación de la elevación final del nivel freático y flujos parciales y totales

Tutorial

Datos de entrada

Puedes descargar los datos de entrada desde este enlace.

Código

#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()