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