HZ model using flopy (Modflow/Modpath)

[nportilloda]

I'm trying to do a Modflow/Modpath model for particle tracking through the hyporheic zone, so not the typical groundwater flow model. In particular, I'm trying to represent an unsteady head boundary that changes according to the movement of the bed (dune sediment). I have an analytical solution for the head at the top row (upper boundary) that I'm currently using to update the CHD package at each timestep. This, however, is not giving me consistent results regarding the time related parameters (nper, perlen, nstp). The pathlines of the injected particles change to much from a setup to another. Any insight would be greatly appreciated.

Below is a minimal setup illustrating how I define the model:

times = np.arange(0, nper * perlen, perlen)

  X, Y = np.meshgrid(x, -y)  # Create a 2D grid of x and y values
  times = np.array(times)[:, None, None]  # Reshape times for broadcasting

  h_analytical = headsolver.head_function_transient(X, Y, times)

  ibound = np.ones((nlay, nrow, ncol), dtype=np.int32)
  ibound[:,  0, :] = -1    # constant head at the first row
  ibound[:, nrow-1, :] =  1    # fixed head at bottom row
  
  
  strt = np.ones((nlay, nrow, ncol), dtype=np.float32)
  strt[:, :, :] = h_analytical[0, nrow-1, :]  # initial head at t=0

  steady = [False] * nper
  
  # Discretization
  dis = flopy.modflow.ModflowDis(
      mf, nlay, nrow, ncol,
      delr=delr, delc=delc,
      top=ztop, botm=zbot,
      nper=nper, perlen=perlen, nstp=nstp,
      steady=steady, itmuni=1, lenuni=2
  )

  # Basic package
  bas = flopy.modflow.ModflowBas(mf, ibound=ibound, strt=strt)

  # --- Transient CHD boundaries ---
  stress_period_data = {}
  for kper in range(nper-1):
      stress_period_data[kper] = []
      for j in range(ncol):
          stress_period_data[kper].append([0, 0, j, h_analytical[kper][0,j], h_analytical[kper+1][0,j]])  # (layer, row, col, head_start, head_end)
          
  chd = flopy.modflow.ModflowChd(mf, stress_period_data=stress_period_data)

  # LPF, OC, PCG
  lpf = flopy.modflow.ModflowLpf(
      mf, hk=KH, vka=1.0, sy=Ss, ss=1.0, laytyp=1, ipakcb=53)
  
  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)
  pcg = flopy.modflow.ModflowPcg(
      mf, mxiter=10000, iter1=10000,
      npcond=1, hclose=1e-10, rclose=1e-10, relax=1.0)

  # --- 9. Run & log ---
  mf.write_input()
  success, buff = mf.run_model(silent=False, report=True)
  if not success:
      raise RuntimeError("MODFLOW failed")

[vcantarella]

Have you compared the head results with different time related parameters? My first impression is that due to your analytical solution the problem might be more unstable than usual. I would first make sure that your time-step is small enough so that your heads do not change with different time-steps.