problem with heads in streamflow routing

[eacuellarq]

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):

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