I have been having some troubles the heads in modflow thus I got this profiles in my river always and I cannot understand why the line is transform to curve. Here I attached the code and image with results
# Asignación directa desde dt (solo las columnas indicadas)
d_bkf = dt["d_bkf"][n]*2 # Profundidad del banco filtrante [m]
w_bkf = dt["w_bkf"][n] # Ancho del banco filtrante [m]
L = dt["L"][n] # Longitud del dominio [m]
B = dt["B"][n] # Ancho del dominio [m]
D = dt["D"][n]*2 # Profundidad activa del acuífero [m]
Ka = dt["Ka"][n] # Conductividad hidráulica horizontal [m/s]
Kr = dt["Kr"][n] # Conductividad hidráulica del lecho [m/s]
dr = dt["dr"][n] # Espesor de la zona (dr) [m]
Sy = dt["Sy"][n] # Porosidad efectiva [-]
Slope = dt["S"][n]
Ss = 1e-7 # Compresibilidad específica [1/m]
n_manning = dt["n"][n] # Rugosidad de Manning [-]
H_0 = dt["H_0"][n] # Nivel base inicial [m]
t_d_days = dt["t_d_days"][n] # Duración del evento [d]
deltay = Slope*L
Slope_trans = 1e-20
Kv_Kh = 1 # Kv/Kh Anisotropia en el ac
# Hydrograph
H_0 = H_0 # Nivel base (m)
eta_0 = (d_bkf-H_0)/H_0 # (d_bkf - H_0)/H_0 # H_p = 1*H_0 = 2m # coeffient, max stage of water. Evita que la creciente no sobrepase la banca llena.
t_d = 1*24*3600 # Duración del evento (s)
eta_d = 0.5 #Asimetria del hidrog. Valores peq. representan crecidas rapidas que bajan lento. explorar valores entre
H_p= eta_0 * H_0 #Elevanción máxima durante la creciente
nrow=int(dt["nrow"][n])
ncol=30
# nrow=15
# ncol=11
delc = L/nrow
delr = B/ncol
nlay = 5 # numero de capas
# Crear el modelo básico
path_results = "../results/simple_model/"
name_model = 'water_sim_bagre'
model = flopy.modflow.Modflow(name_model,
exe_name="MODFLOW-NWT.exe",
version="mfnwt",
model_ws=path_results,
verbose=True
)
hydrom_times = 600
t = np.arange(0, t_d, hydrom_times) # tmin, tmax, t_discretize
T = Ka * d_bkf # Transmisividad del acuífero
# Visualizar hidrograma
H_b = hydrograph_cooper_rorabaugh(t, H_0, eta_0, eta_d, t_d)
Q = (1/n_manning) * (H_b**(5/3)) * (Slope**(1/2)) * (w_bkf) # m3/s
idx_Qmax = np.argmax(Q)
velocity = Q / (w_bkf * d_bkf)
T_arrive = (L/velocity)/(3600*24)
v_mean = Q.mean()/(w_bkf*d_bkf)
v_max = Q.max()/(w_bkf*d_bkf)
t_mean = L/v_mean
t_max = L/v_max
t_canal = L / v_max
tc = phinum.get_tc(Sy=Sy, Ka=Ka, B=B, d_bkf=d_bkf)
a_celda = np.sqrt(delr * delc) # Dimensión de celda representativa
tc3 = Sy * a_celda**2 / (4 * T) # Tiempo característico tc2 tpo de difusión de la perturbación hidrúalica en el modelo.
t_pico = t[idx_Qmax] # Tiempo al pico del evento
# valores = np.array([t_canal, t_lecho, tc,t_pico, tc2]) # Agrupar todos los tiempos para análisis
valores = 1/np.array([t_canal, tc,t_pico, tc3])
TC = len(valores) / np.sum(valores) # Calcular TC como media armónica (pondera procesos rápidos)
delta_t_TC = 0.25 * TC # s # Calcular delta_t con factor de seguridad
reach_len = L / nrow # Evaluar delta_t que cumple Courant
delta_t_courant = .85*reach_len / v_max # Co = 1
delta_t = min(delta_t_TC, delta_t_courant) # Elegir el menor de ambos para respetar procesos rápidos y condición de estabilidad
courant_final = v_max * delta_t / reach_len # Calcular Courant final (referencia)
duracion_transitorio = 20 * 24 * 3600 # segundos # Definir duración total del periodo transitorio (ajustable más adelante)
# Número de pasos para el periodo transitorio
nstp_trans = int(np.ceil(duracion_transitorio / delta_t))
perlen = [1*24*3600, duracion_transitorio] # periodo estacionario y transitorio
nstp = [1, 400]
# Crear Topografía y caja del modelo
topo = generate_topo(nrow, ncol,L=L, delta_h=deltay, slope_t=Slope_trans, min_topo=D+d_bkf+1, plot=True)
botm = np.linspace(D, 0, nlay)
print(nstp)
# perlen = [.001*24*3600, 10*24*3600] # Periodos-longitud temporal
nper = perlen.__len__() # Numero de periodos
# nstp = [1,6000] # step per period
steady = [True, False] # Tipo de periodo
# cte = calculate_tsmult(perlen[1]/24/3600, nstp[1], 500, time_refine=(t_d/24/3600)*2)
# tsmult = [1, cte]
dis = flopy.modflow.ModflowDis(
model,
nlay=nlay, # Número de capas
nrow=nrow, # numero de columnas
ncol=ncol, # numero de filas
delr=delr,#esh_custom_unit.flatten(), # tamaño dy o entre filas
delc=delc, # tamaño dx o entre columnas
top=topo, # topografia o capa top
botm=botm, # Capas inferiores (nlay, nrow, ncol)
itmuni=1, # time unit 0-undefined, 1-seconds, 2-minutes, 3-hours, 4-days, 5-years
lenuni=2, # length unit 0-undefined, 1-feet, 2-meters, 3-centimeters
nper= nper,
perlen=perlen,
nstp=nstp,
steady=steady,
# tsmult=tsmult
)
# dt1 = perlen[1] * (tsmult[1] - 1) / (tsmult[1]**nstp[1] - 1)
# time_steps = np.array([dt1 * tsmult[1]**i for i in range(nstp[1])])
# cumulative_time = np.cumsum(time_steps)
#Entrada del hidrograma o perturbación al sistema
t0 = perlen[0]
t_delay = 1*24*3600
hydrogram = np.c_[t + t0 + t_delay, np.vstack(Q)] # Hydrogram by rating curve
np.savetxt(os.path.join(path_results, "QHidrom.dat"), hydrogram, fmt="%.4f")
# Save in memory place Modflow
model.add_external(
fname='QHidrom.dat', # name file
unit=70, # unit logic
)
# Dictionary for segment
tabfiles_dict = {
1: { # segment number
'numval': len(hydrogram), # len Qs
'inuit': 70 # Unit logic Hydrogram
}
}
#River definition
river_params_basic = {
'A': 0,
'B': 0.5,
'C': 0,
'start_cell': int(ncol/2),
'orientation': "v",
'width_riv': 1,
'plot': True
}
stage_per_time = [.5]
orientation=river_params_basic["orientation"]
all_river = generate_river_stress_period_data(nrow=nrow, ncol=ncol, topo=topo, nper=nper, stage_per_time=stage_per_time,
del_rbot=5, river_params=river_params_basic, type_gen='basic')
i_riv = all_river["i"]
j_riv = all_river["j"]
river=all_river["river"]
nc_riv = all_river["nc_river"]
# Añadir UPW (Upstream Weighting Package)
upw = flopy.modflow.ModflowUpw(
model,
hk=nlay*[Ka], # Conductividad hidráulica horizontal.(m/d) trabajar entre 5 y 50.
vka=nlay*[Kv_Kh*Ka], # Conductividad hidráulica vertical. (m/d) . trabajar entre 0.1 y 0.3
sy=nlay*[Sy], # Rendimiento específico (acuíferos no confinados) sy de 0,15a 0,35 si es más gravoso
ss=nlay*[Ss], # Almacenamiento específico (acuíferos confinados) MF lo exige incluso en libres. valor que sugieren poner 1e-7 en caso de ac. libre..
laytyp=nlay*[1], # 1 = No confinado, 0 = Confinado. El valor de 1 habilita el drenaje gravitacional.
hdry=-1e30, # Valor de referencia para celdas secas
ipakcb=53, # Unidad de salida para balance de agua
iphdry=0 # Flag para celdas secas (opcional)
)
# Parámetros del solver
nwt = flopy.modflow.ModflowNwt(model, iprnwt=1, maxiterout=20000, fluxtol=1000, headtol=1e-6)
#Reach parametrization
reach_data = flopy.modflow.ModflowSfr2.get_empty_reach_data(nreaches=len(i_riv))
reach_data["k"] = 0 # Capa
reach_data["i"] = i_riv # Filas
reach_data["j"] = j_riv # Columna fija para un río recto
reach_data["iseg"] = 1
reach_data["rchlen"] = delc # Longitud del reach
reach_data["strtop"] = topo[i_riv, j_riv]-d_bkf#topo[i_riv, j_riv]-d_bkf#topo[i_riv, j_riv]-4 # Elevación // ISFROPT > 0
reach_data["ireach"] = list(range(1, nc_riv+1))
reach_data["slope"] = 0 # Pendiente ISFROPT > 0
reach_data["strthick"] = 0 # Espesor del lecho ISFROPT > 0
reach_data["strhc1"] = 0 # Conductividad hidráulica del lecho ISFROPT > 0
#Segment definition
segment_data_steady = flopy.modflow.ModflowSfr2.get_empty_segment_data(1)
segment_data_steady[0]["nseg"] = 1 # Número del segmento
segment_data_steady[0]["icalc"] = 0 # Ecuación de Manning
segment_data_steady[0]["outseg"] = 0 # Segmento terminal
segment_data_steady[0]["flow"] = 0 # Flujo inicial en m³/s
segment_data_steady[0]["elevup"] = topo[i_riv, j_riv][0] - d_bkf # Elevación del lecho del rio aguas arriba
segment_data_steady[0]["elevdn"] = topo[i_riv, j_riv][-1] - d_bkf # Elevación del lecho del rio aguas abajo
segment_data_steady[0]["hcond1"] = Kr # Conductividad hidráulica aguas arriba
segment_data_steady[0]["hcond2"] = Kr # Conductividad hidráulica aguas abajo
segment_data_steady[0]["thickm1"] = dr #1 # Espesor del lecho aguas arriba
segment_data_steady[0]["thickm2"] = dr #1 # Espesor del lecho aguas abajo
segment_data_steady[0]["roughch"] = n_manning #0.1
segment_data_steady[0]["roughbk"] = 0#0.005
segment_data_steady[0]["width1"] = w_bkf#1
segment_data_steady[0]["width2"] = w_bkf#10
segment_data_steady[0]["depth1"] = H_0#1
segment_data_steady[0]["depth2"] = H_0#10
segment_data_transient = flopy.modflow.ModflowSfr2.get_empty_segment_data(1)
segment_data_transient[0]["nseg"] = 1 # Número del segmento
segment_data_transient[0]["icalc"] = 1 # Ecuación de Manning
segment_data_transient[0]["outseg"] = 0 # Segmento terminal
segment_data_transient[0]["flow"] = Q.min() # Flujo inicial en m³/s
segment_data_transient[0]["elevup"] = topo[i_riv, j_riv][0] - d_bkf # Elevación del lecho del rio aguas arriba
segment_data_transient[0]["elevdn"] = topo[i_riv, j_riv][-1] - d_bkf # Elevación del lecho del rio aguas abajo
segment_data_transient[0]["hcond1"] = Kr # Conductividad hidráulica aguas arriba
segment_data_transient[0]["hcond2"] = Kr # Conductividad hidráulica aguas abajo
segment_data_transient[0]["thickm1"] = dr #1 # Espesor del lecho aguas arriba
segment_data_transient[0]["thickm2"] = dr #1 # Espesor del lecho aguas abajo
segment_data_transient[0]["roughch"] = n_manning #0.1
segment_data_transient[0]["roughbk"] = 0#0.005
segment_data_transient[0]["width1"] = w_bkf#1
segment_data_transient[0]["width2"] = w_bkf#10
#SFR parametrization
sfr = flopy.modflow.ModflowSfr2(
model,
nstrm=-reach_data["ireach"].size, # Número de reachs
nsfrpar=0,
nss=segment_data_steady["nseg"].size, # Número de segmentos
const=1.0, # Constante de Manning para flujo en
dleak=1e-4, # Tolerancia de fuga
isfropt=0, # 0 y 1: no flujo vertical NS. P.L. per period. Only when start
ipakcb=53,
istcb2=54,
irtflg =3,
numtim = 50, # Vital para la estabilidad de la solucion
weight=.75, # .5 mas dinamico / 1 procesos mas estables
flwtol=1e-5,
reach_data=reach_data,
segment_data={0: segment_data_steady, 1: segment_data_transient},
transroute=True,
# nstrail=150,
# reachinput=True, # Habilitar solamente cuando se define ISFROPT 1
# itmp =1,
tabfiles=True,
tabfiles_dict=tabfiles_dict,
)
sfr.check()
# sfr.get_outlets()
ibound = np.ones((nlay, nrow, ncol), dtype=int) # Todas las celdas son activas
# ibound[:, 0, :] = 1 # aguas arriba (i=0)
# ibound[:, -1, :] = 1 # aguas abajo (i=nrow-1)
# ibound[:, :, 0] = 0 # borde izquierdo (j=0)
# ibound[:, :, -1] = 0 # borde derecho (j=ncol-1)
# ibound[-1] = 0 # no-flow floor (celda inactiva)
strt = np.ones((nlay, nrow, ncol), dtype=float)*H_0#(topo - d_bkf + H_0)
ghb_data = []
connection_thickness = D # Distancia del cuerpo de agua a la celda con la condicion impuesta
# for k in range(nlay): # recorre cada capa
# for j in range(ncol): # recorre cada columna (dirección transversal al río)
# # Borde superior: aguas arriba
# i_top = 0
# topo_top = topo[i_top, j]
# #h_top = topo_top - d_bkf + H_0 + Slope * L
# h_top = topo_top - d_bkf + H_0
# cond_top = Ka * w_bkf * delc / connection_thickness
# ghb_data.append([k, i_top, j, h_top, cond_top, ])
# # Borde inferior: aguas abajo
# i_bot = nrow - 1
# topo_bot = topo[i_bot, j]
# #h_bot = topo_top - d_bkf + H_0 - Slope * 2 * L
# h_bot = topo_top - d_bkf + H_0 - Slope * L
# cond_bot = Ka * w_bkf * delc / connection_thickness
# ghb_data.append([k, i_bot, j, h_bot, cond_bot])
# ghb = flopy.modflow.ModflowGhb(
# model,
# stress_period_data={0: ghb_data, 1: ghb_data},
# ipakcb=53) # Para salida del presupuesto de flujo
for k in range(nlay):
for j in range(ncol):
if j != int(ncol/2): # Excluir columna del río
# Borde superior
i_top = 0
h_top = topo[i_top, j] - d_bkf + H_0
cond_top = Ka * delr * delc / connection_thickness
ghb_data.append([k, i_top, j, h_top, cond_top])
# Borde inferior
i_bot = nrow - 1
h_bot = topo[i_bot, j] - d_bkf + H_0 - Slope * L
cond_bot = Ka * delr * delc / connection_thickness
ghb_data.append([k, i_bot, j, h_bot, cond_bot])
bas = flopy.modflow.ModflowBas(model, ibound=ibound, strt=strt)
#oc - output control, salidas a guardar
stress_period_data = {}
for per in range(len(perlen)): # n-períodos
for step in range(nstp[per]): # pasos en cada período
stress_period_data[(per, step)] = [
'save head',
'save budget'
]
oc = flopy.modflow.ModflowOc(
model,
stress_period_data=stress_period_data,
compact=True)
# Escribir y correr el modelo
model.write_input()
success, buff = model.run_model(silent=False, report=True)`
Hello,
I have been having some troubles the heads in modflow thus I got this profiles in my river always and I cannot understand why the line is transform to curve. Here I attached the code and image with results
`
import os
import flopy
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.gridspec as gridspec
from utils.inputs import (generate_topo,
generate_river_stress_period_data,
plot_sfr_results, hydrograph_cooper_rorabaugh, plot_results_hds_cbc, extract_sfr_data, calculate_tsmult,
analyze_model_timing)
import phi_numbers2 as phinum
import pandas as pd
from flopy.utils.sfroutputfile import SfrFile
PARÁMETROS DEL MODELO. Everything in SI system (m, s, m3/s, m2, etc)
dt = pd.read_csv('sample_points.csv')
for n in range(50, 55, 5):