I’m trying to do run a model using modflow and flopy. I have the steady state model working good, however I would now like to add a transient head profile and see how the model changes. I’m a bit stuck in where to include the stress periods. I want to also do the particle tracking that should be affected by the moving head profile. I’m assuming for the particle tracking I need to change the ReferenceTimeOption to 2 but I do not know how to set the input properly [Period, Step, TimeFraction].

Here is a small replicable example:

import flopy
import numpy as np
import flopy.utils.binaryfile as bf

mf = flopy.modflow.Modflow(modelname, exe_name=mf2005,model_ws=model_ws)

Lx, Ly = 0.5, 0.25 # hz length in x/ y
ztop, zbot = 0., -0.01 # hyporheic zone top/bottom elev.
nlay, nrow, ncol = 1, 100, 100 
delr, delc = Ly/nrow, Lx/ncol  
delv = (ztop - zbot) / nlay 
botm = -0.01 # bot elev. 
kh = 0.0055 

x = np.linspace(0, Lx, ncol) 

def head_f(x, t, v=1.0):
    return np.cos(x/0.025 + v * t)

# Time steps for which heads are computed
time_steps = [0, 5, 10]

# Dictionary to store heads
heads = {}

# Compute and store heads for each time step
for t in time_steps:
    head_values = head_f(x, t)
    heads[t] = head_values

dis = flopy.modflow.ModflowDis(mf, nlay, nrow, ncol, delr=delr, delc=delc, top=ztop, botm=botm)
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

strt = np.ones((nlay, nrow, ncol), dtype=np.float32)
strt[:, 0, :] = heads[0] 
bas = flopy.modflow.ModflowBas(mf, ibound=ibound,strt=strt)

hk_v = np.zeros((nlay, nrow, ncol), dtype=np.float32)   

for j in range(0, ncol, 1):
    hk_v[:, :, j] = kh
    
lpf = flopy.modflow.ModflowLpf(mf, hk=hk_v, vka=1., ipakcb=53)    

spd = {(0, 0): ['print head', 'print budget','save head', 'save budget'],
       (0, 5): ['print head', 'print budget','save head', 'save budget'],
       (0, 10): ['print head', 'print budget','save head', 'save budget']}
oc = flopy.modflow.ModflowOc(mf, stress_period_data=spd,compact=True)

pcg = flopy.modflow.ModflowPcg(mf,mxiter=1000,iter1=1000, npcond=1, hclose=1e-15, rclose=1e-15, relax=1.0)
mf.write_input()
success, buff = mf.run_model() # success should be True
hds = bf.HeadFile(model_ws+'/'+modelname+'.hds') # heads saved in hds
times = hds.get_times() 
head = hds.get_data(totim=times[-1])

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

# Particle tracking

dis_file = '{}.dis'.format(mf.name)
head_file = '{}.hds'.format(mf.name)
bud_file = '{}.cbc'.format(mf.name)

mp = flopy.modpath.Modpath6(modelname="particle", 
                       exe_name='mp6',
                       modflowmodel=mf, 
                       model_ws=model_ws, 
                       dis_file=dis_file,
                       head_file=head_file, 
                       budget_file=bud_file)


hnoflo = -999.99
hdry = -888.88

mpbas = flopy.modpath.Modpath6Bas(mp, hnoflo=hnoflo, hdry=hdry, def_face_ct=0, bud_label= "", def_iface=6, laytyp=0, ibound=ibound, 
                                 prsity=0.38, prsityCB=0.38, extension='mpbas', unitnumber=86)

dis = flopy.modpath.Modpath6Sim(mp, option_flags=[2, 1, 1, 1, 1, 1, 2, 1, 2, 1, 1, 1])
mp.write_input()
mp.run_model()