Re: [MODFLOW] Transient Head Issue in Hyporheic Zone Simulation

[Luther]

I think this is the issue (see attached). Only the final value of h_0[i, j] is being added to the boundary head array. It needs to be appended at each time step and hopefully your issue will be fixed. 

On Thu, May 30, 2024 at 5:19 PM Nerea Karmele Portillo De Arbeloa <nk.portill...@unitn.it> wrote:

I'm trying to run this code where I'm trying to simulate they hyporheic zone. The head boundary is transient (it changes with each timestep). Here is a replicable example, it is however giving me the same head for all time steps, I've checked everything and cant seem to find the issue. Any insight would be great! Thanks! 

#%% Import Packages
import os
import flopy
import numpy as np
import matplotlib.pyplot as plt
import flopy.utils.binaryfile as bf
from mpl_toolkits.axes_grid1 import make_axes_locatable

model_ws = '...'
os.chdir(model_ws)

# Define the model name and create the MODFLOW model object
modelname = "test"
mf = flopy.modflow.Modflow(modelname, exe_name="mac/mf2005")

Ld = 0.5          # length of 1dune (m)
db = 0.25       # depth of the hyporheic zone (m)
la = 12.56637
g = 9.81          
Hd = 0.05        
slope = 0.0      
Y0 = 0.150        
V = 0.36          
K = 110.0        
Ss = 0.0002      
KH = K  

c = 9.547461516439051e-05
hm = 0.0018358574888470564
a = 12.56637061435917
b = 8.679510469490047e-11

# Define the function for the transient head at the top of the boundary
def transient_head(x, y, t):
    return hm * np.exp(a * y) * np.cos(la * x - la * c * t - ((la * c * Ss) / (2 * a * KH)) * y + (0.35/4))
 

# Model domain and grid definition
Lx, Ly = 0.5, 0.25 # hyporheic zone length in x/ y
ztop, zbot = 0., -0.01 # hyporheic zone top/bottom elev.
nlay, nrow, ncol = 1, 250, 500 # number of layers/rows/cols
delr, delc = Ly/nrow, Lx/ncol  # spacing along rows/cols
delv = (ztop - zbot) / nlay # spacing along depth
botm = -0.01 # bottom elevation  


# Define the time series
nper = 1  # stress periods(day)
perlen = 86400  # length of each stress period
nstp = 6  # number of time steps per stress period
steady = [False] * nper  # all transient periods

# Create the time series for the boundary condition
x = np.linspace(0, Lx, ncol)
y = np.linspace(0, Ly, nrow)
times = np.linspace(0, nper * perlen, nper * nstp)

# Calculate the transient head for each time step
h_0 = np.zeros((len(times), ncol))
for i, t in enumerate(times):
    for j, x_val in enumerate(x):
        h_0[i, j] = transient_head(x_val, 0, t)  # Assuming y=0 for the boundary condition

# Create the discretization object
dis = flopy.modflow.ModflowDis(mf, nlay, nrow, ncol, delr=delr, delc=delc, top=ztop, botm=botm, nper=nper, perlen=perlen, nstp=nstp, steady=steady)

# Variables for the BAS package
# IBOUND
ibound = np.ones((nlay, nrow, ncol), dtype=np.int32)
ibound[:, 0, :] = -1  # constant head at the first row
ibound[:, nrow-1, :] = 0  # no flow at the last row

# Starting heads
strt = np.ones((nlay, nrow, ncol), dtype=np.float32)
strt[:, 0, :] = h_0[0, :]  # initial head at t=0

# Define stress period data for the transient head boundary condition
stress_period_data = {}
for kper in range(nper):
    stress_period_data[kper] = []
    for kstp in range(nstp):
        for j in range(ncol):
            time_index = kper * nstp + kstp
            stress_period_data[kper].append([0, 0, j, h_0[time_index, j], h_0[time_index, j]])  # (layer, row, col, head_start, head_end)

# Create the BAS package with the transient head boundary condition
bas = flopy.modflow.ModflowBas(mf, ibound=ibound, strt=strt)

# Create the CHD package for the transient head boundary condition
chd = flopy.modflow.ModflowChd(mf, stress_period_data=stress_period_data)

# Define function for the homogeneous hydraulic conductivity
hk = K

# Add LPF package to the MODFLOW model
lpf = flopy.modflow.ModflowLpf(mf, hk=hk, vka=1., ipakcb=53)

# Set Output control packages
spd = {}
for kper in range(nper):
    for kstp in range(nstp):
        spd[(kper, kstp)] = ['print head', 'print budget', 'save head', 'save budget']
oc = flopy.modflow.ModflowOc(mf, stress_period_data=spd, compact=True)

# Add PCG package to the MODFLOW model
pcg = flopy.modflow.ModflowPcg(mf, mxiter=1000, iter1=1000, npcond=1, hclose=1e-15, rclose=1e-15, relax=1.0)

# Write the MODFLOW model input files
mf.write_input()

# Run the MODFLOW model
success, buff = mf.run_model()  # success should be True

# Retrieve and plot results
hds = bf.HeadFile(model_ws + '/' + modelname + '.hds')  # heads saved in hds
times = hds.get_times()  # get times

head = hds.get_data(totim=times[1])

# Save budget in modelname.cbc
cbb = bf.CellBudgetFile(model_ws + '/' + modelname + '.cbc')

frf = cbb.get_data(text='FLOW RIGHT FACE', totim=times[-1])[0]  # flow in cols
fff = cbb.get_data(text='FLOW FRONT FACE', totim=times[-1])[0]  # flow in rows

head_1 = head[:, 0:(nrow - 1), 0:(ncol - 2)]

fig = plt.figure(figsize=(3, 3))  # w,h tuple in inches
hds = flopy.utils.binaryfile.HeadFile(model_ws + '/' + modelname + '.hds')
h = hds.get_data(kstpkper=(0, 0))
x = np.linspace(0, Lx, (ncol - 2))
y = np.linspace(0, Ly, (nrow - 1))
y = y[::-1]
fig = plt.figure(figsize=(15, 30))
ax = fig.add_subplot(1, 1, 1, aspect="equal")
c1 = plt.contour(x, y, head_1[0], linestyle='dashed', colors="black")
c = ax.contourf(x, y, head_1[0], cmap='Blues', size=24)
ax.clabel(c1, inline=True, inline_spacing=250, rightside_up=True, fontsize=24, colors="black")
ax.set_xlabel('x [m]', size=28)
ax.set_ylabel('z [m]', size=28)
plt.xticks(size=24)
plt.yticks(size=24)
divider = make_axes_locatable(ax)
cax = divider.append_axes("right", size="7%", pad=0.25)
plt.clabel(c)
cbar = plt.colorbar(c, cax=cax, orientation='vertical')  # colorbar
cbar.ax.tick_params(labelsize=24)
cbar.set_label(label='head [m]', size=28, weight='bold')
plt.show()

--
This group was created in 2004 by Mr. C. P. Kumar, Former Scientist 'G', National Institute of Hydrology, Roorkee. Please visit his webpage at https://www.angelfire.com/nh/cpkumar/
---
You received this message because you are subscribed to the Google Groups "MODFLOW Users Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email to modflow+u...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/modflow/82fa46a9-9d27-478d-b214-2ac2b4aa6d5en%40googlegroups.com.

[Abdul Raheem]

Hi dear community,

I was working on calculating spatio-temporal recharge to groundwater based on remote sensing data.

For that I need soil type data in raster format.

Can anybody help  me with this on how to download it?

Thanks

Raheem

To view this discussion on the web visit https://groups.google.com/d/msgid/modflow/CAGqPyrZj37EaM%2BTgvFv%3D7TEZhoyRY7%2BcP1R9ape3Twq2bjLwRg%40mail.gmail.com.